Potential molecular mechanisms of overgrazing-induced dwarfism in sheepgrass (Leymus chinensis) analyzed using proteomic data

Background This study was designed to reveal potential molecular mechanisms of long-term overgrazing-induced dwarfism in sheepgrass (Leymus chinensis). Methods An electrospray ionisation mass spectrometry system was used to generate proteomic data of dwarf sheepgrass from a long-term overgrazed rangeland and normal sheepgrass from a long-term enclosed rangeland. Differentially expressed proteins (DEPs) between dwarf and normal sheepgrass were identified, after which their potential functions and interactions with each other were predicted. The expression of key DEPs was confirmed by high-performance liquid chromatography mass spectrometry (HPLC–MS) using a multiple reaction monitoring method. Results Compared with normal sheepgrass, a total of 51 upregulated and 53 downregulated proteins were identified in dwarf sheepgrass. The amino acids biosynthesis pathway was differentially enriched between the two conditions presenting DEPs, such as SAT5_ARATH and DAPA_MAIZE. The protein–protein interaction (PPI) network revealed a possible interaction between RPOB2_LEPTE, A0A023H9M8_9STRA, ATPB_DIOEL, RBL_AMOTI and DNAK_GRATL. Four modules were also extracted from the PPI network. The HPLC–MS analysis confirmed the upregulation and downregulation of ATPB_DIOEL and DNAK_GRATL, respectively in dwarf samples compared with in the controls. Conclusions The upregulated ATPB_DIOEL and downregulated DNAK_GRATL as well as proteins that interact with them, such as RPOB2_LEPTE, A0A023H9M8_9STRA and RBL_AMOTI, may be associated with the long-term overgrazing-induced dwarfism in sheepgrass. Electronic supplementary material The online version of this article (10.1186/s12870-018-1304-7) contains supplementary material, which is available to authorized users.


Background
Sheepgrass [Leymus chinensis (Trin.) Tzvel.] is a rhizomatous perennial C3 species dominating the Eurasian Steppe grasslands [1]. It adapts well to diverse environmental conditions, such as high alkalinity and salinity, low temperature, drought and various atmospheric N deposition levels [2][3][4][5]. Grazing is the most important economic activity in grassland and is a complex process including touch, defoliation, wounding and bovine serum albumin (BSA) deposition [6]. Moderate herbivory or mowing can stimulate rapid sheepgrass growth [7], whereas long-term overgrazing generally leads to severe decreases in shoot and tiller densities, stem length, plant height and leaf length of sheepgrass, which in turn decreases the aboveground biomass and induces plant dwarfism [8][9][10][11]. Our recent study found that leaf photosynthesis were significantly decreased by grazing of the previous generation, corresponding with the dwarf phenotype of L. chinensis induced by grazing disturbance both in maternal plants in the field and clonal offspring in greenhouse [12]. Grassland productivity is extremely important for proper functioning of the ecosystem and supply of forage for grazing animals [13,14]. Plant dwarfism can lead to significant declines in the aboveground biomass and productivity in grassland, causing changes in the structure and function of the grassland ecosystem through a set of cascade reactions from individual, species and population to ecosystem [15]. Therefore, it is necessary to reveal the mechanisms underlying dwarfism in sheepgrass during grazing.
Major advances have been achieved in previous studies about investigating the molecular mechanisms underlying the response of sheepgrass to grazing. After herbivory, animal saliva significantly increases tiller number, number of buds and biomass of sheepgrass, which are linked to the mobilisation of carbohydrates [7]. Based on RNA sequencing, 2002 genes were identified to be differentially expressed in sheepgrass and to respond to salivary BSA deposition during grazing, indicating that grazing affects plant recovery probably via salivary BSA [16]. Thousands of genes were also identified to be differentially expressed in L. chinensis after wounding and defoliation [17]. However, previous studies mainly focussed on the effects of stress caused by short-term or instantaneous grazing on sheepgrass growth at genome level, and the effect of stress caused by long-term overgrazing, particularly at the protein level, has not been studied extensively.
It is well known that proteins are vital parts of living organisms, with many functions. The term proteomics was coined in 1997 by James [18] in analogy with genomics. Proteomics is more complicated than genomics because the genome of organism is more or less constant, whereas the proteome differs from cell to cell and from time to time. Therefore, in this study we used proteomic profiling to compare the proteins expression in dwarf sheepgrass from a rangeland with a long history of overgrazing and that in normal sheepgrass from a rangeland with a long history of enclosure. Differentially expressed proteins (DEPs) between dwarf and normal sheepgrass were identified and their potential functions and interactions with each other were analysed. In addition, the expression levels of key DEPs were confirmed using high-performance liquid chromatography mass spectrometry (HPLC-MS) using the multiple reaction monitoring (MRM) method. The results may extend our understanding of proteomic changes in sheepgrass in response to long-term overgrazing and the molecular mechanisms underlying plant dwarfism.

Plant materials
Sheepgrass plants were sampled from a long-term overgrazed rangeland (n = 3, GZ group) and from an adjacent long-term enclosed rangeland that had been enclosed since 1983 for long-term ecological observation and research (n = 3, NG group). The two sampling sites were in fact in the same area (distance < 10 m) and were separated by a pasture fence. Photographs of the sheepgrass plants in the GZ and NG groups are depicted in Fig. 1  Sciences. Additionally, the collection of plant materials has been permitted by the landowners.
After trimming the aboveground portion of sheepgrass plants to 2 cm, rhizome samples of sheepgrass were grown in an incubator at 25°C under long-day conditions (16 h light and 8 h dark cycle) to form asexual buds. In this step, modified Hoagland's medium (pH 6.0, adjusted with NaOH) was utilised, consisting of Ca(NO 3  Subsequently, asexual buds were separated from the parent body and grown under the same conditions for 50 days to obtain young seedlings. The seedlings (three biological replicates for each group) were then washed with distilled water; samples were collected from leaves using sterilised scissors and tweezers. The samples were immediately frozen in liquid nitrogen and stored at − 80°C for further analysis.

Protein preparation and analysis
The samples were grinded into a powder in liquid nitrogen, and then the powder was immersed in 30 mL acetone containing 10% (w/v) trichloroacetic acid overnight. Following refrigerated high-speed centrifugation, the sediment was washed three times with acetone and dried. SDT-generated lysates were ultrasonicated and centrifuged. The protein level in the supernatant was quantified using the microbicinchoninic acid method (Pierce) [19]. Then, the proteins were separated using 12.5% sodium dodecyl sulphate polyacrylamide gel electrophoresis. Proteins were hydrolysed into peptides using trypsin [20], and the peptides were quantified at 280 nm. Each peptide sample (70 μg) was labeled using the iTRAQ Reagent-8plex Multiplex Kit (AB SCIEX, Foster City, CA, USA). The labelled peptides were separated on a reversed-phase C18 column (75 μm × 250 mm, 3 μm; Column Technology Inc., Fremont, CA, USA) in a capillary high-performance liquid chromatograph (EASY nLC1000; Proxeon, Denmark) with a linear gradient of 0-55% mobile phase B (0.1% formic acid-84% acetonitrile) in mobile phase A (0.1% formic acid aqueous solution) from 220 to 228 min, followed by a linear gradient of 55-100% mobile phase B from 229 to 240 min. Finally, the obtained fractions were analyzed on a Q-Exactive mass spectrometer (MS) (Ther-moFinnigan, CA, USA) (full scan range: 300-1800 m/ z, detection mode: positive ion, MS1 resolution: 70,000 at 200 m/z, automatic gain control target: 3e6, maximum IT: 20 ms, dynamic exclusion: 25.0 s, and number of scan ranges: 1). Ten fragment maps were collected in each full scan for the MS2 scan (activation type: higher energy collisional dissociation, isolation window: 2 m/z, resolution: 17,500 at 200 m/z, microscans: 1, maximum IT: 60 ms, normalised collision energy: 29 eV, and underfill ratio: 0.1%). The three biological replicates in two groups were analyzed independently.
The proteomic data are available at the EMBL-EBI Proteomics Identifications (PRIDE) database with the accession number PXD006548.

Data preprocessing
Raw proteomic data were converted into mzXML data using ReAdW software (http://tools.proteomecenter.org/wiki/index.php?title=Software:ReAdW). mzXML data was then matched with the Swiss-Prot/ TrEMBL database (release-2015_07) (http://www.uniprot.org/uniprot) [21,22] using the PeptideProphet tool in Trans-Proteomic Pipeline software (http:// www.proteomecenter.org) [23]. The matching results were obtain using pepXMLTab software. The following search parameter settings were used: peptide tolerance, ± 20 ppm and tandem mass spectrometry tolerance, 0.1 Da. A peptide was required to have at least a single assigned fragment and only unique peptides were used for protein identification. Protein with a false discovery rate (FDR) < 0.01 were considered to have high credibility.

Functional annotation and subcellular localisation of proteins
Functional annotation was performed for the proteins with high credibility based on the Gene Ontology (GO) database (http://www.geneontology.org/) [24] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database (http://www.kegg.jp/kegg/pathway. html) [25] using the enrichment analysis tool Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 [26]. The subcellular localisation of the proteins was determined through Cell-PLoc 2.0, a package of web servers for predicting the subcellular localisation of proteins in various organisms [27].

DEPs identification and functional analysis
The t-test method in limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html) was utilised to identify DEPs between GZ and NG groups. The p-value for each protein was adjusted by the Benjamini-Hochberg method [28]. Only proteins meeting the cut-off criteria of fold change > 1.2 and adjusted p-value < 0.05 were identified as DEPs. The DEPs with similar expression patterns were clustered with pheatmap software [29].
To reveal potential functions of DEPs, they were subjected to functional enrichment analysis based on the GO database, and only GO terms with p-value < 0. 05 were considered significant. The pathway location of DEPs was determined using KEGG Mappertool (http:// www.kegg.jp/kegg/mapper.html).

Construction of protein-protein interaction (PPI) network for DEPs
Protein sequences of DEPs in fasta format were downloaded from the UniProt database (http://www. uniprot.org/) [30] and then homologous proteins of these DEPs in Arabidopsis thaliana were entered into the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://string-db. org/) [31], which contains known and predicted protein interactions. PPIs obtained from STRING with a combined score > 0.4 were used to construct a PPI network that was visualised using Cytoscape (http:// www.cytoscape.org/) [32]. In addition, modules from the PPI network were screened using the plug-in MCODE in Cytoscape to identify potential key protein networks.

Determination of the DEP expression using HPLC-MS
Labelled protein powder samples in the GZ (n = 3) and NG groups (n = 3) were dissolved in solution A (H 2 O:acetonitrile:formic acid = 98:2:0.1). From each sample, 1 μL was absorbed and pooled. The pooled sample was centrifuged and the supernatant was enriched in the C18 column (100 μm i.d. × 20 mm, 5 μm) and then graded in the C18 column (150 μm i. d. × 100 mm, 1.9 μm) using a gradient elution pattern. The mobile phase A was 2% acetonitrile with 98% H 2 O and 0.1% formic acid. Mobile phase B was 80% acetonitrile with 20% H 2 O and 0.1% formic acid. The flow rate was set at 0.6 μL/min.
After the gradient elution, total six samples in two groups were entered into a Q ExactiveHPLC system (Thermo Fisher Scientific, Hudson, NH, USA) equipped with an NCS3500 MS system (Dionex, Sunnyvale, CA, USA). Full MS scan with MRM mode used the following settings: scan range = 300-1400 Da, scan time = 90 min, spray voltage = 2.20 KV, capillary temperature = 320°C, normalised collision energy = 27%, first-order scanning resolution = 120,000, AGC = 3e6 and maximum IT = 80 ms. The top 20 ions were chosen for the second scan in accordance with the following conditions: second-order scanning resolution = 150,000, AGC = 5e4 and maximum IT = 45 ms.
Finally, the generated data were analyzed using Skyline software [33] and a reference atlas database was created using the MASCOT (version: 2.1.0) [34] protein identification platform (Matrix Science, London, UK). The peptide fragments with idotp and dotp > 0.9 were considered to be credible.

Protein identification and function annotation
After matching MS/MS data with protein databases, 23,387 peptides were identified, among which approximately 80% had 13-27 amino acids. A total of 6555 proteins were found from six samples, 1022 of which had high credibility (FDR < 0.01) (Additional file 1).
The functional annotation of these 1022 proteins revealed that most of them were cytoplasmic proteins and nucleoproteins related to carbohydrate metabolism (Fig. 2a). According to the subcellular localisation analysis, 147 proteins were located in the plastid, 76 in the nucleus and 48 in the cytoplasm (Fig. 2b).

Functional analysis of DEPs
In total, 104 proteins were differentially expressed between GZ and NG groups, with 51 being upregulated and 53 downregulated in the GZ group (Additional file 1). GO enrichment and KEGG pathway location analyses revealed that the upregulated DEPs were particularly associated with several GO terms related to metabolic pathways, such as pseudouridine synthesis (B8BQ07_THAPS), polyamine catabolic process (PAO3_ARATH) and haem biosynthetic process (J3LTE1_ORYBR), (Table 1). In addition, the downregulated DEPs were markedly associated with functions such as reductive pentose-phosphate cycle (RBL_AMOTI and A0A023HP98_9POAL), cytoplasmic translation (RS72_ARATH) and DNA replicationrelated DNA unwinding (GYRA_ARATH) ( Table 2). The pathway location analysis mainly localised DEPs in metabolic pathways (e.g., SAT5_ARATH, PME21_ ARATH and DAPA_MAIZE), such as, biosynthesis of secondary metabolites (SAT5_ARATH, DAPA_MAIZE and CAS1_ARATH), microbial metabolism in diverse environments (SAT5_ARATH, DAPA_MAIZE and DPNP1_ARATH) and biosynthesis of amino acids (SAT5_ARATH and DAPA_MAIZE) ( Table 3). All results of GO and pathway enrichment analyses are shown in Additional file 1.

PPI network
A total of 89 proteins in A. thaliana were homologous to DEPs and among them, 30 proteins were predicted to interact with each other (Fig. 3). The PPI network revealed a possible interactions between the proteins of A0A023H9M8_9STRA, RPOB2_LEPTE, ATPB_DIOEL, DNAK_GRATL and RBL_AMOTI.

Validation of DEPs expression
The HPLC-MS system with the MRM mode was used to verify proteomic data at the protein level. Information on the target peptides selected for validation is listed in Table 4. In accordance with the criteria of idotp and dotp > 0.9, three peptide fragments were credible, including two fragments of ATPB_ DIOEL and one of DNAK_GRATL. ATPB_DIOEL were confirmed to be upregulated and DNAK_ GRATL was confirmed to be downregulated in the GZ group (Table 5). These findings were consistent with the results of proteomic data in the present study.

Discussion
In this study, a total of 1022 proteins with high credibility were identified from dwarf sheepgrass from the longterm overgrazed rangeland and normal sheepgrass from the long-term enclosed rangeland. Among them, 51 upregulated and 53 downregulated proteins were identified in the dwarf samples. The PPI network revealed a possible interaction between the proteins RPOB2_LEPTE, A0A023H9M8_9STRA, ATPB_DIOEL, RBL_AMOTI and DNAK_GRATL. The HPLC-MS analysis confirmed that ATPB_DIOEL was upregulated and DNAK_GRATL was downregulated in the dwarf samples. ATPB_DIOEL is a ATP synthase subunit beta in the chloroplast. The ATP synthase complex catalyses ATP synthesis in photosynthesis, a process termed as photosynthetic phosphorylation [35]. A previous study has revealed that the activation state of ATP synthase can limit leaf-level photosynthesis [36]. In this study, the expression level of ATPB_DIOEL was confirmed to be upregulated in the dwarf sheepgrass from the long-term overgrazed rangeland. This indicated that long-term overgrazing may upregulate the expression of ATPB_ DIOEL, promoting the formation of ATP synthase and thereby, limiting the photosynthesis of sheepgrass and restricting plant growth.
The expression level of DNAK_GRATL, one of the proteins predicted to interact with ATPB_DIOEL, was confirmed to be downregulated in dwarf sheepgrass, and it was heat shock protein 70 (hsp70). hsp70s can assist in a wide range of protein-folding processes in almost all cellular compartments. hsp70s have critical functions in preventing aggregation and assisting protein refolding  under normal and stress conditions [37]. hsp70s are essential for plant development and hsp70 mutant plants exhibit growth retardation [38]. Therefore, in this study, decreased expression of DNAK_GRATL may be associated with long-term overgrazing-induced dwarfism in sheepgrass. Downregulated expression of RBL_AMOTI, another protein that was predicted to interact with ATPB_DIOEL, was predicted to be associated with reductive pentosephosphate cycle. RBL_AMOTI is the large chain of ribulose bisphosphate carboxylase (RubisCO), which is also a chloroplast enzyme and participates in the reductive pentose-phosphate cycle [39]. RubisCO catalyses the carboxylation of D-ribulose-1,5-bisphosphate, which is the primary event in both carbon dioxide fixation and pentose substrate oxidative fragmentation in the photorespiration process [40,41]. Therefore, RubisCO can support photosynthesis and plant growth [42,43]. In this study, the expression level of RBL_AMOTI was reduced in dwarf sheepgrass from the long-term overgrazed rangeland. Although there is no other evidence to prove the association of RubisCO with overgrazing or plant dwarfism, we speculated that the expression level of RBL_AMOTI was decreased because of long-term overgrazing, thereby reducing RubisCO synthesis and then limiting photosynthesis and slowing down sheepgrass growth.
In addition, RPOB2_LEPTE and A0A023H9M8_9STRA interacted with ATPB_DIOEL as well as DNAK_GRATL and RBL_AMOTI as shown in the PPI network. RPOB2_ LEPTE is the subunit beta C-terminal section of DNAdirected RNA polymerase that catalyses the transcription of DNA into RNA [44]. A0A023H9M8_9STRA is encoded by the chloroplast gene rpoC2, which also encodes the beta subunit of RNA polymerase [45]. Currently, there is no evidence to support the associations of these two proteins with plant growth or overgrazing. We speculated that RPOB2_LEPTE and A0A023H9M8_9STRA are involved in long-term overgrazing-induced dwarfism in sheepgrass through their interactions with other proteins, such as ATPB_DIOEL, DNAK_GRATL and RBL_AMOTI.
Furthermore, the present study showed that metabolic pathways such as the biosynthesis of secondary metabolites Fig. 3 The protein-protein interaction network of the differentially expressed proteins. The orange nodes represent the upregulated proteins in dwarf sheepgrass, the green nodes represent the downregulated proteins and the purple nodes represent the proteins predicted to interact with the differentially expressed proteins and the biosynthesis of amino acids were particularly associated with a series of DEPs, such as SAT5_ARATH and DAPA_MAIZE if applicable. SAT5_ARATH (serine acetyltransferase 5) is a key enzyme in cysteine biosynthesis during sulphur assimilation in higher plants [46]. Sulphur is required for the growth of all organisms, and inorganic sulphate in soil can be assimilated by plants to synthesise sulphur-containing amino acids, such as cysteine and methionine, to support plant growth [47]. However, we found that the expression of SAT5_ARATH was downregulated in dwarf sheepgrass. We speculated that long-term overgrazing decreased the expression of SAT5_ARATH, thereby inhibiting sheepgrass growth. DAPA_MAIZE [4hydroxy-tetrahydrodipicolinate (HTPA) synthase] is a homotetrameric enzyme of lysine biosynthesis that catalyses the condensation of (S)-aspartate-beta-semialdehyde [(S)-ASA] and pyruvate to HTPA [48]. A lysine-rich arabinogalactan protein AtAGP19 functions in multiple processes during plant growth and development, including cell division and expansion, leaf development and reproduction [49]. Protein lysine acetylation plays regulatory roles in the photosynthesis process and Calvin cycle in plants [50][51][52].
These results indicate the critical role of lysine in plant growth. In the current study, the expression level of DAPA_MAIZE was decreased in dwarf sheepgrass. Therefore, long-term overgrazing reduced the expression level of DAPA_MAIZE in sheepgrass, which may suppress photosynthesis, leading to the exhibition of dwarfism.
A previous study of Huang et al. [16] employed nextgeneration sequencing technology to characterize de novo the transcriptome of sheepgrass after defoliation and grazing treatments and to identify differentially expressed genes responding to grazing and BSA deposition. Enrichment analysis of 2002 differentially expressed genes revealed that the effects of grazing and BSA deposition involved cell oxidative changes and apoptosis. However, our present study did not enrich functions associated with cell oxidative changes and apoptosis. Instead, we obtianed several proteins that may be related to photosynthesis, which was in consistent with our recent study as we mentioned above [12].
Despite the aforementioned results, a main limitation of the present sutdy is that the protein interaction and their functions were not confirmed by experiments. Therefore, in future studies, we will experimentally confirm the expression of the proteins (e.g. RPOB2_LEPTE, A0A023H9M8_9STRA, ATPB_DIOEL, RBL_AMOTI and DNAK_GRATL) and their interactions in dwarf  sheepgrass as well as the functions of SAT5_ARATH and DAPA_MAIZE. We plan to further investigate the associations of these proteins with dwarfism in sheepgrass.

Conclusions
In conclusion, based on ESI-MS data, 104 proteins were identified to be differentially expressed between dwarf sheepgrass from the long-term overgrazed rangeland and normal sheepgrass from the long-term enclosed rangeland. HPLC-MS analysis confirmed that the expression levels of ATPB_DIOEL were upregulated in dwarf sheepgrass and those of DNAK_GRATL were downregulated. These two proteins and the proteins with which they interact, as shown in the PPI network, such as RPOB2_LEPTE, A0A023H9M8_9STRA and RBL_AMOTI, may be associated with the long-term overgrazing-induced dwarfism in sheepgrass. The downregulated expression levels of SAT5_ ARATH and DAPA_MAIZE may also play key roles in this dwarfism probably via amino acid synthesis. These results provide novel information for further experimental studies and contribute to a better understanding of the molecular mechanisms underlying dwarfism in sheepgrass.

Additional file
Additional file 1: The identified 1022 proteins that had high credibility (FDR < 0.01); 104 differentially expressed proteins between long-term overgrazed rangeland (GZ) and adjacent long-term enclosed rangeland Availability of data and materials All data generated or analyzed during this study are included in this published article [and its supplementary information files].
Authors' contributions WR and JZ participated in the design of this study. JX and LK performed the statistical analysis. XH carried out the study, together with XL, and collected important background information. ZW drafted the manuscript. HG, NH and CC conceived of this study, and participated in the design and helped to draft the manuscript. All authors read and approved the final manuscript.
Ethics approval and consent to participate This article does not contain any studies with human participants or animals performed by any of the authors. The ethics approval is unnecessary for our study.

Competing interests
The authors declare that they have no competing interests.

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.