Genome-wide transcriptomic comparison of cotton (Gossypium herbaceum) leaf and root under drought stress

In this study, the 454 pyrosequencing platform was used for analyzing the comparative transcriptomic profiles of leaf and root tissues of 1-month-old cotton (Gossypium herbaceum) plants under drought stress. A total of 56,354 and 49,308 reads were generated from leaf and root tissues, respectively, and clustered into 6,313 and 5,858 unigenes. The differentially expressed unigenes that showed up-regulation (≥2-fold) or down-regulation (2≤-fold) were considered for further analysis. A total of 3,517 unigenes were differentially expressed in both tissues. The 1,528 genes specific to leaves and 1,128 specific to roots were obtained. The 28 biological pathways in two tissues were found to respond significantly to drought stress. A total of 289 in leaf and 277 in root unknown (novel) unigenes were found to be remarkably regulated by drought stress. Some key regulatory genes involved in abiotic stress such as WRKY, ERF, AP2 EREBP, MYB, and LEA were highly expressed in leaves. The genes RHD3, LBD, and transcription factor WRKY75, known for root development under various stress conditions, were expressed specifically in root. The genes related to chlorophyll a/b binding protein and photosystem-related proteins showed significant higher expression in roots and as compared to leaves. It can be concluded that cotton leaves are distinct from roots in terms of molecular mechanisms for responses to drought stress. Electronic supplementary material The online version of this article (doi:10.1007/s13205-014-0257-2) contains supplementary material, which is available to authorized users.


Introduction
Drought is one of the major constraints to cotton production worldwide. Cotton plants grown under unsuitable environmental conditions such as temperature extremes, water deficit, and salinity stress face reduced growth and productivity resulting from loss of cotton boll and altered fiber development. Choice of Gossypium herbaceum in these tracts was not preferential but obligatory due to harsh cotton growing conditions. In most regions of India, G. herbaceum is planted for its natural textile fibre. To cope with water stress, plants execute a number of physiological and metabolic responses (Arndt et al. 2001;Dias et al. 2003). The physiological and biological functions of plants are determined by expression pattern of genes which are influenced by several adverse environmental conditions (Shinozaki and Yamaguchi-Shinozaki 2007;Valliyodan and Nguyen 2006). Hence, an importance in stress-related studies is to detect the expression of unique genes or clusters of genes and to determine their expression pattern in specific plant tissue or organs which control the regulatory network. The plant leaves provide an adaptive mechanism for plants in water-deficit condition by reducing the leaf area, decreasing transpiration rate, and maintaining stomatal closure (Caldeira et al. 2014;Sanchez-Blanco et al. 2009). The small leaves or reducing leaf expansion is an important adaptive mechanism for drought tolerance (Caldeira et al. 2014). The proportion of leaf elongation and size has been attributed to root variability, root length and growth (Parent et al. 2010). The root architecture is determined by an endogenous genetic program as well as by external environmental factors (Ishida et al. 2008;Smith and Smet 2012). The few studies that have compared leaf and root tissue transcriptomes have highlighted the organ specificity of drought responses (Cohen et al. 2010;Libault et al. 2010;Milner et al. 2014;Narsai et al. 2010). Roots sense the edaphic water deficit, send chemical signals to the shoots, and maintenance of root growth despite reduced water availability can contribute to drought tolerance through water foraging (Lynch 1995). The comparative analysis of gene expression profiles in the roots of maize (Jansen et al. 2013;Yue et al. 2008), soybean (Libault et al. 2010;You et al. 2011), populous (Cohen et al. 2010) and Arabidopsis has revealed a considerable number of rootspecific genes and these are important for abiotic or biotic tolerance (Smith and Smet 2012). The various research groups using transcriptome studies have revealed several downstream components involved in the complex network of root formation (Ishida et al. 2008;Milner et al. 2014;Sengupta and Reddy 2011). A large number of genes encoding transcription factors (NLP, WRKY75, RAV and REM), osmoprotectants, ion transporters and heat shock proteins and pathways involved in hormone (ABA, ethylene and JA) biosynthesis and signal transduction are known to play an important role under drought stress (Narusaka et al. 2003;Ranjan et al. 2012b;Shinozaki and Yamaguchi-Shinozaki 2007;Tuteja and Sopory 2008). The genes involved in phenylpropanoid and flavonoid biosynthesis, pentose and glucuronate interconversions and starch and sucrose metabolism pathways are able to maintain adaptive condition during drought (Naika et al. 2013;Xiao et al. 2014). Cotton is an important source of natural fibre used in the textile industry and the productivity of the crop is adversely affected by drought stress (Trivedi et al. 2012;Zhu et al. 2013). Molecular processes in response to water-deficit stress have been studied at great length in cotton. Bowman et al. (2013) reported transcriptome analysis of upland cotton (G. hirsutum) root under drought stress. Chen et al. (2013) analyzed the drought-resistibility in two different cultivars of tetraploid upland cotton, using transcriptome profiling of leaf. Park et al. (2012) reported differentially expressed genes under water stress in root and leaf tissues in tetraploid cotton using cDNA-AFLP. Payton et al. (2011) investigated the transcriptional changes in the roots and leaves of diploid cotton under water-stress condition. The recent publication of the whole genome sequence of the cotton diploid relative Gossypium raimondii (Paterson et al. 2012) has expanded the use of NGS as a tool to study cotton development.
In our previous study we analyzed the drought-resistibility in different diploid cotton genotypes, and observed that the drought tolerance is due to several biochemical mechanisms working together (Ranjan et al. 2012a, b). In the present study, a genome-wide comparative analysis of genes in root and leaf tissue was performed to identify root and leaf specific genes and biological pathways that improve tolerance to drought stress. The comparative expression profiles of leaves and roots were useful in consolidating our knowledge of molecular mechanism of cotton pants in response to drought stress.

Drought stress treatment and plant growth
The genotypes of cotton (G. herbaceum), namely, RAHS-14 and GujCot-21 are drought sensitive and drought tolerant, respectively, were used for this study.
The seeds of both genotypes were sterilized and aseptically kept immersed in water for a day at 30°C and then placed for germination in a moist petri dish under the following conditions: 28/25°C as day and night temperature, 12 h of light and dark periods alternatively and relative humidity of 80 % in the first experiment. After 3 days, properly germinated seeds were transferred on paper rafts, which were submerged (3/4th) in Hoagland's media containing different percentages (2, 4, 6 and 8 %) of mannitol and mannitol free Hoagland's media was considered as control ( Fig. 1a, b). Then seedlings were allowed to grow until they had shown healthy growth in control condition then root lengths of grown seedlings were compared. In second experiment, plants were grown in the earthen pot, drought stress was given to the plants by withholding watering till soil moisture reaches below 30 % in pots and drooping effects on plant leaves became prominent (Fig. 1c, e), whereas the control pots were irrigated daily (Fig. 1d, f).

RNA isolation
The earthen pot grown plants were selected for the experiments after drought stress treatment. The first leaves of RAHS-14 and root of GujCot-21 plant were collected for RNA extraction. Total RNA of leaf and root tissues was extracted by Spectrum Plant Total RNA Kit Ò (Sigma-Aldrich) according to the manufacturer's instructions. After DNase I treatment, RNA was purified by QIAquick PCR purification kit (Qiagen). Quality and quantity of the purified RNA were determined by measuring the absorbance at 260/280 nm (A260/A280). RNA integrity was further verified by 1.5 % agarose gel electrophoresis.

Transcriptome sequencing
Total RNA (3 lg) from leaf and root tissue were reverse transcribed using a T7-Oligo (dT) Promoter Primer in the first-strand cDNA synthesis (Affymetrix). The double strand cDNA was synthesized and enriched using in vitro transcription (IVT) reaction (Affymetrix). cRNA (3 lg) was reverse transcribed in the first-strand cDNA synthesis step by using random hexamer primer, followed RNase H-mediated second strand cDNA synthesis and purified in column (QIAquick PCR purification kit) and further used for GS-FLX pyrosequencing. These sequence data were deposited in NCBI Sequence Read Archive (SRA, http:// www.ncbi.nlm.nih.gov) with accession number SRX038499 (leaves) and SRX038496 (roots).

Assembly of transcriptome using CAP3
The analyses of sequences were conducted using custom Perl scripts and publicly available software, R package (www.R.project.org). The sequences for leaves and roots were tagged with the library name and pooled exemplars were queried against the NCBI nucleotide database (NR) and in Arabidopsis TAIR-10 version using BLASTN at e-value 10 -10 and alignment length of more than 50 % of the query sequence for annotation (see additional information-1, Table 1). All the Gossypium ESTs available at NCBI were downloaded and pooled. The pooled exemplars were also queried against all public cotton EST database to identify new transcripts of Gossypium.

Differential expression analyses
The reads of both the tissues were tagged and pooled to create a dataset of 105,662 reads and used for digital profiling. These reads were assembled using the CAP3 programme at an overlap of 100 bp and 80 % identity. These reads assembled into 12,171 contigs which include only those contigs that have more than five reads. We calculated TPM value and R value using the R statistics for contigs and those with R value C3 and fold change C2 were considered significantly differentially expressed contigs. These filtered contigs were annotated using BLASTN against NCBI nucleotide (NT) database and BLASTX against NCBI nonredundant proteins (NR) (see additional information-2).

Heat map analysis
For visualization of the significant comparisons, heat maps of the significant genes after FDR adjustment were produced with the heat maps two package in R. Hierarchical clustering of individual sample with 1,000 bootstrap replications was performed with the R package pvclust and heat map were sorted accordingly. To visualize clusters of genes expression, we grouped the z-transformed expression ratios by using k-means in R.

GO analyses
The GO annotations for the sequences were derived using their UniProt annotation. The UniProt database was used as it has extensive GO mapping. The GO annotation for level 5 was extracted for each library and used for further analysis. Pathway analysis was performed to find out the significant pathways of the differentially expressed genes in roots and leaves according to the KEGG databases. The two-way Fisher's exact test and v 2 -test were used to classify the significant pathways, and the threshold of significance was defined by p \ 0.05 and FDR \0.05. The FDR value was used to correct the 'p value' (Garcia-Alcalde et al. 2011).

MapMan analysis
To adapt MapMan software to cotton (G. herbaceum) the major MapMan BINs and their sub-BINS used for the classification of Arabidopsis genes were transferred to the G. herbaceum system. The classification based on the similarity to Arabidopsis proteins was then automatically checked against domains that have been assigned to a MapMan category (Usadel et al. 2006). Furthermore, cotton tentative consensus sequences that did not meet the prerequisites for an automatic draft assignment were classified manually according to their respective TIGR annotation (TIGR release 10). Finally, all draft assignments were corrected manually for potential mistakes.

Comparison of drought tolerance of GujCot-21 and RAHS-14
The GujCot-21 and RAHS-14 exhibited contrasting difference in their root growth under osmotic stress and control condition. GujCot-21 has longer root length as compared to RAHS-14 at 6 % mannitol concentration (Fig. 1a, b). In earthen pot experiments, RAHS-14 showed prominent effect of drought stress (Fig. 1e) as compared to control (Fig. 1d). However, GujCot-21 showed much better development, less wilting and higher biomass as compared to RAHS-14 (Fig. 1c). In addition, GujCot-21 exhibited better drought tolerance than RAHS-14 as observed by visual comparison of leaf senescence and wilting symptom in two genotypes (Fig. 1c, f.) Overall features of the drought stress-responsive expression profiles in leaves and roots In leaf tissue 1,528 genes were identified; out of these 289 genes showed no hits to any protein or nucleotide in the database (NR-non redundant and NT-nucleotide data base of NCBI) and the 26 genes were hypothetical or unknown proteins were obtained (see additional information 2). The majority of the annotated genes were from the photosynthesis pathway, with high expression in ribulose 1-5 bisphosphate carboxylase/oxygenase activase, photosystem II D and chlorophyll a/b binding proteins. The senescenceassociated proteins were also observed in the differentially expressed genes. In root library 1,128 genes were identified and out of these 275 genes did not have any hit to any of the database (NR, UniProt, NT) and about 8 % genes matched with predicted/hypothetical proteins. The other genes related to drought stress that were upregulated in the root were ascorbate peroxidase, cysteine protease, deltatonoplast intrinsic proteins, Lea proteins, etc. (see additional information 2). Roots behave in a more focused response under drought stress compared to leaves Among the cotton genes, 6,313 and 5,858 from leaves and roots, respectively, were either twofold up-or downregulated with adjusted p B 0.05. Together, 12,171 genes were either significantly up-or down-regulated under drought stress (Table 1). Thus 3,588 genes were differentially regulated in both leaves (1,528) and roots (1,128). However, various tissues behaved differently during the stress response. In terms of the number of genes, there appeared to be a balance between leaves and roots, while in terms of induction high number of genes was primarily observed in the root (see additional information 2). To have a general view of the whole picture, we grouped and classified them into Gene Ontology terms then biological processes were further examined. There were five significantly represented GOs in leaves while roots showed eight indicating the earlier response of roots. It is likely because of the faster dehydration on root represented by shorter root length (Table 2; Fig. 1e); in other words, roots are more water stressed than leaves. In the field, root also is the earlier organ under drought because plants uptake water through their roots. When soils become dry, roots sense this water deficit and produce chemical signals to transport to leaves via the xylems resulting in physiological changes in leaves (Arndt et al. 2001;Sengupta and Reddy 2011;Tuteja and Sopory 2008). Further we analyzed these genes by MapMan, This finding suggested that, the leaves differentially express genes in a more random/scattered manner that lead to less statically significant results. Compared to that, roots seemed to behave in a more focused response. Unlike the responses of the leaves, root tissue showed notable expression of stress related and protein signaling pathway (Fig. 2). The results suggest that, under drought stress the leaves try to maintain normal functional condition and re-direct to the whole system on stress response, whereas, root tissues became more active. The heat map analysis showed that the photosynthesis and related processes were among enhanced processes in the roots, while it was repressed in leaves, In particular, the simultaneous reduction in leaves and induction in root of photosynthesis as well as secondary metabolism suggests the possibility of biological processes functioning in an opposite and tissue-specific manner (Fig. 3).

Differentially regulated cellular metabolism in leaves and roots
To identify the stress affected biochemical network of cells, we retrieved Arabidopsis homologs of cotton genes from PLEXdb (http://www.plexdb.org/) and calculated the effected pathway using KOBAS. Table 3 shows the 28 biological pathways were commonly affected in roots and leaves and the majority of affected pathways were observed in leaves. Some pathways were significantly affected in both leaves and roots; they're involved in energy production, biosynthesis of phytohormones, amino acids, pigments and many others. Not only did each tissue differently regulate genes within the shared pathway, they also possessed tissue specific reactions. Specifically, leaves showed strong induction of choline and trehalose biosyntheses when cellulose and initial step of fatty acid biosynthesis were negatively regulated only in root (Table 3).

Unique and shared transcription factors in roots and leaves in response to drought stress
TFs are the main components to understand the complexity of expression of stress induced genes in their signaling network as suggested in various studies (Agarwal et al. 2006a, b;Ranjan et al. 2012b). In order to further investigate on these molecules, we determined the number of unique and shared TF's expressed in leaves and roots using conserved domain database in NCBI. Among 29 groups that are affected by drought, three and six were uniquely expressed in leaves and roots, respectively. My-bDNA-binding bHLH, HSF and WRKY were the most strongly affected groups even though the expressions were different in both the tissue. Transcription factors, such as WRKY 75, RAV, REM, CAMTA, SBP, NLP and G2 like, Orphan, TCP were uniquely expressed in root and leaf tissue, respectively (Table 4). In general the photosynthesis related genes are expressed at a lower level in roots as compared to the expression in leaves. However, in this study the expression pattern is very different. Photosynthesis genes that are encoded separately in plastid and nuclear genome are responsible for photosynthesis core subunits, and light harvesting complex proteins and pigment biosynthesis; respectively (Gutierrez-Nava et al. 2004;Kino-oka et al. 2001). The genes related to chlorophyll a/b binding protein and photosystem-related proteins showed significant higher expression in roots and as compared to leaves (Fig. 4). It is well known that light and plastid development are both crucial for the expression of nuclear photosynthesis genes. Although in our analysis no biological process or biological pathways shows plastid related biosynthesis in root, we found that biological pathway related to ''ribulose-1,5-bisphosphate carboxylase'' (RuBisCO) was significantly high in leaves.

Defense genes dominantly expressed in roots as compared to leaves
Through long-term evolution and adaptation to extreme conditions, cotton has been found to be rich in resistance genes for abiotic stresses (Wendel and Cronn 2003). When cotton plants face water stress, due to the change in external environment, plants are able to activate a series of mechanisms to respond drought stress (Park et al. 2012;Ranjan et al. 2012a, b). Cells recognize external stresses and induce many signal proteins, including signaling receptor kinases, calcium-dependent protein kinases and G-proteins. In our results showed that the transcripts of    these signal proteins accumulated in both tissues during drought stress. Genes for peroxidase, glutathione peroxidase, redox metabolism and superoxide dismutase, were predominantly expressed in root tissue (Table 5). These genes are involved in scavenging reactive oxygen species (ROS), and play a pivotal role as triggers of gene expression during abiotic stresses (Miller et al. 2010;Ranjan et al. 2012b).

Discussion
The advent of pyrosequencing has enabled to obtain the transcriptional changes of entire genome of a particular plant or specific plant organs in response to various stresses. However, the whole genome sequences are still not available for G. herbaceum species of cotton, leading to a dependence on mRNA sequences by pyrosequencing or collections of ESTs assembled from random cDNA libraries. In this study, we used the 454 pyrosequencing platform for transcriptome sequencing of cotton (G. herbaceum) leaf and root tissues. It was observed that GujCot-21 and RAHS-14 showed substantial differences in root growth under osmotic stress (Fig. 1a, b) and leaf senescence under drought stress condition (Fig. 1c, f). The leaf and root both tissues shared some common gene expression during stress, with the purpose of enhancing protective systems. However, these two tissues appeared to respond differently to drought conditions. Leaves induced more genes than root but those genes were scattered in many processes, most significantly were productions of osmoprotectant and senescence associated genes. The genes expressed in roots mainly associated with antioxidant and photosynthesis related pathways (Table 5; Fig. 4). Chen et al. (2013) recently reported down regulated expression of chloroplast-/plastid-related genes in leaf tissue of drought tolerant cultivar of cotton as compare to drought sensitive cultivar, while Bowman et al. (2013), observed higher expression of genes in root associated with starch and sugar metabolism is similar to our results. Besides, the two tissues can affect each other via the signaling and transportation system. Through long-term selection and adaptation to extreme conditions, cotton has been found to be rich in resistance genes for a range of abiotic stresses, including drought, high temperature and salinity (Ranjan et al. 2012a;Zhou et al. 2014). Chen et al. (2013) results also suggested that high number of genes are expressed temporally and specifically in plant response to drought treatment. Transcription regulation plays a central role in stress signal transduction pathways. In this study, we found 29 group of transcription factor genes differentially regulated at different levels in roots and leaves in response to drought stress. Among these, TF genes some were unique to leaf and root tissues. This finding reveals that a large group of TF genes are involved in the transcription regulation in response to drought stress in leaves and roots in a specific manner when cotton plant face drought stress. Our study also showed that a large proportion (289 in leaves and 278 in roots) of expressed unigenes have no functional annotation or have unknown function. It is also noteworthy that many of these unigenes showed remarkable expression changes in response to drought stress, especially with some having a unique expression in leaf and root tissues. Our results confirm differential expression of transcription factors in response to water-deficit stress, as has been observed in many plant species, including cotton (Park et al. 2012, Bowman et al. 2013). Therefore, it is possible to identify some drought responsive genes unique to cotton. Further analysis of the functions and expression controlling mechanism of these genes in cotton would not only supply the opportunity of isolation and identification of novel genes, but also enhance our further understanding of specific mechanism of drought tolerance in cotton.

Conclusions
To date, a limited number of studies on drought stress mediated gene expression in cotton have been reported. In this study we described an analysis of gene expression in cotton leaves and roots in response to drought stress and intended to carry out a comparative transcript profiling between the leaves and roots. We focused on a set of transcripts that exhibited unique expression in leaves and roots of cotton under drought stress. This is the first trial to comprehensively explore the genome-wide gene expression patterns in leaf and root of drought responsiveness in cotton. Our results show that cotton leaves are distinct from roots in terms of molecular mechanisms for responses to drought stress.