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

Abstract

Background

Ramie fiber extracted from stem bark is one of the most important natural fibers. Drought is a main environment stress which severely inhibits the stem growth of ramie and leads to a decrease of the fiber yield. The drought stress-regulatory mechanism of ramie is poorly understood.

Result

Using Illumina sequencing, approximately 4.8 and 4.7 million (M) 21-nt cDNA tags were respectively sequenced in the cDNA libraries derived from the drought-stressed ramie (DS) and the control ramie under well water condition (CO). The tags generated from the two libraries were aligned with ramie transcriptome to annotate their function and a total of 23,912 and 22,826 ramie genes were matched by these tags of DS and CO library, respectively. Comparison of gene expression level between CO and DS ramie based on the differences of tag frequencies appearing in the two libraries revealed that there were 1516 potential drought stress-responsive genes, in which 24 genes function as transcription factor (TF). Among these 24 TFs, the unigene19721 encoding the DELLA protein which is a key negative regulator in gibberellins (GAs) signal pathway was probably markedly up-regulated under water stress for a increase of tag abundance in DS library, which is possibly responsible for the inhibition of the growth of drought-stressed ramie. In order to validate the change of expression of these potential stress-responsive TFs under water deficit condition, the unigene19721 and another eleven potential stress-responsive TFs were chosen for further expression analysis in well-watered and drought-stressed ramie by real-time quantitative PCR (qRT-PCR) and the result showed that all 12 TFs were authentically involved in the response of drought stress.

Conclusion

In this study, twelve TFs involving in the response of drought stress were first found by Illumina tag-sequencing and qRT-PCR in ramie. The discovery of these drought stress-responsive TFs will be helpful for further understanding the drought stress-regulatory mechanism of ramie and improving the drought tolerance ability of ramie.

Background

Drought is one of the most common environmental stresses that affect the growth and development of plants [1]. The global scarcity of water resources has already become a severe environmental problem worldwide. Poor water management, increased competition for limited water resources, and the uncertain threats associated with global warming all highlight the looming water crisis that threatens agricultural productivity worldwide. It has become urgent to elucidate the responses and adaptation of crops to water stress, and improve the drought tolerance of crops.

Plant response to drought stress is a complex course, and several mechanisms known as drought escape, drought avoidance and drought tolerance are involved in adapting the environment of water deficit [2]. A great number of dynamic responses to water deficit at physiological, biochemical, and molecular levels are presented in plant, thus enabling them to survive under drought environmental conditions [3, 4]. Recently, expanding transcriptome data sets have uncovered a global picture of stress responsive genes in Arabidopsis [5], rice [6], maize [7], wheat [8] and other plants. These transcriptome data revealed that drought stress induced genes not only function to protect cells from drought stress through the production of important enzymes and metabolic proteins (functional proteins), but they also regulate signal transduction and gene expression in the stress response (regulatory proteins). The functional proteins include late embryogenesis abundant (LEA) proteins, a variety of transporters, enzymes involved in osmoprotectant synthesis, fatty acid metabolism, cellular metabolism, carbohydrate metabolism and secondary metabolism. Regulatory proteins that are activated in response to water stress, including transcription factors (TFs) such as DREBs (dehydration-responsive element-binding proteins), AREBs (ABA-responsive element-binding proteins) and NAC proteins, have been identified in plant [4, 9, 10]. Besides, many genes involved in growth and development, such as chloroplast, cell wall and plasma membrane proteins encoded gene, were down-regulated in response to drought stress [10].

Ramie (Boehmeria nivea), popularly called “China grass”, is one of the most important natural fiber crops. Ramie fibers, which are extracted from stem bast, are smooth, long and have excellent tensile strength. This high fiber quality is the major reason that ramie is widely cultivated in China, India, and other Southeast Asian and Pacific Rim countries. In China, ramie is the second most important fiber crop, with its growth acreage and fiber production being second only to those of cotton. Ramie has vigorous vegetative growth and can be harvested three times per year in China, and up to six times per year in well-watered cultivation environments in Philippines, which makes ramie produce a high yield of vegetative fiber. Therefore, enough water supplied by growing environment is essential to meet the requirement of vigorous metabolism for vegetative growth. When ramie suffered from water deficit, there were numbers of morphological and physiological changes in response to drought stress, such as leaf and root shape, malondialdehyde and proline contents, catalase activity and net photosynthetic rate in ramie [11]. However, up to now, none of genes involved in drought tolerance was identified and the potential drought stress-regulatory mechanism is still unknown in ramie. In this study, in order to identify the drought stress-responsive transcription regulator, the potential stress-responsive genes were identified on the basis of Illumina tag-sequencing at first; and then the differentially expressed TFs were screened and further validated by qRT-PCR. This study will be helpful for further elucidating the potential molecular responsive mechanism of ramie to drought stress and improving the drought tolerance ability of ramie.

Result

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.

Figure 1

The changes of ramie fiber yield and stem traits in response to drought stress. The error bar represented the standard error.

Tag identification and quantification

A total of 4,719,982 and 4,804,046 tags were sequenced in control ramie (CO) and drought-stressed ramie (DS) libraries, respectively (Table 1). After filtering out low quality tags (tags containing ‘N’ and adaptor sequences), 4,715,625 and 4,799,759 tags (noted herein as “clean” tags) remained in CO and DS libraries. To increase the robustness of the approach, single-copy tags in the two libraries (305,492 in CO and 319,431 in DS library) were excluded from further analysis. As a result, a total of 4,410,133 (93.44%) and 4,480,328 (93.26%) clean tags remained in the two libraries, from which 328,806 (CO) and 340,187 (DS) unique tags were obtained. Hence, only 6.56% and 6.74% tags in CO and DS libraries respectively were useless, which suggested that the sequence quality was excellent in the two libraries. There were 11,381 more unique tags in DS library than in the CO library, possibly representing genes related to drought tolerance.

Table 1 Illumina tags in the control (CO) and drought stress (DS) libraries

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.

Figure 2

Saturation evaluations for the CO and DS libraries.

Annotation analysis of the unique tags

The ramie transcriptome had been de novo assembled and 43,990 unique genes were identified and annotated for their function [12]. In order to annotate the function of the tags sequenced in DS and CO libraries, the unique tags were aligned with these 43,990 unique genes using BLASTn. Tags with a complete match or one base pair mismatch were further considered to be used for estimating the expression level of gene. The results showed that 43,085 (12.67%) unique tags were matched to 23,912 (54.36%) in DS library and 39,143 (11.90%) unique tags were matched to 22,826 (51.89%) expression genes in CO library (Table 2).

Table 2 Annotation of illumina 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).

Figure 3

Comparision of gene expression level between CO and DS libraries. Red dots represent transcripts more prevalent in the DS library, green dots show those present at a lower frequency in the drought stress ramie and blue dots indicate transcripts that did not change significantly. The parameters “FDR ≤ 0.001” and “log2 Ratio ≥ 1” were used as the threshold to judge the significance of gene expression difference.

Figure 4

Differentially expressed tags in DS tissue library. The “x” axis represents fold-change of differentially expressed unique tags in the DS library. The “y” axis represents the number of unique tags (log10). Differentially accumulating unique tags within 5-fold difference between libraries are shown in the red region (99.35%). The blue (0.36%) and green (0.29%) regions represent unique tags that are up- and down regulated for more than 5 fold in the DS library, respectively.

Potential pathway influenced by drought stress

The possible influence of drought stress on biological pathways was evaluated by enrichment analysis of DEGs. A total of 112 pathways were possibly affected by drought stress (Additional file 3). Pathways with Q value < 0.05 were significantly enriched by DEGs. Nine pathways may be severely influenced by drought stress for the significant enrichment of the DEGs (Q < 0.05) (Table 3). Among 9 enriched pathways, 5 pathways had more up-regulated DEGs; 3 pathways had more down-regulated DEGs; one pathway had a same number of up- and down-regulated potential stress-responsive genes. The Ribosome pathway enriched the most DEGs, followed by Starch and sucrose metabolism, Pentose and glucuronate interconversions, Phagosome, Other glycan degradation, Carbon fixation in photosynthetic organisms, Fructose and mannose metabolism, Ascorbate and aldarate metabolism and Riboflavin metabolism (Table 3).

Table 3 List of pathway significantly enriched in DEGs (Q < 0.05)

Identification of drought stress-responsive TFs

Twenty-four potential drought stress-responsive transcription regulators were identified by Illumina tag-sequencing (Additional file 4). Twenty transcription factors (TFs) showed more and 4 TFs showed less abundance in DS library. Among 20 TFs up-regulated potentially, the unigene19721 encoding DELLA protein which is a key negative regulator in gibberellins (GAs) signal pathway, had more abundance with 335 folds difference in DS library. It is possible that the up-regulation of unigene19721 expression is responsible for the inhibition of the growth of ramie and the decrease of fiber yield under drought stress. Therefore, the expression level of unigene19721 and another eleven potential stress-responsive TFs were further analysis by qRT-PCR (Table 4).

Table 4 The function annotated of TFs validated by qRT-PCR

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.

Figure 5

qRT-PCR analysis of twelve differentially expressed TFs. Data represent fold change of each DEG’s relative quantification in drought stress (DS) vs. control (CO) samples; the error bar represented the standard deviation.

Discussion

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.

Conclusion

In this study, a total of 1516 potential drought stress-responsive genes including 24 TFs were identified in ramie. Twelve TFs were further validated to involve in the response of drought stress by qRT-PCR. Among the 12 stress-responsive TFs, the unigene19721 encoding the DELLA protein which is a key negative regulator in gibberellins signal pathway was markedly up-regulated, which is probably responsible for the inhibition of the growth of ramie under drought stress. The identification of these candidate TFs which may contribute to drought tolerance in ramie will be helpful for further improving ramie drought tolerance ability.

Methods

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)

A rigorous algorithm to identify differentially expressed genes between two samples was developed [38]. P value was used to test differential transcript accumulation. In the formula below the total clean tag number of the CO library is noted as N1, and total clean tag number of DS library as N2; gene A holds x tags in CO and y tags in DS library. The probability of gene A expressed equally between two samples can be calculated with:

P\left(y|x\right)={\left(\frac{N2}{N1}\right)}^y\frac{\left(x+y\right)!}{x!y!{\left(1+\frac{N2}{N1}\right)}^{\left(x+y+1\right)}}

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

Pathway enrichment analysis based on KEGG (Kyoto Encyclopedia of Genes and Genomes pathway database http://www.genome.jp/kegg) was used to identify significantly enriched metabolic pathways or signal transduction pathways in differentially expressed genes comparing with the whole genome background. The calculating formula is:

p=1-{\displaystyle \sum _{i=0}^{m-i}\frac{\left(\begin{array}{c}\hfill M\hfill \\ \hfill i\hfill \end{array}\right)\left(\begin{array}{c}\hfill N-M\hfill \\ \hfill n-i\hfill \end{array}\right)}{\left(\begin{array}{c}\hfill N\hfill \\ \hfill n\hfill \end{array}\right)}},

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].