Skip to main content
  • Research article
  • Open access
  • Published:

Bovine serum albumin in saliva mediates grazing response in Leymus chinensis revealed by RNA sequencing



Sheepgrass (Leymus chinensis) is an important perennial forage grass across the Eurasian Steppe and is adaptable to various environmental conditions, but little is known about its molecular mechanism responding to grazing and BSA deposition. Because it has a large genome, RNA sequencing is expensive and impractical except for the next-generation sequencing (NGS) technology.


In this study, NGS technology was employed to characterize de novo the transcriptome of sheepgrass after defoliation and grazing treatments and to identify differentially expressed genes (DEGs) responding to grazing and BSA deposition. We assembled more than 47 M high-quality reads into 120,426 contigs from seven sequenced libraries. Based on the assembled transcriptome, we detected 2,002 DEGs responding to BSA deposition during grazing. Enrichment analysis of Gene ontology (GO), EuKaryotic Orthologous Groups (KOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways revealed that the effects of grazing and BSA deposition involved more apoptosis and cell oxidative changes compared to defoliation. Analysis of DNA fragments, cell oxidative factors and the lengths of leaf scars after grazing provided physiological and morphological evidence that BSA deposition during grazing alters the oxidative and apoptotic status of cells.


This research greatly enriches sheepgrass transcriptome resources and grazing-stress-related genes, helping us to better understand the molecular mechanism of grazing in sheepgrass. The grazing-stress-related genes and pathways will be a valuable resource for further gene-phenotype studies.


Herbivore feeding is a complex process that includes wounding, defoliation, and BSA deposition [1]. Often, the leaves of a plant are completely or partially removed, affecting the photosynthetic activity, secondary metabolism, and carbohydrate relocation of plants [25]. Many reports have focused on stress-induced gene expression, photosynthetic capacity, root growth, and nutrient uptake after wounding and defoliation [69]. Resistance to herbivores depends largely upon aboveground portion of plants and on leaf-to-leaf wound signaling, which may involve electrical signaling [10]. The propagation of electrical activity leads to the expression of defense genes not only in wounded leaves but also throughout aboveground portion of plants.

Animal saliva deposited on plants (especially on leaves) is another important factor affecting plant recovery after grazing. In 1960, Vittora and Rendina first proposed that interactions between grazers and plants involve the deposition of herbivore saliva during grazing. This hypothesis has since been tested [1114]. To date, many studies of large-herbivore BSA deposition have focused on macroscopic changes in plants, such as biomass accumulation, tiller, and increased bud initiation after grazing [15], rather than on changes in gene expression and plant physiology. On the other hand, many growth regulators in insect salivary systems have been fully researched, such as glucose oxidase, β-glucosidase [16, 17], and various growth factors in certain mammalian submaxillary glands (mainly mouse and human), including thiamine [14], nerve growth factor (NGF), transforming growth factor (TGF), and epidermal growth factor (EGF) [18], are known to have growth-regulating activity. Growth factors can intervene directly in cellular metabolism by promoting differential gene expression and are expected to be active in a variety of organisms [19]. Injecting growth-promoting substances from grasshopper saliva into Bouteloua gracilis stimulates tiller production [11]. Mouse and human EGF can enhance plant growth rates and promote cell division in the epicotyl [20]. However, to date, no study has reported the effects of BSA deposition by large herbivores such as cows, sheep, and camels.

Gene-expression profiling or transcriptome analysis can provide new insights to understand the molecular mechanism of grazing responses in plants. High-throughput next-generation sequencing (NGS) technologies, such as 454 (ROCHE), Solexa (Illumina), and SOLiD (ABI), have been widely and effectively used to generate large-scale transcriptome data in many plant species [2128], including sheepgrass (Leymus chinensis) [29, 30]. Sheepgrass is an important forage species in the genus of Leymus, with good quality, high nutrition value, and various stress resistance [3134]. Its genomic formula was NsNsXmXm and Elymus californicus should be the maternal donor transferred from the genus Elymus to Leymus[35]. One copy of the haploid genome of sheepgrass contains 9.65-Gb, and high-throughput NGS technologies make it possible to generate genome resources at relatively low cost. So far, sheepgrass transcriptome databases have been generated under saline-alkaline treatment [29] and freezing treatment [30] using Roche-454 massive pyrosequencing technology. These databases provide numerous DEGs for two stresses. Recently, a comparative transcriptomics analysis of the Illumina sequencing data was conducted, and the results revealed common and distinct mechanisms for sheepgrass responses to defoliation compared to mechanical wounding [36]. Based these transcriptome databases, some grazing responsive genes were cloned and identified, such as LcSUT1 and LcDREB3[37, 38].

Here, we focus on profiling the effects of herbivore saliva on sheepgrass and distinguishing BSA deposition from defoliation in grazing. In our study, we use bovine serum albumin (BSA) instead of bovine saliva to perform the grazing simulation treatment. The components of herbivore saliva are unstable and it usually contains bacterium. BSA is an important protein in bovine saliva. Its homolog has been found in ovine saliva and probably has interactions with plants [39]. In our study, in order to enrich sheepgrass transcriptome resource, accelerate our understanding of the genetic basis of grazing stress, we used Illumina GAIIx technology to sequence sheepgrass transcriptome after defoliation and grazing. We compared defoliation and grazing treatments, and identified the differentially expressed genes (DEGs) responding to BSA deposition and corresponding pathways involved in saliva effects. We performed further biochemical and morphological experiments to verify these results in transcriptome.


Illumina sequencing and Trinity transcriptome assembly of sheepgrass

We obtained 47,782,901 raw reads from the seven libraries representing different time points or treatments via Illumina GAIIx sequencing (Table 1). The sequence reads generated in this study were deposited in the NCBI sequence-read archive (SRA065691). The raw reads were filtered by a stringent criterion as described in Methods. The remaining reads were considered clean. An additional file show this in more detail (see Additional file 1).

Table 1 Statistics summary of Illumina sequencing data generated for sheepgrass transcriptome

We obtained 120,426 contigs (≥200 bp) using the Trinity assembly software. The mean contig size was 634 bp, and the contig size ranged from 201 to 28,343 bp. About one-third of the contigs were longer than 500 bp, and 20,816 contigs were longer than 1,000 bp. Additional file 2 show the quality of the assembly transcripts in more detail (see Additional file 2).

Contig assembly and gene overview

Using the Trinity de novo assembly program, a total of 120,426 contigs were obtained. 79,459 contigs were detected in Library C (the control). 83,189, 85,184 and 77,786 contigs were detected in three defoliation libraries (D2, D6, and D24), respectively. Excluding these repeated contigs in the three samples, there were 110,955 contigs detected in defoliation libraries. 69,112, 80,829 and 55,874 contigs were detected in three grazing libraries (G2, G6, and G24), respectively, with a total of 99,026 in grazing libraries. The three treatments are summarized in Figure 1. The quality of the contigs in the seven samples is shown in Additional files 3 & 4.

Figure 1
figure 1

The venn diagram of gene counts in the control, defoliation and grazing treatments. The control contains Library C. The defoliation treatment contains Library D2, D6, and D24. The grazing treatment contains Library G2, G6, and G24. Numbers in parentheses are total number of expressed genes in each treatment.

Functional annotation and descriptive profile

Gene ontology (GO) assignments were used to classify the functions of the predicted sheepgrass genes expressed in response to grazing stress. Based on sequence homology, 9,831 genes were assigned at least one GO term, including 49 second-level functional categories (Figure 2). An additional docx file show the summary of WEGO output data in more detail (see Additional file 5). Among the assigned terms, “cell” (7,508 terms, 76.4%), “cell part” (7,508 terms, 76.4%), “organelle” (4,311 terms, 43.9%) and “organelle part” (1,643 terms, 16.7%) were dominant in the cellular component. “Cellular process” (6,690 terms, 68.1%), “metabolic process” (6,378 terms, 64.9%), “biological regulation” (2,011 terms, 20.5%), “pigmentation” (1,902 terms, 19.3%), and “response to stimulus” (1,701 terms, 17.3%) were dominant among biological processes. The absolute majority of molecular-function terms were clustered in “binding” (6,850 terms, 69.7%) and “catalytic activity” (5,927 terms, 60.3%).

Figure 2
figure 2

GO classifications of assembled transcripts using WEGO software ( The genes were assigned to three main categories: biological process, molecular function and cellular component. The right hand y-axis indicates the number of annotated genes. The left hand y-axis indicates the percentage of annotated genes.

To further evaluate the completeness of the de novo transcriptome assembly and to predict the gene functions, all assembled transcripts were compared against the EuKaryotic Orthologous Groups (KOG) database. This comparison revealed 9,985 sequences with significant homology, each of which was assigned to the appropriate KOG cluster. These KOG classifications were grouped into 25 functional categories (Figure 3). The five largest categories were “signal-transduction mechanisms” (16.64%), “general function prediction only” (9.87%), “posttranslational modification, protein turnover, chaperones” (9.29%), “translation, ribosomal structure, and biogenesis” (5.34%), and “intracellular trafficking, secretion and catabolism” (5.14%).

Figure 3
figure 3

KOG function classifications of assembled transcripts. The contigs were assigned to the KOG database to predict possible functions. A total of 9,985 contigs were assigned to 25 categories.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for the systematic understanding of high-level gene functions in terms of biological networks, such as the cell, organism, and ecosystem, from molecular-level information ( The assembled transcripts were searched against the KEGG database using BLASTX with a cut-off E-value of 10−5 to identify the biological pathways related to grazing responses in sheepgrass. We obtained 6,820 matching terms, which were assigned to 275 KEGG pathways in 5 main biological processes. The major pathways were “biosynthesis of amino acids” (ko01230, 191 transcripts), “ribosome” (ko03010, 165 transcripts), “carbon metabolism” (ko01200, 163 transcripts), “purine metabolism” (ko00230, 137 transcripts), “spliceosome” (ko03040, 119 transcripts), “protein processing in endoplasmic reticulum” (ko04141, 118 transcripts), and “RNA transport” (ko03013, 114 transcripts).

We applied gene enrichment analysis between the grazing libraries (G2, G6, and G24) and the control (C) in Table 2. The results showed 10 GO second-level functional categories responded to grazing treatment. These groups were “transporter activity”, “oxidoreductase activity”, “catalytic activity” and “lyase activity” of molecular function, “amino acid and derivative metabolism”, “transport” and other three groups of biological process, “external encapsulating structure” of cellular component. In KOG functional classification, “Lipid transport and metabolism”, “Amino acid transport and metabolism” and “Energy production and conversion” were significantly different (q-value < 0.01). 10 KEGG pathways showed significant difference (q-value < 0.01) between the grazing treatment and the control. Pathways “Two-component system” and “ABC transporters” were two of the most different ones.

Table 2 Gene enrichment analysis of DEGs from grazing vs. control

Differentially expressed genes (DEGs) in response to BSA deposition

To investigate the changes in gene expression and understand the critical genes involved in the response of sheepgrass to BSA deposition, we collected the RPKM means from the grazing (G2, G6, and G24) and defoliation libraries (D2, D6, and D24) and analyzed the significant difference genes using DEGseq package [40]. The differentially expressed genes (DEGs) in response to defoliation (D), grazing (G) and BSA deposition (G vs. D) were determined from the changes in expression of defoliation vs. control, grazing vs. control, and grazing vs. defoliation. Using a screening criterion of FDR < 0.001, we found 2,002 genes that responded to BSA deposition, of which 365 were induced and 1,637 were inhibited. We found 3,156 genes that responded to grazing and 2,759 genes that responded to defoliation (Table 3). Figure 4 show the details of expression changes of DEGs in MA plot. A cluster analysis of the libraries of control, defoliation, and grazing was performed using the heat map shown in Additional file 6.

Table 3 Gene statistics of the differentially expressed genes
Figure 4
figure 4

The expression change of DEGs responded to defoliation, grazing, and BSA deposition in MA plots. The MA plot uses M as the y-axis (log2 fold-change) and A as the x-axis (log2 RPKM mean). The red dots stand for the DEGs which are up-regulated, and the green dots stand for those down-regulated DEGs.

GO functional-enrichment analysis was performed for 377 GO terms from BSA deposition (G vs. D) DEGs compared to the 9,831 GO terms from the full transcriptome at a Bonferroni-corrected P-value ≤ 0.05. The results showed the biological-process categories “cell death” and “transport” and the molecular-function categories “antioxidant activity”, “oxidoreductase activity”, and “transporter activity” were significantly enriched (Table 4). In the KOG enrichment analysis, 185 KOG annotation terms were assigned to the saliva-deposition DEGs. “Energy production and conversion”, “amino acid transport and metabolism” and “lipid transport and metabolism” were the significant enrichment categories (Figure 5, Table 5).

Table 4 Enriched GO categories of DEGs from BSA deposition (G vs. D)
Figure 5
figure 5

KOG function classifications of DEGs responding to BSA deposition. The DEGs were aligned to the KOG database to predict possible functions related to BSA deposition. A total of 185 DEGs were assigned to 19 categories.

Table 5 Enriched KOG categories of DEGs from BSA deposition (G vs. D)

In the KEGG pathway enrichment analysis, 307 annotation terms were obtained for the saliva-deposition DEGs. These 307 terms belonged to 136 pathways. We found that 11 pathways were statistically enriched in the DEGs responding to BSA deposition. As shown in Table 6, the most enriched pathways were related to environmental stress responses. For example, two-component systems play important roles in signal transduction in response to environmental stimuli and growth regulators, light and osmotic stress.

Table 6 Enriched pathways of DEGs from BSA deposition (G vs. D)

By aligning our reads to the reference transcriptome dataset obtained by Roche-454 massive pyrosequencing technology [30], we obtained further details and annotations (Table 1). Among the DEGs responding to BSA deposition were serine/threonine protein kinase, apoptotic ATPase, and aldehyde dehydrogenase-family proteins. Most of these genes showed different modes of expression 24 h after defoliation or grazing, as shown in Table 7.

Table 7 Cell oxidative and apoptosis related genes in DEGs

Programmed cell death after clipping and the effect of BSA

To confirm the reliability of apoptosis in response to BSA deposition during grazing, DNA fragments from cells from the cut ends of leaves treated with water or BSA were analyzed using agarose-gel electrophoresis. As shown in Figure 6, most of the DNA fragments appeared on the third and fourth days after clipping in both the BSA- and water-treated groups. Genomic DNA declined by day 5 and disappeared by day 7 in the water-treated group. In the BSA-treated group, however, genomic DNA decreased slowly and was still observed on day 9.

Figure 6
figure 6

Agarose gel images of the DNA fragments. M refers to the DNA marker 2000, C refers to the control, W refers to the water-treated leaf samples, and B refers to the BSA-treated leaf samples. Lanes from left to right: control, water-treated and BSA-treated leaf samples from days 1 to 9 post-treatment.

Accumulation of oxidative-stress-related factors in grazed sheepgrass

The H2O2 levels increased significantly in the water- and ovalbumin (OVA)-treated groups compared to the unclipped controls (p < 0.01). The OVA treatment was used as a control for BSA to eliminate the interference of proteins daubed on the cut surface. Daubing with BSA marginally affected the H2O2 levels in the cells, but this difference was not significant compared to the controls (Figure 7).

Figure 7
figure 7

The concentration analysis of H2O2, MDA, and SOD in the leaf scar cells. The statistical data were treated by ANOVA (Analysis of Variance) test and S-N-K (Student-Newman-Keuls) test. “H2O2” and “MDA” groups were significant while “SOD” group was not significant by ANOVA test. The same letter over two treatments stands for no significant difference and two treatments under different letters were significantly different in S-N-K test.

The changes in malondialdehyde (MDA) were similar to those in hydrogen dioxide (H2O2). The MDA levels in the water-treated and OVA-treated group were approximately 3- to 4-fold higher than those in the unclipped and BSA-treated groups, and this difference was significant at p < 0.01. The MDA levels in the BSA-treated cells increased slightly compared to the controls (Figure 7).The superoxide dismutase (SOD) levels in the treated groups did not differ significantly from those of the controls (Figure 7). Thus, the grazing treatment using BSA affected the oxidative-stress-related factors in sheepgrass only slightly.

Analysis of leaf-scar lengths

The lengths of the clipping scars on the leaves reflect changes in the sheepgrass transcriptome at the macroscopic level, especially changes in the cell oxidative status and apoptosis. For the first completely expanded leaf, the scars on the leaves daubed with water were approximately 1.3-fold longer than those on the leaves daubed with BSA (Figure 8), and this difference was statistically significant (p < 0.01). The phenotype of the OVA-treated leaves was similar to that of the water-treated leaves but markedly different from that of the BSA-treated leaves. However, the leaf-scar length on the last completely expanded leaf was similar in all treatments.

Figure 8
figure 8

Analysis of the leaf scars. (A) Images of the leaf scar measurements in the first and last complete leaves. A red line showing the position of leaf scars. (B) Statistical analyses. Left: data collected from the first complete leaves. Right: data collected from the last complete leaves. P** refers to a significant difference (p < 0.01) between the BSA and water treatments.


Grazing is an important and frequent stress for pasture and prairie plants. Plant scientists have long studied the effects of grazing on plants as a single process. However, grazing is a complex process that involves wounding effects caused by herbivore feeding, defoliation effects due to leaf-surface loss during grazing, and the deposition of herbivore saliva onto the surface of plants. Some studies have reported that wounding can stimulate plant growth but clearly differs from grazing [17, 4143]. Defoliation affects root development in grasses [44, 45] and alters the carbohydrate-metabolism pathway in rice (unpublished). On the other hand, scientists have examined how plants respond to BSA deposition for decades [11, 12]. Saliva has been found to stimulate plant growth, enhance tiller and increase biomass [14]. However little is known about the molecular mechanisms of grazing responses and the genetic and functional differences among three components of grazing. To investigate the genetic profile of the grazing response in sheepgrass and to elucidate the differences in mechanism between the saliva-deposition response and the responses to other grazing components, we analyzed the transcriptomes of control, defoliated, and grazed plants using RNA-seq.

Sequence quality and annotation

Illumina RNA-seq technology had been widely used in genome-wide analyses of cotton (Gossypium hirsutum), radish (Raphanus sativus), Brassica juncea, and Brassica pekinensis[4649]. Here, we obtained more than five million raw reads from most samples. Using the Trinity transcriptome-assembly software, a total of 120,426 assembled transcripts were obtained from seven sample libraries (untreated control, defoliation, and simulated grazing). Of these transcripts, 14,240 were annotated by BLASTX and functional-bioinformatics analyses, including the GO, KOG, and KEGG databases. No genome of sheepgrass or close relative species was available, so most of our transcripts cannot hit known proteins. In addition, a relative stringent blast parameter (E-value < 1e-5) might discard a part of known hits. We obtained 74,087 GO terms, 9,985 KOG terms and 7,240 KEGG pathway terms for all transcripts combined. The gene-transcription profiles of sheepgrass after grazing were stored in an annotated-gene catalog to provide a molecular understanding of grazing responses. The remaining un-annotated transcripts may represent a sheepgrass-specific gene pool. These results provide a solid foundation for further studies of the molecular mechanisms of grazing responses and for identifying grazing-related genes in this species.

Gene enrichment analysis and the effects of grazing

Grasslands and grass re-growth after grazing are very important for both the ecosystem and human dairy food supply on this planet. Grazing is a processing that have multiple components including at least wound, defoliation, and BSA deposition. Grazing often removes completely or partially the leaf part of plants. After grazing, plants have to transport the carbohydrate such as sucrose and other energy substance due to root sink demand [5052]. In our enrichment analysis (Table 2), “transporter activity” and “amino acid and derivative metabolism” in GO categories, “Amino acid transport and metabolism” and “Lipid transport and metabolism” in KOG categories, and KEGG amino acid metabolism pathways were significant enrichment in grazing treatment. The results indicated that amino acid metabolism involved in plant grazing response. The amino acid metabolism pathways may contribute to protein biosynthesis and plant recovery after grazing, and the related genes are worth further study.

Differential genetic profiles in response to BSA deposition

Sheepgrass is a dominant grassland species in northeastern China and Inner Mongolia and is known for its adaptability to grazing and excellent forage quality among perennial grasses [53]. In this study, 2,759 genes were expressed differently between the control and defoliation libraries, indicating that these genes responded to defoliation. Similarly, 3,156 genes were expressed differently between the control and grazing libraries, indicating that these genes responded to the combined effects of defoliation and BSA deposition. Furthermore, 2,002 genes were expressed differently between the defoliation and grazing libraries, indicating that these genes responded to BSA deposition. The only difference between these two treatments was the liquid deposited on leaves.

As shown in Table 3, most of the DEGs were down-regulated. In the KOG enrichment analysis (Table 5), the down-regulated DEGs were enriched in lipid transport and metabolism; energy production and conversion; and amino-acid transport and metabolism. Thus, several functionally linked metabolic pathways were down-regulated in response to grazing. This result is consistent with a proteomic analysis of rice after ovine BSA deposition [54], in which the authors found that most photosynthesis-related, energy-related, and carbohydrate-metabolism related proteins were down-regulated. BSA deposition on plants is accompanied by a multitude of stresses including oxidative stress, pathogenesis, and wounding [5557].

We examined transcript expression in sheepgrass at three time points following grazing. The gene-expression analysis helped to clarify how the expression of the DEGs adjusts to grazing stress. When G2 (2 h after treatment), G6 (6 h after treatment), and G24 (24 h after treatment) were compared to C (no treatment), 2,367 genes were differentially expressed after 2 h, 2,285 genes were differentially expressed after 6 h, and 1,692 genes were differentially expressed after 24 h. Among these, 1,074 genes were differentially expressed at all three time points, indicating that about half of the DEGs were stable during the first 24 h after grazing. The key pathways involved in grazing-response mechanisms may contain these genes.

Apoptosis-related DEGs

In plants, apoptosis is induced by multiple stresses, including salt, nitric oxide, oxidative stress, and wounding [5861]. This study also shows apoptosis-related DEGs in response to grazing. Based on functional-enrichment analysis, the apoptosis pathway is significantly involved in the saliva-deposition response. In the DNA-fragmentation experiment, we found fewer DNA fragments and delayed DNA fragmentation following BSA deposition compared to defoliation.

The apoptosis-promoting RNA-binding proteins TIA-1 and TIAR (RRM superfamily) was detected among the DEGs. These proteins promote DNA fragmentation in digitonin-permeabilized thymocytes and are pro-apoptotic factors that influence some aspect of RNA metabolism [62]. In our expression-mode analysis, TIA-1/TIAR was significantly down–regulated in the grazing treatment compared to the control and defoliation treatments. Correspondingly, less DNA ladder was seen after BSA deposition.

The expression of the serine/threonine kinase PAK4 increases the phosphorylation of the pro-apoptotic protein BAD and inhibits the activation of caspase, which protects cells against apoptosis [63]. Hsp70 and many other heat-shock proteins can overcome both caspase-dependent and caspase-independent apoptotic stimuli and confer immortality in various cell types [64]. Metacaspases are evolutionarily distant caspase homologs that are found outside the Metazoa and are known to play key roles in programmed cell death (PCD) [65]. However, whether metacaspases in plants function as caspases is controversial [66]. These apoptosis-inhibiting genes were all up-regulated only in the grazing treatment.

ATPases, including Na+, K+, and H+ ATPases, play critical roles in apoptosis [67, 68]. Apoptotic stimuli impair Na+- and K+-ATPase activity as a mechanism of neuronal death mediated by concurrent ATP deficiency and oxidative stress [69]. Several genes were annotated as “apoptotic ATPases” in the KOG functional analysis and their expression modes are shown in Table 7.

Additional apoptosis- or PCD-related genes were detected among the DEGs. Cullin controls non-lysosomal-mediated protein degradation and thus cell death [70]. Aldehyde dehydrogenase-2 (ALDH2) converts acetaldehyde into acetate, and over-expression of an ALDH2 transgene prevents acetaldehyde-induced cell injury and apoptosis [71]. The exosome-complex exonuclease rrp40 forms part of the exosome, which is important to the RNA-processing machinery of eukaryotes and functions in RNA degradation in both the nucleus and the cytoplasm [72]. A Ca2+-dependent cysteine protease (CDP) is associated with anoxia-induced root-tip death in maize [73].

Programmed cell death or apoptosis is an integral part of plant ontogenesis and plays a fundamental role in plant development. According to the above described genes, we suggest that there are some important pathways of apoptosis, in response to BSA deposition during grazing in plants.

Oxidative-stress-related DEGs

Cellular oxidative stress is a common challenge for plants that usually accompanies wounding [74] or senescence [75] and is closely associated with apoptosis [76]. Many apoptosis-inducing agents are either oxidants or stimulators of cellular oxidative metabolism. In our measurements of cellular oxidative status, H2O2, and MDA but not SOD increased substantially after grazing. However, cells from sheepgrass leaves subjected to BSA deposition showed significantly low H2O2 and MDA levels. In addition, we detected cellular oxidative-control genes among the DEGs.

Major ROS-scavenging enzymes in plants include superoxide dismutase (SOD), ascorbate peroxidase (APX), catalase (CAT), glutathione peroxidase (GPX), and peroxiredoxin (PrxR) [77]. We found no SOD genes among the DEGs, indicating that SOD expression was stable following the grazing process. This finding is consistent with the cellular oxidative-status experiment. Furthermore, we found no PrxR DEGs. However, the DEGs included APX, CAT, and GPX genes (Table 7). These three ROS-scavenging enzymes were up-regulated in the grazing treatment compared to the control. Their expression was also relatively higher in the grazing treatment than in the defoliation treatment, possibly explaining the low H2O2 and MDA levels after BSA deposition.

The DEGs included other cellular oxidative-stress-related genes. Ferritin prevents the formation of the highly toxic HO·radical via the metal-dependent Haber–Weiss reaction or the Fenton reaction [78]. Several studies have shown that biotic and abiotic stresses are accompanied by an oxidative burst mediated by NADPH oxidases [79, 80]. The glutaredoxin (Grx) and thioredoxin (Trx) pathways use NADPH to reduce the disulfide bonds that form in some cytoplasmic enzymes during catalysis. The thioredoxin system consists of thioredoxin reductase and thioredoxin, and the glutaredoxin system is composed of glutathione reductase, glutathione, and three glutaredoxins [81]. The differential expression of these genes suggests that grazing stress or BSA deposition was closely related to cellular oxidative changes.The present transcriptome-sequencing results and related biochemical experiments help to elucidate the response of sheepgrass to BSA deposition. After grazing, the plants receive the signal of BSA deposited on the leaves and elevate the expression of ROS-scavenging enzymes and antioxidant pathways to respond to the subsequent oxidative burst in the cells. The decreased cellular oxidative levels result in fewer apoptotic cells in the grazing wound, accelerating the recovery from the grazing stress in plants. A model of sheepgrass responded to BSA deposition was indicated in Figure 9. The co-evolution of plants and herbivores in grazing systems may have led to this antioxidant mechanism in response to BSA deposition.

Figure 9
figure 9

Molecular mechanism of sheepgrass responds to BSA deposition. The DEGs in red oval frames are up-regulated and those in blue oval frames are down-regulated in expression levels in grazing libraries compared with defoliation libraries. The arrow lines stand for the effect of activation. The blunt lines stand for the effects of inhibition. The dotted lines stand for the unknown effects.


To investigate the molecular mechanism of grazing responses, we performed transcriptome sequencing and analysis to identify DEGs in sheepgrass subjected to simulated BSA deposition. Our results show that BSA deposition triggers differential gene expression compared to defoliation and other grazing components. Based on a functional analysis of the saliva-deposition DEGs, the cellular-antioxidant and apoptotic pathways apparently respond to grazing stress. Macroscopic changes confirm the effects of these two pathways in sheepgrass. Although the connection between the two pathways requires further evidence, we believe that the saliva-deposition-induced pathways work together to contribute to plant recovery after grazing.


Plant materials, growth conditions, and treatments

All sheepgrass plants (Zhongke No. 3) were obtained from the field 4 weeks before the initiation of the experiment. Seedlings in good condition were collected and transplanted into trays in the greenhouse. The trays were filled with a mixture of vermiculite and commercial potting soil at a ratio of 1:2. The plates were placed in a greenhouse at 25°C and 70% relative humidity. The aboveground portion of each plant was cut off, and plants were allowed to re-grow to the 5- or 6-leaf stage. To initiate the experimental treatments, two-thirds of the aboveground portion of each plant was cut off. For the defoliation treatment, water was then daubed on the cut ends of the leaves. For the grazing treatment, a BSA solution (1 mM) was daubed on the cut ends of the leaves. The remaining aboveground portions of the plants were collected 2, 6, and 24 h after the clipping and daubing treatments. The corresponding parts of the control (unclipped) seedlings were collected at the same time. All harvested seedlings of each treatment were immediately frozen in liquid nitrogen and stored at −80°C. The clipping and daubing treatments were conducted at 10:00 AM and all the materials collection was conducted in daylight hours to reduce the effort of circadian rhythmicity. In total, seven samples were obtained: C (control); D2, D6, and D24 (2, 6, and 24 h after defoliation, respectively); and G2, G6, and G24 (2, 6, and 24 h after grazing, respectively).

RNA-seq library preparation and Illumina sequencing

Total RNA was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and NucleoSpin® RNA Clean-up kit (CapitalBio Company, China) according to the manufacturer’s instructions. The RNA quality was assessed by agarose-gel electrophoresis, and the RNA absorbance was measured by spectrophotometry. The ratio of the absorbance at 260/280 nm was then used to determine the RNA quality. 6 μg RNA of each sample was used for transcriptome sequencing. The RNA was processed for use on an RNA-seq platform (Illumina, Inc, San Diego, CA, USA) by the Chinese National Human Genome Center at Shanghai (CHGC).

mRNA [poly(A) RNA] was then purified from total RNA using Micropoly(A)Purist™ mRNA purification kit (Ambion, Cat.No.1919, Foster, CA, USA). The mRNA was fragmented and converted into a RNA-seq library using the mRNAseq library construction kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s instructions. 2×100 bp paired-end sequencing was performed using the Illumina Genome Analyzer II x (Illumina GAIIx, San Diego, CA, USA).

Sequence filtering and assembly

Sequence reads from all samples were cleaned using the FASTX toolkit ( First all the reads containing ‘N’ were discarded using a perl script, then adapter sequences were removed using the fastx_clipper program, followed by removal of quality < 5 bases from the 3′ end with fastq_quality_trimmer, requiring a minimum sequence length of 50 bp. Finally the reads with at least 90% bases > quality 20 were chosen using fastq quality filter for further assembly. De novo transcriptome assembly was performed using Trinity RNA-seq assembly v2013-02-25 with default parameters [PMID: 21572440].

Meanwhile, all sequence reads generated by Illumina sequencing were aligned to the reference transcriptome dataset using SOAP2 software [82].

DEG identification

Based on the Trinity assembly results, the number of reads for each contig from each sample (control, defoliation, and grazing) was converted to reads per kilobase per million (RPKM) [83]. The MARS (MA-plot-based method with random-sampling model) module in the DEGseq package was used to calculate the differential expression of each contig between the analyzed samples [40]. The package DEGseq is a free R package for identifying DEGs from RNA-seq data. We used an FDR (false discovery rate) to determine the threshold p-value. An FDR < 0.001 was considered to indicate a significant difference in expression between the control and treated samples.

Functional annotation, classification, and pathway analysis of DEGs

The sheepgrass transcriptome sequencing data had previously undergone Gene ontology (GO) annotation using a BLASTP search against the Swiss-Prot and TrEMBL databases with an E-value ≤ 1e-5 [84]. A GO functional classification was performed using WEGO software ( to understand the distribution of gene functions in grazed sheepgrass [85].

In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation was performed using the KEGG Automatic Annotation Server (KAAS) with the bi-directional best-hit method. For KEGG, KAAS annotates every submitted sequence with a KEGG orthology (KO) identifier representing an orthologous group of genes directly linked to an object in the KEGG pathways and BRITE functional hierarchy [86, 87]. KEGG pathway enrichment analysis was conducted using KOBAS 2.0 (, [88]).

Finally, EuKaryotic Orthologous Groups (KOG) annotation was carried out against the NCBI KOG database with a typical cut-off E-value ≤ 1e-5. The KOG annotations of the DEGs were classified into 25 protein functions and compared between the defoliation and grazing treatments. KOG enrichment analysis was conducted through hypergeometric distribution testing using the Phyper function in the R software package ( The Bonferroni correction was used to adjust the p-values. The significantly enriched functional clusters were selected based on a corrected q-value (< 0.05).

DNA extraction and DNA-fragmentation assay

Using the same plant material, the cut ends of the leaves were collected every day for 9 days after clipping and daubing with water or BSA. Unclipped leaves were used as controls. Each sample for DNA extraction consisted of ten 1-cm-long leaf pieces. The leaves were flash-frozen in liquid nitrogen and stored at −80°C. The leaf tissues were ground to a fine powder in liquid nitrogen using a mortar and pestle, and the DNA was extracted using a plant genomic-DNA extraction kit (TIANGEN, Beijing, China) according to the manufacturer’s protocol. The total DNA was treated with RNase A (TaKaRa Bio Inc., Dalian, China) to remove any contaminating RNA. The isolated DNA, which was mainly derived from apoptotic cell bodies, was electrophoresed on a 2.0% agarose gel at 50 V for 3 h. The DNA fragments, which consisted of 160–200 bp multimers, were visualized under ultraviolet light after staining with ethidium bromide.

Measurement of cell oxidative factors concentrations in the leaves

Using the same plant material, the leaves were clipped and daubed with distilled water, BSA solution or OVA solution. Three independent replicates were performed. Three days later, 1 cm of the cut end of each completely expanded leaf was collected, and approximately 30 cut ends were pooled. We chose the first complete leaf ends for further experiments. The corresponding parts of unclipped seedlings were collected as controls. The samples were frozen in liquid nitrogen and ground to a fine powder using a mortar and pestle. Approximately 0.2 g of each crude extract was added to 5 ml of pre-chilled PBS (50 mM at pH 7.8), thoroughly mixed, and centrifuged for 20 min at 4,000 g at 4°C. The supernatants were simultaneously assayed using the TBA (to measure MDA), KI (to measure H2O2), and NBT (to measure SOD) methods [8991]. The absorbance values were measured using a 2600 UV spectrophotometer (UNICO, Shanghai, China).

Leaf-scar measurements and phenotype comparisons

Using the same plant material, two-thirds of the first and last completely expanded leaves on each seedling were cut off at 10:00 AM. After clipping, the cut surfaces of the leaves were immediately daubed with water, 1 mM BSA or 1 mM OVA. Approximately 120 completely expanded leaves were removed per treatment. Three days later, the leaf scars were measured, and statistical analyses were performed to evaluate the effects of BSA on the first and last completely expanded leaves. Three independent replicates were performed. The statistical analyses were conducted using SAS 9.0 (SAS Institute, Cary, NC, USA) to compare the differences among the three treatments.


  1. Chen S, Li X, Zhao A, Wang L, Li X, Shi Q, Chen M, Guo J, Zhang J, Qi D, Liu G: Genes and pathways induced in early response to defoliation in rice seedlings. Curr Issues Mol Biol. 2009, 11: 81-100.

    CAS  PubMed  Google Scholar 

  2. Gold WG, Caldwell MM: The effects of the spatial pattern of defoliation on regrowth of a tussock grass. Oecologia. 1990, 82: 12-17. 10.1007/BF00318527.

    Article  Google Scholar 

  3. Orodho AB, Trlica MJ: Clipping and long-term grazing effects on biomass and carbohydrate reserves of Indian ricegrass. J Range Manage. 1990, 43: 52-57. 10.2307/3899121.

    Article  Google Scholar 

  4. Visser RD, Vianden H, Schnyder H: Kinetics and relative significance of remobilized and current C and N incorporation in leaf and root growth zones of Lolium perenne after defoliation: assessment by 13C and 15N steady-state labelling. Plant Cell Environ. 1997, 20: 37-46. 10.1046/j.1365-3040.1997.d01-9.x.

    Article  Google Scholar 

  5. Macduff JH, Humphreys MO, Thomas H: Effects of a stay-green mutation on plant nitrogen relations in Lolium perenne during N starvation and after defoliation. Ann Bot. 2002, 89: 11-21. 10.1093/aob/mcf001.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  6. Christopher ME, Miranda M, Major IT, Constabel CP: Gene expression profiling of systemically wound-induced defenses in hybrid poplar. Planta. 2004, 219: 936-947. 10.1007/s00425-004-1297-3.

    Article  CAS  PubMed  Google Scholar 

  7. Parsons AJ, Leafe EL, Collett B, Stiles W: The physiology of grass production under grazing. I. Characteristics of leaf and canopy photosynthesis of continuously- grazed swards (perennial ryegrass Lolium perenne). Ying Yong Sheng Tai Xue Bao. 1983, 20: 117-126.

    Google Scholar 

  8. Jarvis SC, Macduff JH: Nitrate nutrition of grasses from steady-state supplies in flowing solution culture following nitrate deprivation and/or defoliation. I. Recovery of uptake and growth and their interactions. J Exp Bot. 1989, 40: 965-975. 10.1093/jxb/40.9.965.

    Article  Google Scholar 

  9. Clement CR, Hopper MJ, Jones LHP, Leafe EL: The uptake of nitrate by Lolium perenne from flowing nutrient solution: II. Effect of light, defoliation, and relationship to co2 flux. J Exp Bot. 1978, 29: 1173-1183. 10.1093/jxb/29.5.1173.

    Article  CAS  Google Scholar 

  10. Mousavi SA, Chauvin A, Pascaud F, Kellenberger S, Farmer EE: GLUTAMATE RECEPTOR-LIKE genes mediate leaf-to-leaf wound signalling. Nature. 2013, 500: 422-426. 10.1038/nature12478.

    Article  CAS  PubMed  Google Scholar 

  11. Vittoria A, Rendina N: Fattori condizionanti la funzionalita tiaminica in piante superiori e cenni sugli effetti dell bocca dei runinanti sull erbe pascolative. Acta Medica Vet (Naples). 1960, 6: 379-405.

    Google Scholar 

  12. Johnston A, Bailey CD: Influence of bovine saliva on grass regrowth in the greenhouse. Can J Anim Sci. 1972, 52: 573-574. 10.4141/cjas72-070.

    Article  Google Scholar 

  13. Reardon PO, Leinweber CL, Merrill LB: The effect of bovine saliva on grasses. J Anim Sci. 1972, 34: 897-898.

    Google Scholar 

  14. Reardon PO, Leinweber CL, Merrill LB: Response of sideoats grama to animal saliva and thiamine. J Range Manage. 1974, 27: 400-401. 10.2307/3896502.

    Article  Google Scholar 

  15. Liu J, Wang L, Wang D, Bonser SP, Sun F, Zhou Y, Gao Y, Teng X: Plants can benefit from herbivory: stimulatory effects of sheep saliva on growth of Leymus chinensis. PLoS ONE. 2012, 7: e29259-10.1371/journal.pone.0029259.

    Article  PubMed Central  PubMed  Google Scholar 

  16. Musser RO, Farmer E, Peiffer M, Williams SA, Felton GW: Ablation of caterpillar labial salivary glands: technique for determining the role of saliva in insect-plant interactions. J Chem Ecol. 2006, 32: 981-992. 10.1007/s10886-006-9049-4.

    Article  CAS  PubMed  Google Scholar 

  17. Mattiacci L, Dicke M, Posthumus MA: beta-Glucosidase: an elicitor of herbivore-induced plant odor that attracts host-searching parasitic wasps. Proc Natl Acad Sci U S A. 1995, 92: 2036-2040. 10.1073/pnas.92.6.2036.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  18. Cohen S: Isolation of a mouse submaxillary gland protein accelerating incisor eruption and eyelid opening in the new-born animal. J Biol Chem. 1962, 237: 1555-1562.

    CAS  PubMed  Google Scholar 

  19. McNaughton SJ: Interactive regulation of grass yield and chemical properties by defoliation, a salivary chemical, and inorganic nutrition. Oecologia. 1985, 65: 478-486. 10.1007/BF00379660.

    Article  Google Scholar 

  20. Kato R, Nagayama E, Suzuki T, Uchida K, Shimomura TH, Harada Y: Promotion of plant cell division by an epidermal growth factor. Plant Cell Physiol. 1993, 34: 789-793.

    CAS  Google Scholar 

  21. Cheung F, Haas BJ, Goldberg SMD, May GD, Xiao Y, Town CD: Sequencing Medicago truncatula expressed sequenced tags using 454 Life Sciences technology. BMC Genomics. 2006, 7: 272-10.1186/1471-2164-7-272.

    Article  PubMed Central  PubMed  Google Scholar 

  22. Weber APM, Weber KL, Carr K, Wilkerson C, Ohlrogge JB: Sampling the Arabidopsis transcriptome with massively parallel pyrosequencing. Plant Physiol. 2007, 144: 32-42. 10.1104/pp.107.096677.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  23. Emrich SJ, Barbazuk WB, Li L, Schnable PS: Gene discovery and annotation using LCM-454 transcriptome sequencing. Genome Res. 2007, 17: 69-73.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  24. Steuernagel B, Taudien S, Gundlach H, Seidel M, Ariyadasa R, Schulte D, Petzold A, Felder M, Graner A, Scholz U, Mayer KF, Platzer M, Stein N: De novo 454 sequencing of barcoded BAC pools for comprehensive gene survey and genome analysis in the complex genome of barley. BMC Genomics. 2009, 10: 547-10.1186/1471-2164-10-547.

    Article  PubMed Central  PubMed  Google Scholar 

  25. Deschamps S, Rota M, Ratashak JP, Biddle P, Thureen D, Famer A, Luck S, Beatty M, Nagasawa N, Michael L, Liaca V, Sakai H, May G, Lightner J, Campbell MA: Rapid genome-wide single nucleotide polymorphism discovery in soybean and rice via deep resequencing of reduced representation libraries with the Illumina Genome Analyzer. Plant Genome. 2010, 3: 53-68. 10.3835/plantgenome2009.09.0026.

    Article  CAS  Google Scholar 

  26. Garg R, Patel RK, Jhanwar S, Priya P, Bhattacharjee A, Yadav G, Bhatia S, Chattopadhyay D, Tyaqi AK, Jain M: Gene discovery and tissue-specific transcriptome analysis in chickpea with massively parallel pyrosequencing and web resource development. Plant Physiol. 2011, 156: 1661-1678. 10.1104/pp.111.178616.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  27. Hiremath PJ, Farmer A, Cannon SB, Woodward J, Kudapa H, Tuteja R, Knmar A, Bhanuprakash A, Mulaosmanovic B, Gujaria N, Krishnamurthy L, Gaur PM, Kavikishor PB, Shah T, Srinivasan R, Lohse M, Xiao Y, Town CD, Cook DR, May GD, Varshney RK: Large-scale transcriptome analysis in chickpea (Cicer arietinum L.), an orphan legume crop of the semi-arid tropics of Asia and Africa. Plant Biotechnol J. 2011, 9: 922-931. 10.1111/j.1467-7652.2011.00625.x.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  28. Troncoso-Ponce MA, Kilaru A, Cao X, Durrett TP, Fan J, Jensen JK, Thrower NA, Pauly M, Wilkerson C, Ohlrogge JB: Comparative deep transcriptional profiling of four developing oilseeds. Plant J. 2011, 68: 1014-1027. 10.1111/j.1365-313X.2011.04751.x.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  29. Sun Y, Wang F, Wang N, Dong Y, Liu Q, Zhao L, Chen H, Liu W, Yin H, Zhang X, Yuan Y, Li H: Transcriptome exploration in Leymus chinensis under saline-alkaline treatment using 454 pyrosequencing. PLoS ONE. 2013, 8: e53632-10.1371/journal.pone.0053632.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  30. Chen S, Huang X, Yan X, Liang Y, Wang Y, Li X, Peng X, Ma X, Zhang L, Cai Y, Ma T, Cheng L, Qi D, Zheng H, Yang X, Li X, Liu G: Transcriptome analysis in sheepgrass (Leymus chinensis): a dominant perennial grass of the eurasian steppe. PLoS ONE. 2013, 8: e67974-10.1371/journal.pone.0067974.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  31. Li X, Gao Q, Liang Y, Ma T, Cheng L, Qi D, Liu H, Xu X, Chen S, Liu G: A novel salt-induced gene from sheepgrass, LcSAIN2, enhances salt tolerance in transgenic Arabidopsis. Plant Physiol Biochem. 2013, 64: 52-59.

    Article  PubMed  Google Scholar 

  32. Li X, Hou S, Gao Q, Zhao P, Chen S, Qi D, Lee B, Cheng L, Liu G: LcSAIN1, a novel salt-induced gene from sheepgrass, confers salt stress tolerance in transgenic Arabidopsis and rice. Plant Cell Physiol. 2013, 54: 1172-1185. 10.1093/pcp/pct069.

    Article  CAS  PubMed  Google Scholar 

  33. Cheng L, Li X, Huang X, Ma T, Liang Y, Ma X, Peng X, Jia J, Chen S, Chen Y, Deng B, Gongshe Liu G: Overexpression of Sheepgrass R1-MYB transcription factor LcMYB1 confers salt tolerance in transgenic arabidopsis. Plant Physiol Biochem. 2013, 70: 252-260.

    Article  CAS  PubMed  Google Scholar 

  34. Peng X, Zhang L, Zhang L, Liu Z, Cheng L, Yang Y, Shen S, Chen S, Liu G: The transcriptional factor LcDREB2 cooperates with LcSAMDC2 to contribute to salt tolerance in Leymus chinensis. Plant Cell Tissue Organ Cult. 2013, 113: 245-256. 10.1007/s11240-012-0264-0.

    Article  CAS  Google Scholar 

  35. Liu Z, Chen Z, Pan J, Li X, Su M, Wang L, Li H, Liu G: Mol Phylogenet Evol. 2008, 46: 278-289. 10.1016/j.ympev.2007.10.009.

    Article  CAS  PubMed  Google Scholar 

  36. Chen S, Cai Y, Zhang L, Yan X, Cheng L, Qi D, Zhou Q, Li X, Liu G: Transcriptome analysis reveals common and distinct mechanisms for sheepgrass (Leymus chinensis) responses to defoliation compared to mechanical wounding. PLoS ONE. 2014, 9: e89495-10.1371/journal.pone.0089495.

    Article  PubMed Central  PubMed  Google Scholar 

  37. Su M, Li X, Li X, Cheng L, Qi M, Chen S, Liu G: Molecular characterization and expression analysis of a Sheepgrass sucrose transporter LcSUT1 after defoliation. Plant Mol Biol Report. 2013, 31: 1184-1191. 10.1007/s11105-013-0582-3.

    Article  CAS  Google Scholar 

  38. Peng X, Ma X, Fan W, Su M, Cheng L, Alam I, Lee B, Qi D, Shen S, Liu G: Improved drought and salt tolerance of Arabidopsis thaliana by transgenic expression of a novel DREB gene from Leymus chinensis. Plant Cell Rep. 2011, 30: 1493-1502. 10.1007/s00299-011-1058-2.

    Article  CAS  Google Scholar 

  39. Lamy E, Costa G, Santos R, Silva FC, Potes J, Pereira A, Coelho AV, Baptista ES: Effect of condensed tannin ingestion in sheep and goat parotid saliva proteome. J Anim Physiol Anim Nutr (Berl). 2010, 95: 304-312.

    Article  Google Scholar 

  40. Wang L, Feng Z, Wang X, Wang X, Zhang X: DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010, 26: 136-138. 10.1093/bioinformatics/btp612.

    Article  PubMed  Google Scholar 

  41. White LM: Carbohydrate reserves of grasses: a review. J Range Manage. 1973, 26: 13-18. 10.2307/3896873.

    Article  CAS  Google Scholar 

  42. Baldwin IT: Herbivory simulations in ecological research. Trends Ecol Evol. 1990, 5: 91-93. 10.1016/0169-5347(90)90237-8.

    Article  CAS  PubMed  Google Scholar 

  43. Capinera JL, Roltsch WJ: Response of wheat seedlings to actual and simulated migratory grasshopper defoliation. J Econ Entomol. 1980, 73: 258-261.

    Article  Google Scholar 

  44. Thornton B, Millard P: Effects of severity of defoliation on root functioning in grasses. J Range Manage. 1996, 49: 443-447. 10.2307/4002927.

    Article  Google Scholar 

  45. Mackie-Dawson LA: Nitrogen uptake and root morphological responses of defoliated Loliumperenne (L.) to a heterogeneous nitrogen supply. Plant Soil. 1999, 209: 111-118. 10.1023/A:1004534609280.

    Article  CAS  Google Scholar 

  46. Wang F, Li H, Zhang Y, Li J, Li L, Liu L, Wang L, Wang C, Gao J: MicroRNA expression analysis of rosette and folding leaves in Chinese cabbage using high-throughput Solexa sequencing. Gene. 2013, 532: 222-229. 10.1016/j.gene.2013.09.039.

    Article  CAS  PubMed  Google Scholar 

  47. Liu X, Lu Y, Yuan Y, Liu S, Guan C, Chen S, Liu Z: De novo transcriptome of Brassica juncea seed coat and identification of genes for the biosynthesis of flavonoids. PLoS ONE. 2013, 8: e71110-10.1371/journal.pone.0071110.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  48. Kong X, Luo Z, Dong H, Eneji AE, Li W, Lu H: Gene expression profiles deciphering leaf senescence variation between early- and late-senescence cotton lines. PLoS ONE. 2013, 8: e69847-10.1371/journal.pone.0069847.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  49. Wang Y, Xu L, Chen Y, Shen H, Gong Y, Limera C, Liu L: Transcriptome profiling of radish (Raphanus sativus L.) root and identification of genes involved in response to Lead (Pb) stress with next generation sequencing. PLoS ONE. 2013, 8: e66539-10.1371/journal.pone.0066539.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  50. Amiard V, Morvan-Bertrand A, Billard JP, Huault C, Prud’homme MP: Fate of fructose supplied to leaf sheaths after defoliation of Lolium perenne L.: assessment by 13C-fructose labelling. J Exp Bot. 2003, 54: 1231-1243. 10.1093/jxb/erg125.

    Article  CAS  PubMed  Google Scholar 

  51. Morvan-Bertrand A, Boucaud J, Prud’homme MP: Influence of initial levels of carbohydrates, fructans, nitrogen, and soluble proteins on regrowth of Lolium perenne L. cv. Bravo following defoliation. J Exp Bot. 1999, 50: 1817-1826. 10.1093/jxb/50.341.1817.

    Article  CAS  Google Scholar 

  52. Morvan-Bertrand A, Pavis N, Boucaud J, Prud’Homme MP: Partitioning of reserve and newly assimilated carbon in roots and leaf tissues of Lolium perenne during regrowth after defoliation: assessment by 13C steady-state labelling and carbohydrate analysis. Plant Cell Environ. 1999, 22: 1097-1108. 10.1046/j.1365-3040.1999.00485.x.

    Article  CAS  Google Scholar 

  53. Bai YF, Han XG, Wu JG, Chen ZZ, Li LH: Ecosystem stability and compensatory effects in the Inner Mongolia grassland. Nature. 2004, 431: 181-184. 10.1038/nature02850.

    Article  CAS  PubMed  Google Scholar 

  54. Fan WH, Cui W, Li X, Chen S, Liu S, Shen S: Proteomics analysis of rice seedlingresponses to ovine saliva. J Plant Physiol. 2010, 168: 500-509.

    Article  PubMed  Google Scholar 

  55. Agrawal GK, Rakwal R, Yonekura M, Kubo A, Saji H: Proteome analysis of differentially displayed proteins as a tool for investigating ozone stress in rice (Oryza sativa L.) seedlings. Proteomic. 2002, 2: 947-959. 10.1002/1615-9861(200208)2:8<947::AID-PROT947>3.0.CO;2-J.

    Article  CAS  Google Scholar 

  56. Mahmood T, Jan A, Kakishima M, Komatsu S: Proteomic analysis of baeterial-blight defense-responsive proteins in rice leaf blades. Proteomics. 2006, 6: 6053-6065. 10.1002/pmic.200600470.

    Article  CAS  PubMed  Google Scholar 

  57. Wei Z, Hu W, Lin Q, Cheng X, Tong M, Zhu L, Chen R, He G: Understanding rice plant resistance to the Brown Plant hopper (Nilaparvata lugens): a proteomic approach. Proteomics. 2009, 9: 2798-2808. 10.1002/pmic.200800840.

    Article  CAS  PubMed  Google Scholar 

  58. Maki K: Apoptosis-like cell death in barley roots under salt stress. Plant Cell Physiol. 1997, 38: 1091-1093. 10.1093/oxfordjournals.pcp.a029277.

    Article  Google Scholar 

  59. Pedroso MC, Magalhaes JR, Durzan D: A nitric oxide burst precedes apoptosis in angiosperm and gymnosperm callus cells and foliar tissues. J Exp Bot. 2000, 51: 1027-1036. 10.1093/jexbot/51.347.1027.

    Article  CAS  PubMed  Google Scholar 

  60. Averyanov A: Oxidative burst and plant disease resistance. Front Biosci (Elite Ed). 2009, 1: 142-152. 10.2741/s14.

    Article  Google Scholar 

  61. Sanchez P, de Torres Zabala M, Grant M: AtBI-1, a plant homologue of Baxinhibitor-1, suppresses Bax-induced cell death in yeast and is rapidly upregulated during wounding and pathogen challenge. Plant J. 2000, 21: 393-399. 10.1046/j.1365-313x.2000.00690.x.

    Article  CAS  PubMed  Google Scholar 

  62. Förch P, Valcárcel J: Molecular mechanisms of gene expression regulation by the apoptosis-promoting protein TIA-1. Apoptosis. 2001, 6: 463-468. 10.1023/A:1012441824719.

    Article  PubMed  Google Scholar 

  63. Nerina G, Qu J, Audrey M: The Serine/Threonine kinase PAK4 prevents caspase activation and protects cells from apoptosis. J Biol Chem. 2001, 276: 14414-14419.

    Google Scholar 

  64. Sõti C, Sreedhar AS, Csermely P: Apoptosis, necrosis and cellular senescence: chaperone occupancy as a potential switch. Aging Cell. 2003, 2: 39-45. 10.1046/j.1474-9728.2003.00031.x.

    Article  PubMed  Google Scholar 

  65. Choi CJ, Berges JA: New types of metacaspases in phytoplankton reveal diverse origins of cell death proteases. Cell Death Dis. 2013, 4: e490-10.1038/cddis.2013.21.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  66. Tsiatsiani L, Van Breusegem F, Gallois P, Zavialov A, Lam E, Bozhkov PV: Metacaspases. Cell Death Differ. 2011, 18: 1279-1288. 10.1038/cdd.2011.66.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  67. Okahashi N, Nakamura I, Jimi E, Koide M, Suda T, Nishihara T: Specific inhibitors of vacuolar H(+)-ATPase trigger apoptotic cell death of osteoclasts. J Bone Miner Res. 1997, 12: 1116-1123. 10.1359/jbmr.1997.12.7.1116.

    Article  CAS  PubMed  Google Scholar 

  68. Aravind L, Dixit VM, Koonin EV: The domains of death: evolution of the apoptosis machinery. Trends Biochem Sci. 1999, 24: 47-53. 10.1016/S0968-0004(98)01341-3.

    Article  CAS  PubMed  Google Scholar 

  69. Wang XQ, Xiao AY, Sheline C, Hyrc K, Yang A, Goldberg MP, Choi DW, Yu SP: Apoptotic insults impair Na+, K+−ATPase activity as a mechanism of neuronal death mediated by concurrent ATP deficiency and oxidant stress. J Cell Sci. 2003, 116: 2099-2110. 10.1242/jcs.00420.

    Article  CAS  PubMed  Google Scholar 

  70. Rai V, Dey N: Identification of programmed cell death related genes in bamboo. Gene. 2012, 497: 243-248. 10.1016/j.gene.2012.01.018.

    Article  CAS  PubMed  Google Scholar 

  71. Li SY, Gomelsky M, Duan J, Zhang Z, Gomelsky L, Zhang X, Epstein PN, Ren J: Overexpression of aldehyde dehydrogenase-2 (ALDH2) transgene prevents acetaldehyde-induced cell injury in human umbilical vein endothelial cells: role of ERK and p38 mitogen-activated protein kinase. J Biol Chem. 2004, 279: 11244-11252. 10.1074/jbc.M308011200.

    Article  CAS  PubMed  Google Scholar 

  72. Houseley J, LaCava J, Tollervey D: RNA-quality control by the exosome. Nat Rev Mol Cell Biol. 2006, 7: 529-539. 10.1038/nrm1964.

    Article  CAS  PubMed  Google Scholar 

  73. Subbaiah CC, Kollipara KP, Sachs MM: A Ca(2+)-dependent cysteine protease is associated with anoxia-induced root tip death in maize. J Exp Bot. 2000, 51: 721-730. 10.1093/jexbot/51.345.721.

    Article  CAS  PubMed  Google Scholar 

  74. Bradley DJ, Kjellbom P, Lamb CJ: Elicitor- and wound-induced oxidative cross-linking of a proline-rich plant cell wall protein: a novel, rapid defense response. Cell. 1992, 70: 21-30. 10.1016/0092-8674(92)90530-P.

    Article  CAS  PubMed  Google Scholar 

  75. Prigione A, Fauler B, Lurz R, Lehrach H, Adjaye J: The senescence-related mitochondrial/oxidative stress pathway is repressed in human induced pluripotent stem cells. Stem Cells. 2010, 28: 721-733. 10.1002/stem.404.

    Article  CAS  PubMed  Google Scholar 

  76. Buttke TM, Sandstrom PA: Oxidative stress as a mediator of apoptosis. Immunol Today. 1994, 15: 7-10. 10.1016/0167-5699(94)90018-3.

    Article  CAS  PubMed  Google Scholar 

  77. Mittler R, Vanderauwera S, Gollery M, Van Breusegem F: Reactive oxygen gene network of plants. Trends Plant Sci. 2004, 9: 490-498. 10.1016/j.tplants.2004.08.009.

    Article  CAS  PubMed  Google Scholar 

  78. Bowler C, Slooten L, Vandenbranden S, De Rycke R, Botterman J, Sybesma C, Van Montagu M, Inzé DC: Manganese superoxide dismutase can reduce cellular damage mediated by oxygen radicals in transgenic plants. EMBO J. 1991, 10: 1723-1732.

    CAS  PubMed Central  PubMed  Google Scholar 

  79. Park HJ, Miura Y, Kawakita K, Yoshioka H, Doke N: Physiological mechanisms of a sub-systemic oxidative burst triggered by elicitor-induced local oxidative burst in potato tuber slices. Plant Cell Physiol. 1998, 39: 1218-1225. 10.1093/oxfordjournals.pcp.a029323.

    Article  CAS  Google Scholar 

  80. Cazale’ AC, Rouet-Mayer MA, Barboer-Brygoo H, Mathieu Y, Laurière C: Oxidative burst and hypoosmotic stress in tobacco cell suspensions. Plant Physiol. 1998, 116: 659-669. 10.1104/pp.116.2.659.

    Article  Google Scholar 

  81. Prinz WA, Aslund F, Holmgren A, Beckwith J: The role of the thioredoxin and glutaredoxin pathways in reducing protein disulfide bonds in the Escherichia coli cytoplasm. J Biol Chem. 1997, 272: 15661-15667. 10.1074/jbc.272.25.15661.

    Article  CAS  PubMed  Google Scholar 

  82. Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, Wang J: SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009, 25: 1966-1967. 10.1093/bioinformatics/btp336.

    Article  CAS  PubMed  Google Scholar 

  83. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5: 621-628. 10.1038/nmeth.1226.

    Article  CAS  PubMed  Google Scholar 

  84. Chen ZZ, Xue CH, Zhu S, Zhou FF, Xenfeng BL, Liu GP, Chen LB: GoPipe: streamlined gene ontology annotation for batch anonymous sequences with statistics. Prog Biochem Biophys. 2005, 32: 187-190.

    CAS  Google Scholar 

  85. Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, Wang J, Li S, Li R, Bolund L, Wang J: WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006, 34: W293-W297. 10.1093/nar/gkl031.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  86. Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M: KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007, 35: W182-W185. 10.1093/nar/gkm321.

    Article  PubMed Central  PubMed  Google Scholar 

  87. Mao X, Cai T, Olyarchuk JG, Wei L: Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005, 21: 3787-3793. 10.1093/bioinformatics/bti430.

    Article  CAS  PubMed  Google Scholar 

  88. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li CY, Wei L: KOBAS 2.0: a web server for annotation and identification of enriched pathways and disease. Nucleic Acids Res. 2011, 39: W316-W322. 10.1093/nar/gkr483.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  89. Liu CZ, Lei B: Effect of acupuncture on serum malonaldehyde content, superoxide dismutase and glutathione peroxidase activity in chronic fatigue syndrome rats. Zhen Ci Yan Jiu. 2012, 37: 38-40. 58

    PubMed  Google Scholar 

  90. Zhang Y, Lin L, Su B: Determination of trace H2O2 by KI iodine blue spectrophotometry. Chinese J Anal Lab. 2001, 20: 41-42.

    Google Scholar 

  91. Siddiqi ZA, Khalid M, Kumar S, Shahid M, Noor S: Antimicrobial and SOD activities of novel transition metal complexes of pyridine-2,6-dicarboxylic acid containing 4-picoline as auxiliary ligand. Eur J Med Chem. 2010, 45: 264-269. 10.1016/j.ejmech.2009.10.005.

    Article  CAS  PubMed  Google Scholar 

Download references


We are grateful to Huajun Zheng for his assistance in RNA sequencing and data support. This work was supported by the National Basic Research Program of China (“973”, 2014CB138704), the National High Technology Research and Development Program of China (“863”, 2011AA100209), and Project of Ningxia Agricultural Comprehensive Development Office (NTKJ −2013-04; NTKJ −2014-04).

Author information

Authors and Affiliations


Corresponding authors

Correspondence to Shuangyan Chen, Liqin Cheng or Gongshe Liu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

HX, ZL, CL, and CS participated in the design of the study and performed the molecular biology experiments of RNA-seq. PX have made substantial contributions to acquisition of RNA-seq data. HX participated in the design of the non-sequencing experiments and performed the statistical analysis. HX, LG, CS and PX conceived of the study, and helped to draft the manuscript. All authors read and approved the final manuscript.

Xin Huang, Xianjun Peng contributed equally to this work.

Electronic supplementary material

Additional file 1: Clean reads in seven sample libraries.(JPEG 327 KB)

Additional file 2: The quality of the assembly transcripts.(JPEG 683 KB)


Additional file 3: The analysis of sequencing saturation in seven samples. The horizontal axis stands for the number of reads. The vertical axis stands for the contig obtained by contig assembly. (JPEG 315 KB)

Additional file 4: The number of reads that contigs contain in seven samples.(JPEG 718 KB)


Additional file 5: The summary of WEGO output data of sheepgrass genes expressed in response to grazing stress.(DOCX 14 KB)


Additional file 6: The cluster analysis of differently expression genes among the control, defoliation and grazing treatments.(PNG 100 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit

The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Huang, X., Peng, X., Zhang, L. et al. Bovine serum albumin in saliva mediates grazing response in Leymus chinensis revealed by RNA sequencing. BMC Genomics 15, 1126 (2014).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: