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

RNA-seq for gene identification and transcript profiling in relation to root growth of bermudagrass (Cynodon dactylon) under salinity stress

BMC Genomics volume 16, Article number: 575 (2015) Cite this article

Abstract

Background

Soil salinity is one of the most significant abiotic stresses affecting plant shoots and roots growth. The adjustment of root architecture to spatio-temporal heterogeneity in salinity is particularly critical for plant growth and survival. Bermudagrass (Cynodon dactylon) is a widely used turf and forage perennial grass with a high degree of salinity tolerance. Salinity appears to stimulate the growth of roots and decrease their mortality in tolerant bermudagrass. To estimate a broad spectrum of genes related to root elongation affected by salt stress and the molecular mechanisms that control the positive response of root architecture to salinity, we analyzed the transcriptome of bermudagrass root tips in response to salinity.

Results

RNA-sequencing was performed in root tips of two bermudagrass genotypes contrasting in salt tolerance. A total of 237,850,130 high quality clean reads were generated and 250,359 transcripts were assembled with an average length of 1115 bp. Totally, 103,324 unigenes obtained with 53,765 unigenes (52 %) successfully annotated in databases. Bioinformatics analysis indicated that major transcription factor (TF) families linked to stress responses and growth regulation (MYB, bHLH, WRKY) were differentially expressed in root tips of bermudagrass under salinity. In addition, genes related to cell wall loosening and stiffening (xyloglucan endotransglucosylase/hydrolases, peroxidases) were identified.

Conclusions

RNA-seq analysis identified candidate genes encoding TFs involved in the regulation of lignin synthesis, reactive oxygen species (ROS) homeostasis controlled by peroxidases, and the regulation of phytohormone signaling that promote cell wall loosening and therefore root growth under salinity.

Background

Soil salinity is one of the most significant abiotic stresses affecting plant shoots and roots growth. Generally, a significant growth reduction in plant shoots is observed under salinity stress [1] as a result of limited photosynthesis [2]. Consequently, most morphological and transcriptional studies on the effects of salinity have majorly focused on shoots and leaves. However, the effect of salinity stress on roots should be more vivid as the root is the organ directly exposed to the salinity [3]. In addition, roots serve numerous key functions, including plant anchourage in the soil, water and nutrients absorption, hormones and metabolites synthesis [4, 5]. Plant root systems display an array of responses to saline conditions, and root architecture alteration is one of such phenotypic responses [6]. These responses are dynamically adjusted in morphology in relation to spatial and temporal salinity heterogeneity to enhance root adaptability to saline conditions [7]. Previous studies on salinity have primarily focused more on the negative effect of salinity than the positive effects on root growth. In most plants, root growth was inhibited under salinity stress [2, 8, 9]. However, for halophytes and a few salt tolerant plants species, root growth was maintained or even promoted under salinity stress. Such responses have been considered as one of the strategies for acclimation to salt stress [10, 11].

Plant roots comprise different regions involved in specific functions for the overall plant maintenance. Root tip is one of the essential regions, which encompasses the cap, the apical meristem, the cell elongation zone and the maturation zone [12]. Root tip is a region of active cell division, elongation. It is also the region of synthesizing vital phytohormones that support continuous growth and development under normal or stressed-conditions [13]. The molecular mechanisms of root growth inhibition under salinity has been extensively investigated [9, 14, 15], while little is known about the molecular mechanisms that control the physiologically well documented processes that are involved in the growth maintenance or promotion of roots in response to salinity stress.

Next generation sequencing (NGS) is a moniker used to represent multiple high throughput or massively parallel nucleic acid sequencing technologies that have emerged since the mid 2000s [16]. In contrast to traditional Sanger sequencing, where typically a few sequential or parallel sequencing reactions of relatively long read lengths (700–1000 bp) generate a modest amount of data, the shared basis of most NGS technologies is the simultaneous execution of millions of sequencing reactions of relatively short read length (30–500 bp) in parallel, and generation of gigabases (Gb) of sequence data per run [16]. The high throughput and low sequencing costs provided by NGS systems has been demonstrated that RNA sequencing is a powerful tool for comparing gene expression, discovering novel transcripts and rare transcripts in plants [17]. RNA-Seq results reveal high levels of reproducibility, for both technical and biological replicates [18]. Currently, NGS is a method of choice to unravel a diversity of stress responses on a transcriptome-wide scale in non-model plant species, where the complete genome sequence and annotation are not yet available. Using two different genotypes contrasting in their salinity stress response facilitated the coverage of a broad spectrum of genes influenced by salt stress, including those involved in a general stress response network, in susceptibility to NaCl and in salt adaptation.

The adjustment of root architecture to spatio-temporal heterogeneity in salinity is particularly crucial for plants growth and survival. Bermudagrass (Cynodon dactylon) is a widely used turf and forage perennial grass, and widely adapted in a number of climatic zones around the world. In addition, bermudagrass has a high degree of salinity tolerance and it is well adapted to salt spray and coastal settings. A moderate salinity appears to stimulate growth of roots and decrease their mortality in this botanical species [11]. This stimulatory effect was found to be genotype-specific. However, information about the genetic and genomic underpinnings associated with these responses is limited. Moreover, molecular mechanisms that control the positive response of root architecture to salinity are still undocumented.

In this study, we performed a transcriptome analysis using RNA-seq to analyze the transcriptome of roots focusing on the molecular basis of the physiological processes during root elongation under salinity stress in two distinct bermudagrass genotypes: C198, a salt-susceptible genotype, and C43, a salt tolerant genotype [11]. The transcriptome analysis identified a complex network of reactive oxygen species (ROS) homeostasis, cell wall loosening and transcription factors involved in the regulation of root tip elongation under salinity. In addition, the comprehensive analysis of bermudagrass root transcriptome can further be utilized as a reference to conduct future research on plant roots.

Results and discussions

Sequencing output and assembly

The bermudagrass salt/control RNA samples described above were used for Illumina Genome Analyzer deep sequencing. An approximate of 0.123 billion (C43) and 0.128 billion (C198) raw reads were generated, and approximately 0.116 billion and 0.121 billion clean reads were generated from C43 and C198 respectively, with a total of 0.25 billion raw reads and 0.24 billion cleans in this sequencing (Table 1). Among all the clean reads, more than 95 % had Phred-like quality scores at the Q20 level (an error probability of 1 %). The four set clean reads of the two genotypes were de novo assembled into one reference transcriptome with the “Trinity” program. After assembly, the four set of clean reads were mapped to the reference transcriptome. Approximately 91 million (C43) and 94 million (C198) read were mapped to the reference transcriptome, which account for 78 % and 77 % of the total clean reads for C43 and C198, respectively (Table 1).

Table 1 Summary of the sequencing results
Full size table

Using the “Trinity” assembler, totally 250,359 transcripts were generated with an average length of 1115 bp and an N50 of 1742 bp, and totally 103,324 unigenes with an average length of 764 bp and an N50 of 1311 bp were obtained by combining C43 and C198 clean reads (Table 2). The length distribution of the transcripts and unigenes from combined C43 and C198 were in the 200–500 bp range making up 35.6 % and 57.5 % of the total, respectively (Table 2).

Table 2 Summary of the assembly results
Full size table

Functional annotation of assembled unigenes

All the assembled high-quality unigenes were first blasted against the National Centre for Biotechnology Information (NCBI) non-redundant (NR) database using BLASTX with a cut-off E-value of 10-5 [19]. Of the 103,324 all unigenes, 45,957 (44.47 %) returned at least one match at the E-value < 10–5. 55.5 % of the unigenes did not match to known genes in the database due to the absence of genome and EST information for C. dactylon or closely related taxa (Table 3). In addition, all the assembled unigenes were annotated by aligning with the other six public databases, including NCBI nucleotide sequences (Nt), Protein family (Pfam), euKaryotic Ortholog Groups (KOG), a manually annotated and reviewed protein sequence database (Swiss-Prot) and Gene Ontology (GO) database with a cut-off E-value of 10−5. Analyses showed that 26,243 unigenes (25.4 % of all unigenes) were annotated with a significant BLAST result in the Nr database; 30,741 unigenes (29.8 % of all unigenes) were annotated in Swiss-Prot database; and 35,960 (34.8 % of all unigenes) unigenes were annotated in the Pfam, and 40,483 (39.2 %) were annotated in GO databases (Table 3). In total, there were 53,765 unigenes (52 %) successfully annotated in at least one of the NR, Nt, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), GO, KOG and Pfam databases, with 2936 unigenes (2.8 %) in all seven databases (Table 3). The majority of transcripts had a significant level of sequence identity to Sorghum bicolor, Zea mays and Oryza sativa proteins, which account for 30.8 %, 22.6 % and 21.1 % of the total transcripts, respectively (Fig. 1).

Table 3 Summary of the annotation results.
Full size table
Fig. 1
figure 1

Species distribution of the top BLAST hits for the bermudagrass sequences

Full size image

Gene ontology (GO) classification

Based on the sequence homology, there were 40,483 unigenes annotated into three ontologies with 47 functional groups (Fig. 2, Additional file 1). Among these groups, genes involved in ‘cellular process’ (23,647), ‘metabolic process’ (22,453) and ‘single-organism process’ (11,836) were highly represented in the biological process (BP) category. The cellular component (CC) category mainly comprised proteins involved in ‘cell’ (17,029), ‘cell part’ (17,002) and ‘organelle’ (13,006). Within the molecular function (MF) category, ‘binding’ (23,866), ‘catalytic activity’ (19,006) and ‘transporter activity’ (2521) were highly represented (Fig. 2).

Fig. 2
figure 2

Histogram of gene ontology (GO) classification. The results are summarized in three main categories: biological process, cellular component and molecular function

Full size image

In addition, all unigenes were subjected to KOG classification for functional prediction. Out of 45,957 nr hits, there were 18,003 unigenes assigned to KOG classification and divided into 26 specific categories (Fig. 3, Additional file 2). Among the 26 KOG categories, the ‘general functional prediction only’ (3110, 17.3 %) was the largest group, followed by post-translational modification, protein turnover, chaperon’ (2400, 13.3 %), ‘translation’ (1915, 10.6 %), ‘signal transduction’ (1664, 9.2 %). The categories ‘extracellular structures (30, 0.17 %) and ‘cell motility’ (18, 0.1 %) had the fewest corresponding genes.

Fig. 3
figure 3

The euKaryotic Ortholog Groups (KOG) annotation of putative proteins. All 18,003 putative proteins showing significant homology to those in KOG database were functionally classified into 26 molecular families

Full size image

To identify further the active biological pathways in C. dactylon, the 45,957 annotated unigenes were mapped to the reference canonical pathways in the KEGG (Fig. 4, Additional file 3). Among those, 12,343 unigenes were assigned to 248 KEGG pathways (Additional file 3). The pathways involving the highest number of unique transcripts were ‘translation’ (1761 unigenes), followed by ‘carbohydrate metabolism’ (999 unigenes), ‘signal transduction’ (978 unigenes), ‘folding, sorting and degradation’ (963 unigenes) and ‘energy metabolism’ (869 unigenes), indicating that these were the active pathways in C. dactylon.

Fig. 4
figure 4

Pathway assignment based on Kyoto Encyclopedia of Genes and Genomes (KEGG). a Cellular Processes; b Environmental Information Processing; c Genetic Information Processing; d Metabolism

Full size image

Differential expression analysis of assembled transcripts

We used the normalized RPKM (reads per kilobase per million) to quantify the transcript level in reads, which facilitated the comparison of mRNA levels both within and between genotypes. Differential expressed genes (DEGs) (q-value < 0.005 and |log2 (fold change)| >1) were defined as genes that were significantly enriched or depleted in one genotype and/or treatment relative to the other genotype and/or treatment (Additional file 4). A hierarchical clustering of the differentially expressed genes based on the four sample’s log10RPKM was made, so we could observe the overall gene expression pattern. The blue bands identify low quantity gene expression, while the red represent the high quantity gene expression (Fig. 5). Four groups of DEGs with specific expression patterns were delineated from the clustering. DEGs showed down-regulation in Group A and up-regulation in group C under salinity for both genotypes when compared to the control level. DEGs in Group B and group D showed no changes in expression levels between control and salinity in both genotypes with contrary expression pattern between the two genotypes (Fig. 5).

Fig. 5
figure 5

Hierarchical clustering analysis of salinity-induced changes in gene expression in root tips of bermudagrass (C43_salt, 200 mM NaCl treated C43; C198_salt, 200 mM NaCl treated C198; C43_control, 0 mM NaCl treated C43; C198_control, 0 mM NaCl treated C198)

Full size image

Comparison of changes in gene expression between control and salt-stressed plants in two genotypes had shown similarities and considerable differences. The Venn diagram indicates the distribution of expressed genes among the four samples. Under control (0 mM NaCl) and salt stress (200 mM NaCl), there were 277 and 314 DEGs between C43 and C198, respectively. Under salt-stressed conditions, there were 848 and 536 DEGs between the control and the salt-stressed roots for C43 and C198, respectively (Fig. 6). There were more genes exclusively differentially expressed in C43 (383) than that of C198 (95) under salinity.

Fig. 6
figure 6

Comparison between the amounts of differential expressed genes (DEGs) found in the two genotypes by mapping reads to the de novo assembled TCs. A Venn diagram depicting the number of statistically significant (>2-fold) DEGs when the de novo transcriptome (green and red ovals, FDR <0.025). The numbers of DEGs exclusively expressed in one sample are shown in each circle. The numbers of DEGs with a common or opposite tendency of expression changes between the two treatments are shown in the overlapping regions. The total numbers of DEGs in each treatment are shown outside the circles (C43_salt, 200 mM NaCl treated C43; C198_salt, 200 mM NaCl treated C198; C43_control, 0 mM NaCl treated C43; C198_control, 0 mM NaCl treated C198)

Full size image

Validate the DEGs by real-time RT-PCR analysis

To validate the data from RNA-sequencing, 46 DEGs were randomly selected for real-time RT-PCR analysis in both genotypes in response to salt stress. The primers of selected genes are listed in Additional file 5. Actin was used as reference gene for data normalization according to Hu et al. [11]. The qRT-PCR results showed a strong correlation with the RNA-seq-generated data (Pearson correlation coefficients r = 0.87; Fig. 7, Additional file 5).

Fig. 7
figure 7

qRT-PCR validations of differentially expressed genes in root of bermudagrass under salinity. Correlation of fold change analyzed by RNA-Seq platform (x axis) with data obtained using real-time PCR (y axis)

Full size image

Transcription factors in relation to salinity stress and root growth regulation

The assembled transcriptome of salt-stressed root tips demonstrated that totally 24 unique TFs were differentially expressed under salinity stress in both genotypes, with 17 exclusively in C43 (Table 4, Additional file 6). They included major TF families linked to stress responses and growth regulation, such as AP2/ERF, bZIP, MYB, MYC, NAC, MADS-box, WRKY, bHLH, zinc finger family (Table 4). Several salt-induced TF genes also respond to other abiotic stresses such as osmotic stress, cold and heat [15], suggesting that they generally participate in stress response, and the spatial differences of TF gene regulation by environmental stresses in root tips may be crucial for the adaptation of their growth to specific soil environments. TF genes encoding AP2/ERF, MYB, MYB, NAC, WRKY showed a considerable enhancement of expression in root apexes in response to salt stress when compared to whole roots in M. truncatula [15]. The regulation of specific members of TF families in Medicago root tips supports the hypothesis that these genes may intersect root developmental pathways and salt-related transcriptional networks [15].

Table 4 Transcription factors differentially expressed in the two genotypes under salinity stress
Full size table

Among the differentially expressed TF genes in root tip, there were several classes (MYB, NAC and bHLH) involved in growth regulation, and other metabolic processes. The differential expression of many of these genes in the root tip is not striking as they participate in active processes in the root tip such as expansion and cell division, as well as the encountering a variety of biotic/abiotic stresses [19]. It has been previously demonstrated that bHLH-type TFs are linked to the adaptation of Medicago truncatula to saline soil environments [20], and a group of bHLH-type TFs were found to be involved in the root growth and development of Medicago sativa under salinity [15]. A bHLH-type TF (comp117508_c0) was up-regulated, while a number of peroxidases were down-regulated in C43 genotype, this observation indicated that the bHLH TF might modulate the ROS balance by directly regulating the expression of a set of peroxidases, consequently, regulating the root cell proliferation and differentiation [21].

The MYB TFs is one of the largest superfamily of plant TFs. Apart from metabolic, signal transduction and defense-related pathways, the induced or repressed expressions of TFs affected plant growth and development under salinity [22, 23]. The over expression of MYB31 and MYB42 repressed the expression of genes related to lignin synthesis in maize [24], which exhibited up to 45 % reduction in lignin content and substantially increased leaf, root, and stem growth [25]. A MYB-like TF (KUA1) modulates leaf cell expansion and final organ size by controlling the expression of peroxidases and ROS homeostasis [26]. Three MYB TFs were differentially expressed in C43 with two down-regulated (comp134391_c0, comp112426_c0) and one MYB TF (comp120905_c3) up-regulated under salinity. This observation suggests that the induction or repression of MYB TFs may be participating in the lignin synthesis and/or ROS homeostasis controlled by peroxidases.

The WRKY transcription factors are considered to be repressors of the GA signaling pathway [27], activators of the ABA signaling pathway [28] and regulators of many other signaling pathways in plants [29]. WRKY31 gene was found to enhance disease resistance and reduced lateral root formation and elongation with induced constitutive expression of auxin-response genes, such as OsIAA4 and OsCrl1 genes [30]. Two WRKY TFs (comp115018_c1 and comp120384_c0) were down-regulated and three WRKY TFs (comp112466_c0, comp122277_c0 and comp117715_c0) up-regulated in C43. One WRKY TFs (comp123500_c0) were up-regulated in C198 (Additional file 6). These differentially expressed WRKY TFs may be involved in the regulation of root growth under salinity through the regulation of phytohormone signaling.

The expression pattern of genes involved in cell wall loosening

In contrast to mammalian cells, plant cells are encased by a cell wall that gives structural support, and cell expansion is affected by alterations in cell wall and architecture [31]. Several proteins have been directly implicated in cell wall loosening, including xyloglucan endotransglucosylase/hydrolases (XETs) [32]. In addition, non-enzymatic processes involving ROS that produce wall polysaccharide scission also participate in cell wall loosening [33]. Cell walls as conditioners of cell growth under salt stress have been investigated in detail [34]. The balance between cell wall loosening and stiffening activities defines the regions of accelerated and decelerated root growth in the elongation zone [35]. The root system comprises different regions that are involved in specific functions for the overall plant sustenance. Root tips are essential regions that encompass the root apical meristem and elongation zone. They produce pall types of cells in a highly defined pattern of cell division to assist growth during saline conditions [36].

Root tip transcriptome analyses indicate that the expression of genes related to cell wall loosening and stiffening is modified by salinity in bermudagrass. One XET transcript (comp116615_c0) was found to be up-regulated in C43 (Additional file 7). The XET transcript expression pattern following the distribution of growth rates in the growing zone of Festuca pratensis, and the resulting XET activity was proposed to be involved in cell wall modification processes during cell elongation [37]. XET activity was enhanced in the apical region of maize roots from plants grown under low water potential [38]. These results indicated that the increased expression level of XET transcript in the root tips of bermudagrass under salinity might be necessary for maintaining root elongation under these conditions.

Salinity may promote cell wall stiffening possibly mediated by peroxidase [39]. Peroxidases are considered bifunctional enzymes that not only oxidize various substrates in the presence of H2O2, but also generate H2O2 [40]. The relationship between plant peroxidase and cell expansion restriction has been shown in several studies [41]. Apoplastic peroxidases are known to either restrict or promote cell expansion [21]. In C. gayana leaves, the peroxidase activity increase, and the phenolic compounds of the cell walls caused a reduction in the length of the elongation zone grown under saline conditions [42]. However, Cramer et al. [43] did not observe increased wall stiffening or stimulated peroxidase activity under salinity. This result suggests that the apparent contradictory effects are linked to the regulatory modes under which the peroxidases are working. The repressed expression of peroxidases by KUA1 (a MYB-like TFs) in leaves promoted cell expansion, which is clearly linked to changed levels of apoplastic H2O2 [26]. In our study, 21 peroxidases transcript were found to be differentially expressed with 19 down-regulated in C43 and 17 down-regulated in C198. The down-regulation of those peroxidases transcript may favor root elongation by the reduction of apoplastic H2O2, and the generation of oxygen radicals to cleave the cell wall polymers thus promote cell wall loosening and therefore root growth [26].

Conclusions

Root growth maintenance or promotion under salinity is a complex process that integrates spatio-temporal developmental events from the sensing of osmotic and ionic stress. We performed RNA-seq in the root tips of two genotypes of bermudagrass, which had different root growth characters and salt tolerance. The aim of the sequencing was to identify the causes for the observed changes in the spatial pattern of root elongation under salinity. In total, 250,359 transcripts with an average length of 1115 bp and totally 103,324 unigenes obtained with 53,765 unigenes (52 %) were successfully annotated in databases. Moreover, the 848 and 536 differentially expressed unigenes in C43 and C198 were useful to identify the genes related to the root growth regulation under salinity. The RNA-seq identified candidate genes encoding TFs involved in the regulation of lignin synthesis, ROS homeostasis controlled by peroxidases, and the regulation of phytohormone signaling that promote cell wall loosening and hence root growth under salinity. Our results support and add detailed molecular information to the root growth maintenance under stress via increased wall loosening. The data suggest need for control of intracellular ROS content by peroxidase under the regulation of transcription factors for apoplastic hydrogen peroxide production and other cell wall proteins for wall loosening induction. Collectively, the data indicate that the regulation of root growth under salinity involves changes in many different aspects of cell metabolism, signaling, and transport.

Methods

Plant growth, salt treatment and sampling

Uniform stolons (5 cm long) of bermudagrass, ‘C43’ (salt tolerant) and ‘C198’ (salt sensitive) were planted in solid growth substances (peat soil: sand = 2:1, v/v) [11]. After 14 d of establishment, equal amount of plantlets were transfer to plastic pots (7 cm diameter and 9 cm height) filled with coarse silica sand as the plant anchor medium. Pots were suspended over tubs containing 46 L of constantly aerated half-strength Hoagland’s solution [44]. The tubs were refilled every other day and renewed weekly. Pot bottoms consisted of a coarse nylon screen allowing roots to freely grow into the solutions. Plants were grown in an environmentally controlled walk-in growth room with the temperature regime of 30 /25 °C (day/night), photosynthetically active radiation (PAR) levels of 800 μmol · m−2 · s−1 at canopy height for 14 h. Plants were allowed to adapt to this nutrient solution for 2 weeks. During this period, the plant shoots were hand-clipped weekly at 6 cm height, and roots were clipped back to the bottoms of the pots at the beginning of the salt treatment in order to allow the plants to reach full maturity and develop uniform and equal size roots and shoots.

After 2 weeks of hydroculture, the half-strength Hoagland’s solution was supplemented with 200 mM NaCl. After 7 d of treatment, plants were removed from the nutrient solution and gently washed for 20 s with distilled water for the roots. Roots elongated from nylon screen into nutrient solution were sampled after 4 h light. The root tips, encompassing the meristem and the elongation zone, were excised with a scalpel from the remaining root system, and immediately frozen in liquid nitrogen. Multiple independent biological replicates, each containing a pool of twenty different plants, were sampled for mRNA-Seq (two biological replicates) and reverse transcription-PCR (RT-PCR) or quantitative RT-PCR (qRT-PCR) validation (three biological replicates).

The salt treatments and grass genotypes were arranged in a randomized complete block design with six replicates. Physiological measuring to one set of plants and taking tissue samples for RNA isolation and metabolite analysis to another set of plants were performed simultaneously at the time period of 1000–1500 h.

Total RNA, mRNA isolation and library preparation for transcriptome sequencing

Total RNA was extracted using Trizol reagent (Invitrogen) and purified using the RNeasy Plant Mini kit (Qiagen) according to the Handbook, with an additional sonication step after addition of RLT buffer (Qiagen). The quality and integrity of RNA was checked by Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) and agarose gel electrophoresis. A total amount of 3 μg RNA per sample from root tips was used for mRNA-Seq library construction using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendations. The mRNA was isolated using oligo(dT)-attached magnetic beads and subsequently fragmented using divalent cations under elevated temperature, and the cleaved RNA fragments were copied into first-strand cDNA using random hexamer primer and M-MuLV Reverse Transcriptase. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H, and the cDNA fragment were processed for end repair, an addition of a single “A” base, and ligation of the adapters. After second-strand cDNA synthesis and adaptor ligation, cDNA fragments of 150 ~ 200 bp in length were isolated with AMPure XP system (Beckman Coulter, Beverly, USA). These products were then purified and enriched by PCR to create the final cDNA library. After cluster generation on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia), the multiplexed library was sequenced on an Illumina Hiseq 2000 platform according to the manufacturer’s recommendations (Illumina) at Novogene Bioinformatics Institute, Beijing, China. RNA-seq read data were deposited to the NCBI Sequence Read Archive (NCBI SRA) under accession number SRR2075766.

Preprocessing and de novo assembly

In total, 25.18 G bases (Gb) raw reads were generated by Illumina Hi-seq platform, with a total of ~ 5.77 Gb, ~6.05, ~ 6.58 Gb and ~ 6.78 Gb in C43-control, C198-control, C43-salt and C198-salt, respectively. The raw reads were initially processed through in-house perl scripts to remove the adapter sequences, reads containing ploy-N, and low-quality bases, and finally get the clean data (clean reads). All the downstream analyses were based on clean data with high quality. After preprocessing, we obtained four set clean reads, with a total of ~ 5.44 Gb, ~ 5.72 Gb, ~6.22 Gb and ~6.4 Gb quality filtered short reads for C43-control, C198-control, C43-salt and C198-salt, respectively. The four set clean reads of the two species were de novo assembled with the “Trinity” program (v2012-10-05) after merged following the protocol documented in [45] with the min_kmer_cov set to 2 and all other parameters set default. The base calling and base quality assignment were evaluated by using PHRED [46].

Unigene annotation and classification

Unigenes (≥100 bp) were used to search against the NR [19], SwissProt [47], KEGG (version 58) [48] and KOG [49] databases by BLASTALL package (release 2.2.22) with the significant threshold of E-value ≤ 10−5. Each known gene from the best BLASTx hit was parsed and assigned. The ORF of assembled transcripts was determined based on the results of BLASTx search in the following order: NR, KEGG and KOG. Extending from both sides of the aligned region, the coding region sequences were translated into amino acid sequences with the standard codon table using custom PERL scripts. For unigenes that did not align to any of the above databases were scanned by ESTScan to determine the nucleotide and amino acid sequences of the coding regions. The predicted amino sequences were submitted to search against the Pfam database (version 26.0) [50] for domain/family annotation using HMMER 3.0, with the ‘Best Match Cascade’ protocol. Gene ontology (GO) [51] terms for each unigenes were assigned based on the best BLASTx hit from the NR and Pfam database using Blast2GO software (version 2.5) [52] with an E-value threshold of 10−5.

Quantification of gene expression levels and differential expression analysis

Gene expression levels were estimated by RSEM [53] for each sample. Clean data were mapped back onto the assembled transcriptome. Readcount for each gene was obtained from the mapping results and normalized to reads per kb of exon model per million mapped reads (RPKM) [53]. Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor [54]. Differential expression analysis of two samples was performed using the DEGseq (2010) R package. P-value was adjusted using q value. q value < 0.005 and the absolute value of log2 > 1 was set as the threshold for significant difference in gene expression. Gene ontology (GO) analysis of the differentially expressed genes (DEGs) was done for biological process, cellular components and molecular function by BGI WEGO (Web Gene Ontology Annotation Plotting, http://wego.genomics.org.cn/cgi-bin/wego/index.pl) [55] and agriGO (GO Analysis Toolkit and Database for Agricultural Community, http://bioinfo.cau.edu.cn/agriGO/index.php) [56], and pathways that were statistically significantly (Q-value < 0.05) were enriched with KEGG [57]. The differential expression profiles of the genes have been colored based on the Venn diagram [58].

Quantitative real-time PCR analysis

The expression of selected candidate genes was validated by quantitative real time PCR (qRT-PCR) using the same RNA samples as in the RNA-seq library construction. The first strand cDNA fragments were synthesized from 2 μg of total RNA using oligo(dT)12–18 primer using cDNA synthesis kit (Fermentas, Burlington, Ontario, Canada) according to the user manual. Gene-specific primers (Table 1) were designed based on the target gene sequences using Primer 5 software. UBI gene was used as internal standard. The qRT-PCRs were performed with ABI7500 in a final volume of 20 μL, with each containing 2 μL of cDNA, 10 μL of 2 × SYBR Green qPCR Mix (Takara, Otsu, Shiga, Japan) and 2 μM of the forward and reverse primers. Three independent biological replicates of each sample and two technical replicates of each biological replicate were used for real-time PCR analysis. The thermal cycling conditions were as follows: 40 cycles of 95 °C denaturation for 5 s, and 52 ~ 55 °C annealing and extension for 20 s. After the PCR, a melting curve was generated by gradually increasing the temperature to 95 °C to test the amplicon specificity. To determine relative fold differences for each sample, the CT value for each gene was normalized to the CT value for the reference gene and was calculated relative to a calibrator using the DDCT method as described by Livak and Schmittgen [59].

Abbreviations

BLAST:

Basic local alignment search tool

bHLH:

Basic Helix-Loop-Helix

bZIP:

Basic leucine zipper

BP:

Biological process

CC:

Cellular component

DEGs:

Differential expressed genes

GO:

Gene ontology

KEGG:

The Kyoto encyclopedia of genes and genomes database

KOG:

euKaryotic Ortholog Groups

MF:

Molecular function

NAC:

NAM/ATAF/CUC domain protein

NCBI:

National center for biotechnology information

NGS:

Next-generation sequencing

NR:

NCBI non-redundant protein database

Nt:

Nucleotide sequences database

PAR:

Photosynthetically active radiation

PERL:

Practical extraction and report language

Pfam:

Protein family, Swiss-Prot, a manually annotated and reviewed protein sequence database

q-PCR:

Real-time quantitative PCR

RNA-Seq:

RNA-sequencing

ROS:

Reactive oxygen species

RPKM:

Reads per kb per million reads

TF:

Transcription factor

XET:

Xyloglucan endotransglucosylase

References

  1. Katsuhara M, Kawasaki T. Salt stress induced nuclear and DNA degradation in meristematic cells of barley roots. Plant Cell Physiol. 1996;37(2):169–73.

    Article  CAS  Google Scholar 

  2. Geilfus CM, Zörb C, Mühling KH. Salt stress differentially affects growth-mediating β-expansin in resistant and sensitive maize (Zea mays L.) cultivars. Plant Physiol Biochem. 2010;48(12):993–8.

    Article  CAS  PubMed  Google Scholar 

  3. Hajibagheri MA, Yeo AR, Flowers TJ. Salt tolerance in Suaeda maritima (L.) Dum. Fine structure and ion concentrations in the apical region of roots. New Phytol. 1985;99(3):331–43.

    Article  Google Scholar 

  4. Waisel Y, Eshel A, Kafkafi U. Plant roots-the hidden half. 3rd ed. New York, NY, USA and Basel, Switzerland: Marcel Dekker, Inc.; 2002.

    Google Scholar 

  5. Zhai R, Feng Y, Wang H, Zhan X, Shen X, Wu W, et al. Transcriptome analysis of rice root heterosis by RNA-Seq. BMC Genomics. 2013;14:19.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  6. Potters G, Pasternak TP, Guisez Y, Palme KJ, Jansen MAK. Stress-induced morphogenic responses, growing out of trouble? Trends Plant Sci. 2007;12(3):98–105.

    Article  CAS  PubMed  Google Scholar 

  7. Lynch JP, Brown KM. Topsoil foraging–an architectural adaptation of plants to low phosphorus availability. Plant Soil. 2001;237(2):225–37.

    Article  CAS  Google Scholar 

  8. Zolla G, Heimer YM, Barak S. Mild salinity stimulates a stress-induced morphogenic response in Arabidopsis thaliana roots. J Exp Bot. 2010;61:211–24.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  9. Cotsaftis O, Plett D, Johnson AAT, Walia H, Wilson C, Ismail AM, et al. Root-specific transcript profiling of contrasting rice genotypes in response to salinity stress. Mol Plant. 2011;4(1):25–41.

    Article  CAS  PubMed  Google Scholar 

  10. Pessarakli M, Touchane H. Growth responses of bermudagrass and seashore paspalum under various levels of sodium chloride stress. J Food Agric Environ. 2006;4(3&4):240–3.

    Google Scholar 

  11. Hu LX, Huang ZH, Liu SQ, Fu JM. Growth response and gene expression in antioxidant-related enzymes in two bermudagrass genotypes differing in salt tolerance. J Amer Soc Horticul Sci. 2012;137(3):134–43.

    CAS  Google Scholar 

  12. Dolan L, Janmaat K, Willemsen V, Linstead PJ, Poethig RS, Roberts K, et al. Cellular organisation of the Arabidopsis root. Development. 1993;119(1):71–84.

    CAS  PubMed  Google Scholar 

  13. Bielach A, Podlesáková K, Marhavý P, Duclercq J, Cuesta C, Muller B, et al. Spatiotemporal regulation of lateral root organogenesis in Arabidopsis by cytokinin. Plant Cell. 2012;24(10):3967–81.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  14. Neves-Piestun BG, Bernstein N. Salinity-induced changes in the nutritional status of expanding cells may impact leaf growth inhibition in maize. Funct Plant Biol. 2005;32(2):141–52.

    Article  CAS  Google Scholar 

  15. Postnikova OA, Shao J, Nemchinov LG. Analysis of the alfalfa root transcriptome in response to salinity stress. Plant Cell Physiol. 2013;54(7):1041–55.

    Article  CAS  PubMed  Google Scholar 

  16. Shendure J, Ji H. Next-generation DNA sequencing. Nat Biotechnol. 2008;26:1135–45.

    Article  CAS  PubMed  Google Scholar 

  17. Zhang N, Liu B, Ma C, Zhang G, Chang J, Si H, et al. Transcriptome characterization and sequencing-based identification of drought-responsive genes in potato. Mol Biol Rep. 2014;41(1):505–17.

    Article  CAS  PubMed  Google Scholar 

  18. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, et al. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science. 2008;320(5881):1344–9.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  19. Pereira SS, Guimarães FCM, Carvalho JFC, Stolf-Moreira R, Oliveira MCN, Rolla AAP, et al. Transcription factors expressed in soybean roots under drought stress. Genet Mol Res. 2011;10(4):3689–701.

    Article  CAS  PubMed  Google Scholar 

  20. Zahaf O, Blanchet S, De Zelicourt A, Alunni B, Plet J, Laffont C, et al. Comparative transcriptomic analysis of salt adaptation in roots of contrasting Medicago truncatula genotypes. Mol Plant. 2012;5(5):1068–81.

    Article  CAS  PubMed  Google Scholar 

  21. Tsukagoshi H, Busch W, Benfey PN. Transcriptional regulation of ROS controls transition from proliferation to differentiation in the root. Cell. 2010;143(4):606–16.

    Article  CAS  PubMed  Google Scholar 

  22. Cominelli E, Tonelli C. A new role for plant R2R3-MYB transcription factors in cell cycle regulation. Cell Res. 2009;19:1231–2.

    Article  PubMed  Google Scholar 

  23. Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010;15(10):573–81.

    Article  CAS  PubMed  Google Scholar 

  24. Fornalé S, Sonbol FM, Maes T, Capellades M, Puigdomènech P, Rigau J, et al. Down-regulation of the maize and Arabidopsis thaliana caffeic acid O-methyl-transferase genes by two new maize R2R3-MYB transcription factors. Plant Mol Biol. 2006;62(6):809–23.

    Article  PubMed  Google Scholar 

  25. Hu WJ, Harding SA, Popko JL, Ralph J, Stokke DD, Tsai CJ, et al. Repression of lignin biosynthesis promotes cellulose accumulation and growth in transgenic trees. Nature Biotech. 1999;17:808–12.

    Article  CAS  Google Scholar 

  26. Lu D, Wang T, Persson S, Mueller-Roeber B, Schippers JH. Transcriptional control of ROS homeostasis by KUODA1 regulates cell expansion during leaf development. Nat Commun. 2014;5:3767.

    CAS  PubMed Central  PubMed  Google Scholar 

  27. Zhang ZL, Xie Z, Zou X, Casaretto J, Ho TH, Shen QJ. A rice WRKY gene encodes a transcriptional repressor of the gibberellin signaling pathway in aleurone cells. Plant Physiol. 2004;134(4):1500–13.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  28. Zou X, Seemann JR, Neuman D, Shen QJ. A WRKY gene from creosote bush encodes an activator of the abscisic acid signaling pathway. J Biol Chem. 2004;279:55770–9.

    Article  CAS  PubMed  Google Scholar 

  29. Li J, Brader G, Palva ET. The WRKY70 transcription factor: a node of convergence for jasmonate-mediated and salicylate-mediated signals in plant defense. Plant Cell. 2004;16(2):319–31.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  30. Zhang J, Peng YL, Guo ZJ. Constitutive expression of pathogen-inducible OsWRKY31 enhances disease resistance and affects root growth and auxin response in transgenic rice plants. Cell Res. 2008;18(4):508–21.

    Article  CAS  PubMed  Google Scholar 

  31. Rubio-Díaz S, Pérez-Pérez JM, González-Bayón R, Muñoz-Viana R, Borrega N, Mouille G, et al. Cell expansion-mediated organ growth is affected by mutations in three EXIGUA genes. PLoS ONE. 2010;7(5):e36500.

    Article  Google Scholar 

  32. Nishitani K. Division of roles among members of the XTH gene family in plants. Plant Biosyst. 2005;139(1):98–101.

    Article  Google Scholar 

  33. Schweikert C, Liszkay A, Schopfer P. Scission of polysaccharides by peroxidase-generated hydroxyl radicals. Phytochemistry. 2000;53(5):565–70.

    Article  CAS  PubMed  Google Scholar 

  34. Zhong H, Läuchli A. Changes of cell wall composition and polymer size in primary roots of cotton seedlings under high salinity. J Exp Bot. 1993;44(4):773–8.

    Article  CAS  Google Scholar 

  35. Muller B, Reymond M, Tardieu F. The elongation rate at the base of a maize leaf shows an invariant pattern during both the steady-stage elongation and the establishment of the elongation zone. J Exp Bot. 2001;52(359):1259–68.

    Article  CAS  PubMed  Google Scholar 

  36. Brady SM, Orlando DA, Lee JY, Wang JY, Koch J, Dinneny JR, et al. A high-resolution root spatiotemporal map reveals dominant expression patterns. Science. 2007;318(5851):801–6.

    Article  CAS  PubMed  Google Scholar 

  37. Reidy B, Nösberger J, Fleming A. Differential expression of XET-related genes in the leaf elongation zone of F. pratensis. J Exp Bot. 2001;52(362):1847–56.

    Article  CAS  PubMed  Google Scholar 

  38. Wu Y, Jeong BR, Fry SC, Boyer JS. Change in XET activities, cell wall extensibility and hypocotyl elongation of soybean seedlings at low water potential. Planta. 2005;220(4):593–601.

    Article  CAS  PubMed  Google Scholar 

  39. Neumann PM, Azaizeh H, Leon D. Hardening of root cell walls: a growth inhibitory response to salinity stress. Plant Cell Environ. 1994;17(3):303–9.

    Article  Google Scholar 

  40. Passardi F, Pennel C, Dunand C. Performing the paradoxical: how plant peroxidases modify the cell wall. Trends Plant Sci. 2004;9(11):534–40.

    Article  CAS  PubMed  Google Scholar 

  41. Fry SC, Willis SC, Paterson AE. Intraprotoplasmic and wall-localised formation of arabinoxylan-bound diferulates and larger ferulate coupling-products in maize cell-suspension cultures. Planta. 2000;211(5):679–92.

    Article  CAS  PubMed  Google Scholar 

  42. Ortega L, Fry SC, Taleisnik E. Why are Chloris gayana leaves shorter in salt-affected plants? Analyses in the elongation zone. J Exp Bot. 2006;57(14):3945–52.

    Article  CAS  PubMed  Google Scholar 

  43. Cramer G, Schimdt C, Bidart C. Analysis of cell wall hardening and cell wall enzymes of salt-stressed maize (Zea mays) leaves. Funct Plant Biol. 2001;28(2):101–9.

    Article  CAS  Google Scholar 

  44. Hoagland DR, Arnon DI. The water-culture method for growing plants without soil. California Agri Exp Stat Cir. 1950;347(2):1–32.

    Google Scholar 

  45. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  46. Ewing B, Green P. Base-calling of automated sequencer traces using Phred. II Error probabilities Genome Res. 1998;8(3):186–94.

    CAS  PubMed  Google Scholar 

  47. Wu CH, Apweiler R, Bairoch A, Natale DA, Barker WC, Boeckmann B, et al. The universal protein resource (UniProt): an expanding universe of protein information. Nucleic Acids Res. 2006;34(Database issue):D187–91.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  48. Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The KEGG resource for deciphering the genome. Nucleic Acids Res. 2004;32(Database issue):D277–80.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  49. Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR, Kiryutin B, Koonin EV, et al. The COG database: an updated version includes eukaryotes. BMC Bioinformatics. 2003;4:41.

    Article  PubMed Central  PubMed  Google Scholar 

  50. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40(Database issue):D290–301.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  51. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  52. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the blast2go suite. Nucleic Acids Res. 2008;36(10):3420–35.

    Article  PubMed Central  PubMed  Google Scholar 

  53. Li B, Dewey CN. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  54. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  55. Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, et al. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;4(Web Server issue):W293–7.

    Article  Google Scholar 

  56. Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38 (Web Server issue): W64–70.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  57. Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res 2012;40(Databaseissue):D109-12.

  58. Pirooznia M, Nagarajan V, Deng Y. GeneVenn-A web application for comparing gene lists using Venn diagrams. Bioinformation. 2007;1(10):420–2.

    Article  PubMed Central  PubMed  Google Scholar 

  59. Livaka KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2 − ΔΔCT method. Methods. 2001;25(4):402–8.

    Article  Google Scholar 

Download references

Acknowledgement

This work was financially supported by the Young Scientists Fund of the National Natural Science Foundation of China (Grant No. 31101563), China-Africa Center for Research and Education, CAS (SAJC201325) and National Natural Science Foundation of China (Grant No. 31272194).

Author information

Authors and Affiliations

  1. Key Laboratory of Plant Germplasm Enhancement and Specialty Agriculture, Wuhan Botanical Garden, The Chinese Academy of Sciences, Wuhan, Hubei, 430074, PR China

    Longxing Hu, Huiying Li, Liang Chen, Yanhong Lou, Erick Amombo & Jinmin Fu

Authors
  1. Longxing Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Huiying Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Liang Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yanhong Lou
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Erick Amombo
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Jinmin Fu
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Jinmin Fu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

HL performed the salinity experiment, prepared the mRNA for sequencing and participated in the writing of the manuscript. LH analyzed the data. AE and CL revised the manuscript. LY provided the plants and contributed to the sample collection. FJ designed the project and revised the manuscript. All authors read and approved the final manuscript.

Additional files

Additional file 1:

GO_classification. (XLS 10826 kb)

Additional file 2:

KOG_classification. (XLS 270 kb)

Additional file 3:

KEGG_classification. (XLS 264 kb)

Additional file 4:

Differentially expressed genes in two genotypes of bermudagrass. (XLS 1612 kb)

Additional file 5:

Pimers for qRT-PCR analysis. (XLS 31 kb)

Additional file 6:

Selected differentially expressed transcription factors. (XLS 20 kb)

Additional file 7:

Selected cell wall related DEGs. (XLS 25 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hu, L., Li, H., Chen, L. et al. RNA-seq for gene identification and transcript profiling in relation to root growth of bermudagrass (Cynodon dactylon) under salinity stress. BMC Genomics 16, 575 (2015). https://doi.org/10.1186/s12864-015-1799-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-015-1799-3

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us