Identification of drought stress-responsive transcription factors in ramie (Boehmeria nivea L. Gaud)

Stem traits and fiber yield of ramie in response to drought stress

Under well water condition, the stem length, diameter and bark thickness were 128.9 cm, 11.79 mm and 0.987 mm, respectively (Figure 1); whereas significant decreases in stem length, diameter and bark thickness (98.6 cm, 9.70 mm and 0.793 mm, respectively) were observed when ramie suffered from drought stress (Figure 1). Besides, the fiber yield of drought-stressed ramie (6.62 g per plant) was far lower than that of well-watered ramie (8.99 g per plant) (Figure 1). Therefore, drought environment severely inhibits the stem growth of ramie and leads to a decrease of the fiber yield.

Tag identification and quantification

Depth of sampling

Saturation of the library is determined by checking the number of detected genes. Sequencing reaches saturation when no new genes are detected. The results showed that when sequencing amount reached 2 M or higher, there were few new genes detected in both libraries (Figure 2). The number of detected genes reached a plateau when 4 M tags were sequenced. No new genes were identified as the total tag number approached 4.7 M in CO library and 4.8 M in DS library. Therefore, the CO and DS libraries were sequenced to saturation, producing a full representation of the transcripts in the conditions tested.

Annotation analysis of the unique tags

Comparison of gene expression level between two libraries

Tag abundance appearing in library was used for estimating the expression level of gene mapped by tag. Differences of tag frequencies appearing in CO and DS libraries were used to determine the expression changes of genes in response to drought stress. The transcripts detected with at least two-fold differences in two libraries were shown in Figure 3 (FDR ≤ 0.001). The red dots (1,011) and green dots (505) respectively represent more and less abundant transcripts with more than two folds difference in DS library, designated as differentially expressed genes (DEGs, i.e. potential drought stress-responsive genes); while the blue dots represent transcripts with less than two-fold abundant difference between two libraries, which were designated as “no difference in expression”. In other words, a total of 1011 and 505 genes were probably up- and down- regulated under drought stress, respectively. The differentially expressed unique tags with more than five folds difference were shown in Figure 4. A total of 427 genes which were matched by about 0.36% total unique tags had a more than five-fold increase in expression abundance, and 123 genes matched by about 0.29% total unique tags had a decrease of abundance with more than 5 folds in the DS library, while the expression difference of 99.35% unique tags was within five-fold between two samples. Among 1516 potential drought stress-responsive genes, there were 157 genes up-regulated and 27 genes down-regulated with greater than hundred folds difference in DS library (Additional file 1) and 1258 genes showed significant similarity with known proteins in Nr database (Additional file 2).

Potential pathway influenced by drought stress

Identification of drought stress-responsive TFs

The ramie actin gene (CL5463.Contig2) with a similar value of transcripts per million clean tags (TPM) in DS and CO libraries was selected as the endogenous control of qRT-PCR. The t-test showed that the qRT-PCR Ct value of actin in CO and DS ramie had no difference (P > 0.05). Thus, the actin expression did not have differences between the DS and CO ramie. The qRT-PCR result was presented in Figure 5. Six and three TFs were up- and down- regulated with 2~6 folds under water deficit condition. The unigene9044 and unigene19721 were up-regulated with more than 40 and 80 folds under drought stress, respectively; whereas unigene19209 was down-regulated with greater than 80 folds. These result suggested that all these 12 TFs were authentically drought stress-responsive TFs.

Identification of 1516 potential drought stress-responsive genes by Illumina tag-sequencing

In China, almost 90% ramie distributes in Yangtze valley, which indicates that ramie has a poor eco-adaptability. The correlation between environment factors of ramie cultivation region and ramie fiber yield showed that the fiber yield severely depended on the rainfall of ramie growth region [13]. Ramie fiber extracted from stem bark is vegetative tissue and its yield is determined by the stem growth. In this study, severe inhibition of stem growth and significant decrease of fiber yield were observed in drought-stressed ramie. Except these morphological traits such as stem traits and fiber yield, a large number of physiological characters are easily influenced by water stress in ramie. Previous study showed that significant decreases in contents of chlorophyll a, carotenoid and endogenous GAs and relative water content, and increases in the activities of peroxidase, superoxide dismutase, and catalase and the contents of proline, malondialdehyde and soluble sugar were observed under drought stress [14]. Hence, sufficient water supply from growth environment is essential for ramie high production. However, in order to ensure the food security in China, the irrigable cultivated land was used to produce foodstuff and the ramie was mainly grown in un-irrigable dry land such as hill sloping land. Therefore, elucidating potential molecular responsive mechanism of ramie to water stress and improving its drought tolerance ability have important significance for ramie producing in China. In present study, on the basis of Illumina tag-sequencing technology, a total of 1516 potential drought stress-responsive genes were identified. The identification of these potential stress-responsive genes will be helpful for further understanding of ramie drought tolerance.

The Illumina sequencing technology is the next generation sequencing (NGS) which is a powerful tool utilized in many researching areas, including re-sequencing, micro-RNA expression profiling, DNA methylation, de novo transcriptome sequencing and whole sequencing [15–21], especially in the analysis of whole-genome expression profiling [22–25]. In current study, the NGS was used to identify the potential drought stress-responsive genes based on principle that tag frequencies of given gene can be used to estimate its transcript abundance. Theoretically, tags should be generated by NlaIII from the 3′-most ends of transcripts. Other tags may also be generated because of incomplete enzyme digestion in practice. Since only one tag could be generated in each transcript from any NlaIII site in a cDNA, one tag represented a transcript of a given gene. In other word, the total of all NlaIII tags’ copy number in a gene represented the transcript abundance of this gene in library. For alternative splicing, it was possible that there were some genes spliced as multiple transcripts. In our study, when one tag matched to multiple transcripts, this tag would not be used to estimate the expression level. Out of 43,990 reference genes [12], about 54.4% and 51.9% genes were matched by tags of DS and CO library. Several potential reasons were responsible for the residual genes which were not matched by tags. First, there were about 20% genes without restriction enzyme cutting site of NlaIII, which led to a fact that these 20% genes could not generate their tag. Second, the reference transcriptome used in this study were sequenced based on the RNA of several growth stages, including seedling, vigorous vegetative growth stage and fiber ripeness stage [12], while the RNA used to construct library in this study was only extracted from the tissue of 30-day-old ramie. Probably, some genes which only express in seedling and fiber ripeness stage will not appear in the library constructed in this study. Moreover, a large number of ESTs in reference transcriptome are partial sequence of genes. Hence, of these 43,990 ESTs in reference transcriptome, it is likely that several ESTs are sequenced from a common gene which only generates a type of tag by NlaIII from the 3′-most ends of gene. In other words, only one of several ESTs sequenced from a gene can be aligned with the tag of this gene and the others can not be annotated by tags, which is another major reason for the presence of about 40% genes un-annotated.

Transcription factors responding to drought stress in ramie

Transcriptome analyses using microarray technology, together with conventional approaches, have revealed that dozens of transcription factors (TFs) are involved in the plant response to drought stress [4, 9, 10]. Most of these TFs fall into several large TF families, such as AP2/ERF, bZIP, NAC, MYB, MYC, Cys2His2 zinc-finger and WRKY. The expression of TFs regulates the expression of downstream target genes which are involved in the drought stress response and tolerance. Up to now, hundreds of TFs were validated in its ability of drought tolerance by further study. Taking NAC TFs as an example, scores of NAC TFs in rice, Arabidopsis, wheat, maize and so on were found to respond to abiotic and biotic stress and over-expression of these TFs can improve the drought tolerance ability of transgenic plant [26–31]. In this study, a total of 24 TFs involving in several families such as NAC, MYB, HD-Zip, AP2/ERF and so on were found that they are probably differential expression (20 TFs up-regulated and 4 TFs down-regulated) between well-watered and drought-stressed ramie by Illumina tag-sequencing. Twelve TFs (3 NAC TFs, 2 HD-Zip TFs, 1 bHLH TF, 1 DOF TF, 1 C2H2L TF, 1 MYB TF, 1 AP2 TF, 1 ARF TF and 1 GRAS TF) were chosen for further expression analysis in well-watered and drought-stressed ramie by qRT-PCR and the result validated that all these 12 TFs were drought stress-responsive TFs.

Up-regulation of unigene19721 expression is probably responsible for the inhibition of the growth of drought-stressed ramie

Under drought condition, a major morphological characteristic of plants is dwarfism, which is considered as an adaptive change of plants to help them avoid high energy costs under unfavorable conditions [32]. Opposite to drought inhibiting the growth, GAs can stimulate stem elongation and promote the growth of plants [33]. There is a potential crosstalk between drought stress signal and GAs signal resulting in antagonic interaction to regulate the plant growth. DELLA protein is a negative regulator of GAs signal pathway and can inhibit the growth of plants. GAs signal can induce the destruction of DELLA protein and then relieve the repression function of DELLA [33]. Previous studies showed that DELLA protein can respond to abiotic and biotic stress [34], and the accumulation of DELLA protein markedly improves the ability of stress tolerance [35–37]. Therefore, DELLA protein not only functions as negative regulator to repress the growth of plants but also enhances the ability of stress tolerance of plants. In this study, a DELLA protein encoded gene, unigene19721, was found up-regulated expression with 335 folds under drought stress. The up-regulated expression of this gene in drought-stressed ramie was further confirmed by qRT-PCR. Probably, ramie increases its DELLA protein to enhance drought tolerance ability by up-regulating the expression of unigene19721 under water deficit, whereas the increase of DELLA protein leads to a corresponding inhibition of stem elongation and decrease of fiber yield in ramie.

Plant material, treatment of water stress and RNA extraction

Elite ramie variety Zhongzhu 1 was used in this study. Zhongzhu 1 is an elite variety with characteristics of high yield, good fiber quality and strong drought resistance and it has the largest growth area in China during recent years. The cuttage seedlings of Zhongzhu 1 were transplanted to pot in May 2011. In April 2012, the potted 30-day-old ramies were transferred to a movable rain-off shelter and were parted into two groups. One group (CO) where the ramie was grown under well water condition by daily watering was used as control, and the other group (DS) was treated with drought stress by controlling the relative water content of soil at a level of no more than 35%. Each group was planted with three replicates. After seven days the ramie suffering from drought stress, the CO and DS ramie tissues of three replicates including leaf, root, stem bast, stem xylem and stem shoot were individually collected. The sampled tissues were immediately frozen in liquid nitrogen and stored at −80° until use. The same tissue sample of three replicates of each group was mixed to extract RNA. Total RNAs of two treatment ramie were extracted from each tissue sample using TRIzol reagent (Transgene Company, Illkirch Graffenstaden Cedex, France) according to the manufacturer’s protocol. The RNA was equally pooled from the five tissues for cDNA library preparation.

Preparation of digital expression libraries

Sequence tag preparation was done with the Digital Gene Expression Tag Profiling Kit (Illumina Inc; San Diego, CA, USA) according to the manufacturer’s protocol. Six micrograms of total RNA for CO and DS ramie was individually purified using biotin-Oligo (dT) magnetic bead adsorption. First- and second-strand cDNA synthesis was performed after the RNA was bound to the beads. While on the beads, double strand cDNA was digested with NlaIII endonuclease to produce a bead-bound cDNA fragment containing sequence from the 3′-most CATG to the poly (A)-tail. These 3′ cDNA fragments were purified by using magnetic bead precipitation and the Illumina adapter 1 (GEX adapter 1) was added to new 5′ end. The junction of Illumina adapter 1 and CATG site was recognized by MmeI, which is a Type I endonuclease (with separated recognition sites and digestion sites). The enzyme cuts 17 bp downstream of the CATG site, producing 17 bp cDNA sequence tags with adapter 1. After removing 3′ fragments with magnetic bead precipitation, the Illumina adapter 2 (GEX adapter 2) was ligated to 3′ end of the cDNA tag. These cDNA fragments represented the tag library.

Illumina sequencing

Illumina sequencing using the HiSeq^TM 2000 platform was performed at Beijing Genomics Institute (BGI)-Shenzhen, Shenzhen, China (http://www.genomics.cn/) with the method of sequencing by synthesis. Briefly, PCR amplification with 15 cycles using Phusion polymerase (Finnzymes, Espoo, Finland) was performed with primers complementary to the adapter sequences to enrich the samples for the desired fragments. The resulting 105 base strips were purified by 6% TBE PAGE Gel electrophoresis. These strips were then digested, and the single chain molecules were fixed onto the Illumina Sequencing Chip (flow cell). Each molecule grew into a single-molecule cluster sequencing template through in situ amplification. Four color-labeled nucleotides were added, and sequencing was performed with the method of sequencing by synthesis. Image analysis and base calling were performed by using the Illumina Pipeline, and cDNA sequence tags were revealed after purity filtering. The tags passing initial quality tests were sorted and counted. Each tunnel generates millions of raw reads with sequencing length of 49 bp (target tags plus 3′adaptor). Each molecule in the library represented a single tag derived from a single transcript.

Gene annotation

“Clean Tags” were obtained by filtering off adaptor-only tags and low-quality tags (containing ambiguous bases). Comparison of the Clean Tags sequences with ramie transcriptome sequence [12] by BLASTN was carried out. All clean tags were annotated based on ramie reference genes. The number of annotated clean tags for each gene was calculated and then normalized to TPM (number of transcripts per million clean tags) [24, 25].

Identification of differentially expressed genes (DEGs)

FDR (False Discovery Rate) was applied to determine the threshold of P Value in multiple tests and analyses [39]. An “FDR < 0.001 and the absolute value of log2-Ratio ≥ 1” was used as the threshold to judge the significance of gene expression difference.

Real-time quantitative PCR (qRT-PCR) analysis

The CO and DS ramie were used for qRT-PCR analysis. Entire plants of six individuals (three CO plants and three DS plants) were individually sampled. The sampled tissues were immediately frozen in liquid nitrogen and used to extract RNA. For each sample, first-strand cDNAs were reverse-transcribed from RNAs treated with DNase I (Fermentas, Canada) by using M-MuLV Reverse Transcriptase (Fermentas, Canada) according to the manufacturer’s instructions. qRT-PCR was performed using an optical 96-well plate with an iQ5 multicolor real time PCR system (Bio-RAD, USA). Each reaction contained 1 μL of cDNA template, 10 nM gene-specific primers, 10 μL of SYBR Premix Ex Taq, and 0.4 μL of ROX Reference Dye (FINNZYMES, Finland) in a final volume of 20 μL. The ramie actin gene was selected as the endogenous control. The primer sequence of DEGs and actin gene were listed in Additional file 5. The thermal cycle used was as follows: 95°C for 15 min, followed by 40 cycles of 95°C for 10 s, 55°C for 20s and 72°C for 30 s. qRT-PCR was performed in triplicate for each sample. Relative expression levels were determined as described previously [40].

Pathway enrichment analysis of DEGs

where N is the number of all genes that with KEGG annotation, n is the number of DEGs in N, M is the number of all genes annotated to specific pathways, and m is number of DEGs in M. Q value was used for determining the threshold of P Value in multiple test and analysis [41]. Pathways with Q value < 0.05 are significantly enriched in DEGs.

Availability of supporting data

The raw data (tag sequences, counts and TPM value) has been deposited at the Gene Expression Omnibus [NCBI GEO] with accession number GSE46253, [http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE46253].