Transcriptome Analysis Reveals Common and Distinct Mechanisms for Sheepgrass (Leymus chinensis) Responses to Defoliation Compared to Mechanical Wounding

Background Herbivore grazing is a multiple-component process that includes wounding, defoliation, and saliva deposition. Despite the extensive published research on mechanical wounding and defoliation, no analysis to identify the genes that specify defoliation and mechanical wounding has been performed. Moreover, the influence of the expression of these genes on plant regrowth after defoliation remains poorly understood. Results Seven cDNA libraries for RNA samples collected from stubble tissues that had been mechanically wounded or defoliated at 2, 6 and 24 h along with the control were sequenced using the Illumina/Solexa platform. A comparative transcriptomic analysis of the sequencing data was conducted. In total, 1,836 and 3,238 genes were detected with significant differential expression levels after wounding and defoliation, respectively, during one day. GO, KOG and pathway-based enrichment analyses were performed to determine and further understand the biological functions of those differentially expressed genes (DEGs). The results demonstrated that both wounding and defoliation activated the systemic synthesis of jasmonate (JA). However, defoliation specifically reduced the expression levels of ribosomal protein genes, cell division or cell expansion-related genes, and lignin biosynthesis genes and may have negatively affected plant growth. Further analysis revealed that the regrowth of elongating leaves was significantly retarded after defoliation at 6 h through the following 7 days of measurement, suggesting that the gene expression pattern and phenotype are consistent. Fifteen genes were selected, and their expression levels were confirmed by quantitative RT-PCR (qRT-PCR). Thirteen of them exhibited expression patterns consistent with the digital gene expression (DGE) data. Conclusions These sequencing datasets allowed us to elucidate the common and distinct mechanisms of plant responses to defoliation and wounding. Additionally, the distinct DEGs represent a valuable resource for novel gene discovery that may improve plant resistance to defoliation from various processes.


Introduction
Grasslands are a major part of the global ecosystem, covering 37% of the Earth's terrestrial area [1]. Herbivory is a fundamental driver of grassland plant community composition. Nearly all grasslands worldwide have been grazed by populations of wild and domestic herbivores, such as cattle [2], horses [3], sheep [4], deer [5], rabbits [6], and kangaroos [7]. Herbivore grazing is a multiple-component process that includes wounding, defoliation, and saliva deposition [8].
Wounding responses have been demonstrated in species throughout the plant kingdom, and previous research has focused on identifying systemic wound signals and the mechanisms by which they are generated, transported, and perceived [9]. Leaf wounding is a well-characterized stimulus for the rapid activation of jasmonate (JA) synthesis through distinct autonomous and nonautonomous cellular pathways [9][10][11][12].
Defoliation not only induces serious wounding but also removes a large fraction of photosynthetically active leaves from grass plants. Therefore, plant photosynthesis is reduced; carbohydrate reserves stored in the plant are mobilized; and the transcript levels of photosynthesis-related and carbohydrate metabolism-related genes, such as genes encoding the small subunits of Rubisco, Rubisco activase, light harvesting chlorophyll a/b binding proteins, fructan exohydrolase (FEH), sucrose transporter (SUT), invertase, b-amylase, and starch synthase, are changed [13][14][15][16][17][18].
We previously identified 466 genes that responded to simulated grazing using the Affymetrix Rice GeneChips. These genes encode proteins that participate in signal transduction, miRNA regulation, cell wall modification, carbohydrate metabolism, hormone synthesis, and molecular transport [8]. However, it is unclear whether the responses were due to wounding, defoliation, or other factors.
Currently, the genes that specify defoliation and mechanical wounding remain unknown, and little is understood about how expression of these genes influences plant regrowth after defoliation. In this study, the transcriptomes of sheepgrass stubble tissues collected from defoliated and mechanically wounded plants were compared using the Solexa sequencing system. By comparing the differentially expressed gene (DEG) responses to defoliation and wounding, the common and distinct pathways and genes between defoliation and wounding were identified. These sequencing datasets allowed us to elucidate the common and distinct mechanisms of the plant response to defoliation and wounding, and the distinct DEGs represent a valuable resource for novel gene discovery that may improve plant tolerance toward defoliation caused by herbivores, insects, pathogens, or mechanical stress.

Plant Materials, Growing Conditions, and Treatments
The sheepgrass variety Zhongke 3 was grown in a soil mix of peat moss and vermiculite (2:1, v/v) in a greenhouse at 23uC with a photoperiod of 16 h light/8 h dark. Wounding and defoliation treatments were performed on eight-week-old sheepgrass seedlings. To induce defoliation, approximately two-thirds of the aboveground biomass was removed, and to induce wounding, tweezers were used to mechanically wound a site at the same position as the defoliation site. The remaining one-third of the aboveground biomass was collected 2, 6, and 24 h after defoliation, and the corresponding parts were collected for control seedlings and wounding at 2, 6, and 24 h. Each sample included 3 replicate pots and 21 seedlings. The samples were collected and immediately frozen in liquid nitrogen and stored at 280uC.

RNA Extraction and Solexa/Illumina Sequencing
Total RNA was extracted from individual control, wounded and defoliated plants at 2, 6, and 24 h using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). After precipitation, the RNA was purified using Qiagen's RNeasy Kit with on-column DNase I digestion according to the manufacturer's instructions. RNA concentrations were measured using the Bioanalyzer 2100 (Agilent Technologies, Inc., Waldbronn, Germany), and the RNA integrity was analyzed on a 1.0% (w/v) agarose gel. Approximately 5 to 8 mg total RNA was used to construct each RNA-seq library. Poly (A) mRNAs were isolated from the total RNA using oligo (dT) magnetic beads (Illumina, San Diego, CA). RNA fragmentation, cDNA synthesis, and PCR amplification were performed according to the Illumina RNA-seq protocol (Cat # RS-100-0801). Seven cDNA libraries, which were isolated from mechanically wounded and defoliated material at 2, 6, and 24 h and the control, were sequenced using the Illumina HiSeq 2000 System at the Chinese National Human Genome Center (Shanghai, China). The sequence reads generated in this study have been deposited in the NCBI sequence read archive (SRA065691).

Analysis and Screening of Differentially Expressed Genes (DEGs)
The reads for each sample were counted and mapped back to previously generated sheepgrass reference genes [19] using the software tool Bowtie [20]. The read number of each gene was transformed into RPKM (Reads Per Kilobases per Million reads) [21], and the differently expressed genes were identified by a DEGseq package using the MARS method (MA-plot-based method with Random Sampling model) [22]. The p-value was obtained from the differential gene expression test. FDR (False Discovery Rate) manipulation was used to determine the p-value threshold in multiple tests and analyses. Both FDR #0.001 and the absolute value of | the log2Ratio $1 were used as thresholds to identify the significant DEGs [23].

Gene Ontology, EuKaryotic Orthologous Groups and Pathway Enrichment Analysis
Gene Ontology (GO) is an international standardized functional gene classification system that describes the properties of genes and their products in any organism using three ontologies: cellular component, molecular function, and biological process. GO enrichment analysis was performed using agriGO (http:// bioinfo.cau.edu.cn/agriGO/, [24]). EuKaryotic Orthologous Groups (KOG) enrichment analysis was conducted through hypergeometric distribution testing using the Phyper function in the R software package (http://www.rproject.org/). The Bonferroni correction was used to adjust the p-values. The significantly enriched functional clusters were selected based on a corrected qvalue less than 0.05. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted using KOBAS 2.0 (http://kobas.cbi.pku.edu.cn/, [25]).

Real-time Quantitative RT-PCR (qRT-PCR) Analysis
Real-time quantitative RT-PCR (qRT-PCR) analysis was used to verify the DGE results. The RNA samples used for the qRT-PCR assays were the same as in the DGE experiments. Genespecific primers were designed according to the reference unigene sequences using Primer Premier 6.0. Fifteen genes were selected from the DEGs under defoliation at 2 h for quantitative RT-PCR assays, and these genes were involved in JA biosynthesis, the gibberellin signaling pathway, cell wall biogenesis/degradation, and starch and sucrose metabolism (Table S1). qRT-PCR was performed according to the manufacturer's specifications (Ta-KaRa SYBR Premix Ex Taq II Kit, Dalian, China). The following SYBR Green PCR cycling conditions were used: denaturation at 95uC for 10 s, followed by 40 cycles of 94uC for 5 s and 60uC for 20 s. The PCR experiments were performed using an iQ 5 Multicolor Real-time PCR Detection System (Bio-Rad, USA).
The results were normalized to the sheepgrass actin gene, and the relative gene expression levels were calculated using 2 2DDCT .

Wounding and Defoliation Exhibit Distinct DEGs
To identify the DEGs involved in wounding and defoliation, we sequenced seven cDNA libraries and generated 192 million clean reads ranged from 4,811,093 to 6,411,192 per library after filtering the low-quality and adapter sequences. These clean reads were mapped to the sheepgrass reference genes as described in the Materials and methods. In total, 36 (Table S2). The reference unigenes mapped by all clean reads were enhanced as the sequencing amount (total read number) increased. However, when the sequencing counts reached 1.5 million reads or higher, the number of detected genes became saturated ( Figure S1), indicating that the sequencing depth was sufficient to cover the transcriptome.
By comparing the treated and control Solexa libraries, a great number of differentially expressed transcripts were identified. To study the subset of genes associated with wounding and defoliation, we analyzed the most differentially regulated genes with a |log2Ratio|$1 and a FDR #0.001. There were 706 (356 up-regulated and 350 down-regulated), 891 (442 up-regulated and 449 down-regulated), and 803 (362 up-regulated and 441 down-  Figure 1A). These significantly differentially wounding-and defoliation-responsive genes were further divided into genes expressed at one time point and those expressed at overlapping time points. For example, 1,836 and 3,238 genes were differentially expressed between the wounding and defoliation treatments, respectively, during one day ( Figure 1B, C), and 954 genes were uniquely expressed under only defoliation. These DEGs are listed in Table S3 and Table S4, respectively.
To reveal significantly enriched GO terms in DEGs compared to the reference transcriptome, GO enrichment analysis of the functional significance was performed. The GO terms with p# 0.05 were considered significantly enriched among the DEGs. This analysis allowed us to determine the major molecular functions, biological processes, and cellular components that were represented among the DEGs ( Figure S2). For the molecular functions, the significantly enriched GO terms of the DEGs for wounding and defoliation included DNA binding, serine-type endopeptidase inhibitor activity, oxidoreductase activity, and structural ribosome constituents specific to defoliation. For the enriched biological processes, there were 53 GO terms of DEGs for wounding and 40 GO terms of DEGs for defoliation, which represented cellular component organization and cellular component biogenesis (e.g., chromatin assembly, protein-DNA complex assembly, and nucleosome assembly), metabolic processes (fatty acid and oxylipin biosynthesis), lipid localization, and responses to stimuli (e.g. response to water and defense response to insects). For the cellular components, significantly enriched GO terms of DEGs for wounding and defoliation included chromatin, nucleosome, nucleus, extracellular region, and ribosome, and cell wall was specific to defoliation.

KOG and KEGG Pathway Enrichment Analysis of DEGs
A total of 832 wounding-responsive DEGs and 1,605 defoliation-responsive DEGs were classified functionally into 24 clusters of KOG. KOG enrichment analyses were performed by hypergeometric distribution testing, and the Bonferroni correction was applied. Significantly enriched functional clusters were selected based on a q-value #0.05. The results demonstrated that genes associated with chromatin structure and dynamics, energy production and conversion, and amino acid transport and metabolism were over-represented in the KOG classification of wounding (Table 1), and chromatin structure and dynamics, carbohydrate transport and metabolism, translation, ribosomal structure and biogenesis, and the cytoskeleton were over-represented in the KOG classification of defoliation (Table 2).
A total of 460 wound-responsive DEGs and 881 defoliationresponsive DEGs were classified into the KEGG pathway. For wounding, the significantly enriched KEGG pathways included valine, leucine, and isoleucine degradation; linoleic acid metabolism; ascorbate and aldarate metabolism; starch and sucrose metabolism; fatty acid elongation; arginine and proline metabolism; and alpha-linolenic acid metabolism (Table 3). For defoliation, starch and sucrose metabolism; ribosome; alpha-linolenic acid metabolism; phenylpropanoid biosynthesis; flavonoid biosynthesis; stilbenoid, diarylheptanoid and gingerol biosynthesis; linoleic acid metabolism; and phenylalanine metabolism were over-represented in the KEGG pathway classification ( Table 4).

Enrichment of KEGG Pathways Common to Wounding and Defoliation
From the above KEGG pathway enrichment analysis of the DEGs, we observed that linoleic acid metabolism, alpha-linolenic acid metabolism, and starch and sucrose metabolism were enriched KEGG pathways common to both wounding and defoliation. In the alpha-linolenic acid metabolism and linoleic acid metabolism pathways, we identified genes encoding lipoxygenase (LOX), allene oxide synthase (AOS), allene oxide cyclase (AOC), 12-oxophytodienoic acid reductase, acyl-CoA oxidase (ACOX1, ACOX3), acetyl-CoA acyltransferase 1, and enoyl-CoA hydratase/3-hydroxyacyl-CoA dehydrogenase, all enzymes involved in JA biosynthesis. The gene expression levels were significantly increased for wounding and defoliation at 2 h ( Figure 3). In starch and sucrose metabolism, genes involved in sucrose metabolism, including genes encoding sucrose synthase, fructokinase, and beta-fructofuranosidase (cell wall invertase), were down-regulated, while genes related to starch biosynthesis and metabolism and trehalose biosynthesis, such as genes encoding starch synthase, starch phosphorylase, beta-amylase, and trehalose 6-phosphate synthase, were up-regulated (Figure 3).

Confirmation of Read-mapped Genes by qRT-PCR
To confirm the reliability of the Solexa/Illumina sequencing technology, fifteen genes were selected as described in the Materials and methods for quantitative RT-PCR assays. The expression patterns for thirteen genes were consistent between the qRT-PCR and the DGE analyses ( Figure 5). The corresponding primers are listed in Table S1. The inconsistent expression between the qPCR and DGE analyses for the two remaining genes was most likely because DGE was more sensitive in detecting lowabundance transcripts and small changes in gene expression than qRT-PCR.

Effect of Defoliation on Regrowth of Elongating Leaves
Eight-week-old plants were defoliated to a stubble height of 5 cm. The regrowth of the elongating leaves emerging above the cutting level was measured ( Figure 6). The regrowth of elongating leaves was significantly restrained at 6 h after defoliation compared with the growth of the control plants, and the inhibition continued through the following seven days of measurement.

Distinct Global Gene Transcription Changes for Sheepgrass under Wounding and Defoliation Treatments
In this study, a transcriptome profiling analysis was performed to identify genes that are differentially expressed in wounding-and defoliation-treated plants. A sequencing depth of 4.8 to 6.4 million reads per library was reached (Table S2). The results of the saturation sequencing analyzed in the seven libraries indicated that the sequencing depth was sufficient to cover the transcriptome ( Figure S1). We identified 1,836 and 3,238 wounding-responsive and defoliation-responsive DEGs, respectively, which suggested that defoliation led to nearly twice as many DEGs as wounding. Among these DEGs, 954 were uniquely expressed under defoliation and thus represent candidate genes that merit further study. Further analysis indicated that the defoliation-responsive gene number was highest at 2 h, during which the number of down-regulated defoliation-responsive DEGs (1,526) were over four times the number of down-regulated wounding-responsive DEGs (350), suggesting that plants could rapidly respond to the enormous damage of defoliation by regulating global gene expression.

Wounding and Defoliation Commonly Activate Systemic JA Synthesis
JA is derived from the unsaturated fatty acid a-linolenic acid (18:3), an abundant octadecanoid in the cellular membranes of higher plants. Lipoxygenase (LOX) is the initial step in JA biosynthesis and catalyzes the oxygenation of a-linolenic acid, leading to the formation of hydroperoxy octadecadienoic acids, which are then converted into 12-oxo phytodienoic acid  in the plastids via enzymatic reactions catalyzed first by allene oxide synthase (AOS) and subsequently by allene oxide cyclase (AOC). OPDA is then converted to JA by 12-oxophytodienoic acid reductase (OPR) followed by three cycles of b- oxidation in the peroxisome [26,27]. In this study, we observed that the DEGs involved in JA biosynthesis after wounding and defoliation were significantly enriched according to the GO terms and KEGG pathway enrichment analysis ( Figure S2, Table 3, 4). The transcription levels of key genes in JA biosynthesis, including genes encoding lipoxygenase (LOX), allene oxide synthase (AOS),  allene oxide cyclase (AOC), 12-oxophytodienoic acid reductase (OPR), acyl-CoA oxidase (ACOX1, ACOX3), acetyl-CoA acyltransferase 1, and enoyl-CoA hydratase/3-hydroxyacyl-CoA dehydrogenase (Figure 3), were significantly increased after wounding and defoliation. Our results are consistent with the previous reports, demonstrating that wound-induced responses are triggered by the de novo synthesis of the plant hormone JA when plant tissues are injured by various means, including herbivores, pathogens, and mechanical stress [9,[28][29].
Defoliation Specifically Reduces the Expression Levels of Ribosomal Protein Genes, Cell Division or Cell Expansionrelated Genes, and Lignin Biosynthesis Genes and may Negatively Affect Plant Growth Enrichment analyses of GO terms, KOG, and the KEGG pathway suggested that the ribosome was specifically overrepresented after defoliation. The ribosome is an essential ribonucleoprotein complex that is engaged in translation and is indispensable for growth [30]. In this study, many ribosomal protein genes, including the large 60S and small 40S subunits of the cytosolic ribosome, were identified, and their expression levels mainly decreased at 2 h after defoliation ( Figure 4A). In both plants and animals, the loss of ribosomal proteins leads to reduced growth, which most likely correlated with reduced ribosome production and lower rates of protein synthesis [31]. For example, in Arabidopsis thaliana, several RP loss-of-function mutations that affect cell division or cell expansion and consequently result in a deformed leaf size and shape have been identified, indicating cellor development-specific roles of RPs during leaf growth [32]. Notably, in this study, the expression levels of numerous genes involved in cell division or cell expansion were down-regulated after defoliation, and these genes include those encoding cyclins, cyclin-dependent kinases, and expansion (Table S4).
Lignin is a phenolic heteropolymer of the secondary cell walls that plays a major role in plant development and defense against pathogens. The monolignols represent the main component of lignin, and their synthesis involves many intermediates and enzymes. In this study, we identified several key genes involved in lignin biosynthesis that were down-regulated, including genes encoding phenylalanine ammonia-lyase (PAL), caffeoyl-CoA Omethyltransferase (CCOMT), shikimate O-hydroxycinnamoyltransferase (HCT), and cinnamyl alcohol dehydrogenase (CAD) ( Figure 4B). The down-regulation of these genes has been demonstrated to reduce lignin biosynthesis and cell-wall thickness, alter xylem organization, and retard growth [33][34][35][36].
Notably, the regrowth of elongating leaves was significantly retarded after defoliation at 6 h through the following 7 days of measurement ( Figure 6), suggesting that the above gene expression patterns correlated consistently with the phenotype. Castrillón-Arbeláez et al. demonstrated that although defoliation led to a rapid and transient reduction of non-structural carbohydrates (NSC) in grain amaranth, only a few changes in gene expression and enzyme activity could be associated with the NSC changes in an analysis of 25 genes. They also revealed that the rapid mobilization of foliar starch reserves followed by an efficient recovery of all NSC reserves after defoliation did not alter plant growth or reproductive fitness [17]. Vargas-Ortiz et al. demon-strated that defoliation reduced all NSC levels in the stems and roots of grain amaranth, and this reduction was associated with reduced sucrose synthase and cell wall invertase activity [37]. However, in this study, we observed that the expression levels of genes encoding sucrose synthase and cell wall invertase were reduced under both mechanical wounding and defoliation. Ida et al. also suggested that although defoliation in Oxytropis sericea reduced photosynthesis and nectar production, it did not alter photosynthate allocation or fruit or seed production [38]. However, Akiyama and Å gren demonstrated that defoliation in Arabidopsis thaliana reduced seed production, and this reduction correlated negatively with the removed leaf area [39]. The above contradictory results suggest that defoliation is a complex process, which includes defoliation timing, intensity, and species tolerance, and may lead to different responses. The highly complex nature of the NSC metabolic flux in plants also indicates that defoliation effects cannot be fully explained at the transcriptional level, and its effects on plant growth and reproductive fitness might be elucidated by systems biology in the future.