Transcriptional profiling of sugarcane leaves and roots under progressive osmotic stress reveals a regulated coordination of gene expression in a spatiotemporal manner

Sugarcane is one of the most important crops worldwide and is a key plant for the global production of sucrose. Sugarcane cultivation is severely affected by drought stress and it is considered as the major limiting factor for their productivity. In recent years, this plant has been subjected to intensive research focused on improving its resilience against water scarcity; particularly the molecular mechanisms in response to drought stress have become an underlying issue for its improvement. To better understand water stress and the molecular mechanisms we performed a de novo transcriptomic assembly of sugarcane (var. Mex 69–290). A total of 16 libraries were sequenced in a 2x100 bp configuration on a HiSeq-Illumina platform. A total of 536 and 750 genes were differentially up-regulated along with the stress treatments for leave and root tissues respectively, while 1093 and 531 genes were differentially down-regulated in leaves and roots respectively. Gene Ontology functional analysis showed that genes related to response of water deprivation, heat, abscisic acid, and flavonoid biosynthesis were enriched during stress treatment in our study. The reliability of the observed expression patterns was confirmed by RT-qPCR. Additionally, several physiological parameters of sugarcane were significantly affected due to stress imposition. The results of this study may help identify useful target genes and provide tissue-specific data set of genes that are differentially expressed in response to osmotic stress, as well as a complete analysis of the main groups is significantly enriched under this condition. This study provides a useful benchmark for improving drought tolerance in sugarcane and other economically important grass species.


Introduction
Saccharum officinarum L. (Sugarcane) belongs to the Poaceae family, an economically important group of tropical and subtropical grasses, comprising plants such as maize, wheat, rice, and sorghum.Sugarcane is a key plant for the global production of sucrose; as well to the bioethanol industry, animal feed and molasses [1,2].Moreover, sugarcane is one of the most important productive crops around the world, yielding an average of 1.66 billion of tons per year [3].
Sugarcane cultivation is severely affected by drought stress, and this is the major limiting factor on its productivity worldwide [4].The global climate change projections in the last years points to a bleak scenario where drought is the major problem for sugarcane crop production [5,6].Drought has affected large agricultural areas and rural communities, resulting in imbalance at micro and macroeconomics levels [7].Even more, drought projections point to a 25% decrease in production by 2080 [8,9].
In recent years, sugarcane has been subjected to intensive research focused on improving its resilience against water scarcity; particularly the molecular mechanisms in response to drought stress have become an underlying issue for its improvement.Relatively few droughtresponsive genes have been described so far, and these studies have focused only on leaves [10][11][12][13].Thereby, novel methodologies and techniques have been useful for the identification of genes that respond to drought stress and they have been used to improve stress tolerance and raise crop productivity [14].RNA-Seq can analyse species without a sequenced genome on a genomic scale and can provide useful transcriptional data for a deeper understanding of molecular processes governing the plant response to a specific stress [15,16].Nowadays, transcriptomic approaches on sugarcane by using RNA-seq have been centred on sugar accumulation in the cane and leaves development [17,18].
In the present study, we characterized the global transcriptional changes in leaves and roots of sugarcane var.Mex 69-290 under osmotic stress.De novo transcriptome analysis was conducted by using plantlets subjected to osmotic stress conditions induced by PEG 8000 treatment.The sugarcane transcriptomes were analysed on leaves and roots after 24, 48, and 72 h of stress induction, and using a non-stressed group as a treatment control.This work represents the most detailed transcriptomic description of sugarcane under water stress condition, up-to-date.We are providing a comprehensive catalog of up and down-regulated genes for leaves and roots, which could be a useful resource for future genetic plant breeding programs in sugarcane as well as other economically important grass species.

Physiological characterization of sugarcane plantlets under osmotic stress
Morphological differences were observed between T0 and T1, differences in treatments T2 and T3 included: curled leaves with a generalized loss of turgor and became chlorotic with some degree of apical dead tissue.Additionally, in T3 the root system showed some degree of apparent necrosis (Fig 1A).
These morphological observations are correlated with the physiological measured parameters.For T3 treatment a statically difference in RWC was found (Tuckey's HSD, p<0.05) in comparison to T0, T1, and T2, with a reduction of 31.17% with respect to T0, indicating a loss of the cell turgor (Fig 1B).The osmotic potential showed a reduction of approximately 60% in T3 in relation to the control (Fig 1C).Statistical differences in the CO 2 assimilation were found in T2 and T3, with reductions of 62.7% and 76%, respectively (Fig 1D).Transpiration rate was significantly different in T2 and T3, with reductions of 30% and 39%, respectively (Fig 1E).Proline content was tested and showed increased in both leaves and roots during treatment (Fig 1F).Leaves showed the biggest increment in proline during T3 treatment, increasing from 3.5 to 6 nmol mg -1 FW in T3.Higher accumulation of proline was detected in roots with 3 and 6 fold increments for T2 and T3, respectively.

Transcriptome analysis
The sugarcane transcriptome consisted in 252,702 transcripts (all assembled sequences containing isoforms) and 140,339 unigenes (the longer fragments without N bases).The summary of our transcriptome is shown in Table 1.After de novo assembly each unigene was compared by BLASTX (E-value 1e -20 ) against several databases.A total of 31,130 unigenes (22.18%) were annotated against the PlantRef protein database from NCBI.
All unigenes were mapped against the GO terms database (S1 File).A total of 19,372 unigenes were classified in 29 subcategories belonging to three major GO categories: Biological Process, Molecular Function, and Cellular Component (Fig 2A).Among the Biological Process category, "Organic substance metabolic process (16.1%)", "Cellular metabolic process  (15.8%)", and "Primary metabolic process (15%)" were the main functional groups, followed by "Single-organism cellular process (10.8%)" and "Nitrogen compound metabolic process (10.6%)".For the Cellular Component, "Intracellular (37.9%)" and "Membrane-bounded organelle (13.2%)" were the most abundant subcategory in this group, "Organic cycling compound binding (20%)" and "Heterocyclic compound binding (20%)" were the two main groups into the Molecular Function category.The length of unigenes ranged from 201 bp to 14,847 bp, with an average of 1,084 bp and a N50 length of 1,213 bp (N50 means that 50% of the entire assembly is contained in Transcripts equal to or larger than this value; Fig 2B).Among the annotated unigenes, BLAST result had homology to Zea mays, Sorghum bicolor, and Setaria italica grass species (Fig 2C ), higher identity of the matches was against S. bicolor (Fig 2D).All unigenes were annotated by mapping them to several public databases showed in Table 1, and results were depicted by using Venn diagrams (Fig 2E) and also, different Venn diagrams were plotted for the easiest comparison among databases (Fig 2F and 2G).

Sample clustering analysis for sugarcane libraries
We used principal component analysis (PCA), gene expression patterns during osmotic stress, and hierarchical clustering analysis (HCA) to obtain time-and tissue-specific gene profiling.PCA shows similarity among the non-stressed samples for both leaf and root tissues and was clustered as a neighbouring group, showing a similar transcript level in comparison to stressed samples.
The stressed leaf samples (24, 48, and 72 h after stress induction (ASI)) were grouped together as a negative correlated cluster, in comparison to the root stressed samples, that were grouped apart from leaves, as a tissue specific cluster (Fig 4A ).
From here on we will refer to the following codes: L = leaf, R = root, L0 and R0 = non-stress treatment, L1 and R1 = 24 h ASI, L2 and R2 = 48 h treatment ASI, L3 and R3 = 72 h treatment ASI.Similar to PCA, the HCA shows that L0 was assigned to the adjacent node of R0, which corresponds to non-stress treatments.L1, L2, and L3 were grouped into the same cluster, corresponding to stress treatment in leaf tissue; while R1, R2, and R3 were grouped in a cluster corresponding to stress treatment in root tissue.This finding demonstrates the gene expression specificity for each tissue during osmotic stress.During the stress treatments, similarities in gene expression pattern were observed between L1, L2, and, L3, and R1, R2, and R3.On the other hand, L0 and R0 were clustered together showing high r values above the mean, suggesting a similar pattern in gene expression for a non-stress condition for leaf and root (Fig 4B).The transcriptome data shows that gene expression is significantly altered from the first time of stress imposition (L1 and R1) for both leaf and root tissues of sugarcane sprouts.

Global analysis of Differentially Expressed Genes (DEGs) during osmotic stress
Identification of sugarcane DEGs associated with osmotic stress was done by analysis of gene expression levels and it was quantified by RSEM software v1.2.27 [20], all read counts were normalized to FPKM values.Pearson's correlation coefficients were higher than 0.  partitioned into six expression clusters (log2-transformed) with similar expression profiles (Fig 5B -5G).Each profile represents a group of genes that exhibited similar expression trend.The six clusters were generated by cutting the hierarchically clustered gene tree (5A) at 60% of the tree height (recommend; see methods).In cluster 5B, where most of the genes were grouped, the mean suggests few significant changes among treatments.In cluster 5C, the 4,630 clustered showed a down-regulation in the stress treatments for both leaf and root tissues.In contrast, cluster 5D and cluster 5E showed an opposite trend, where genes in cluster 5D were down-regulated on stress treatments and genes in cluster 5E were up-regulated on the stress treatments.The most drastical regulation trend was shown in clusters 5F and 5G, where genes were quickly up-regulated by stress treatments.Only genes belonging to root tissues in cluster 5G showed a low increment in the transcriptional responses.

DEGs in response to osmotic stress treatments in sugarcane
DEGs of leaves and roots grouped out of all the stress-time series.The volcano plots show that more genes were up-and down-regulated in T2 (48h vs. control) than in T1 (24 h vs. control), or T3 (72 h vs. control) for leaf tissues (Fig 6A).Red dots inside volcano plots show DEG with a Fold Change > 2 and a FDR<0.001.The volcano plots show that large number of genes are differentially expressed across the stress-treatments.According to this, we have that T2>T1 >T3 for down-regulated genes and T2>T3>T1 for up-regulated genes in leaf tissues.For root tissues a lower amount of DEGs was observed compared to leaves, the relationship for downregulated genes was T1>T3>T2 and T3>T1>T2 for up-regulated genes in root tissues.A total of 11,796 genes were differentially regulated in leaves, of which 7,441 were down-regulated, and 4,355 were up-regulated.For root tissues, a total of 8,611 genes were differentially regulated, of which 3,970 were down-regulated, and 4,641 were up-regulated.(Fig 6B ).A total of 95 down-regulated genes were shared among leaf and root tissues for all three stress treatments and 231 up-regulated genes were shared for both leaf and root tissues, in all the timestress treatments.
For leaf tissues, 3,539 transcripts were found to be differentially expressed in T1 treatment (24 h ASI), which included 1,220 up-regulated genes ( In the case of the up-regulated genes in leaf tissues, 887 genes were shared between T1 and T2, 802 genes were shared between T2 and T3, and 617 genes were shared between T1 and T3.In addition, 536 genes were shared between the three treatments (Fig 7A ) and we decided to call it the "CORE of Up-regulated Leaf Genes" (CULG).For the down-regulated genes in leaf tissues, 1,765 genes were shared between T1 and T2, 1,393 genes were shared between T2 and T3, and 1,223 genes were shared between T1 and T3.In the same way, 1,093 genes were shared Among the up-regulated genes in root tissues, 982 genes were shared between T1 and T2, 912 genes were shared between T2 and T3, and 878 genes were shared between T1 and T3 (Fig 7C ).Also, 750 genes were shared between the three treatments (Fig 7C) defined as "CORE of Up-regulated Root genes" (CURG).For the down-regulated genes in root tissues, 683 genes were shared between T1 and T2, 742 genes were shared between T2 and T3 treatments and 713 genes were shared between T1 and T3.In the same way, we determined 531 genes shared between the three treatments (Fig 7F) and was defined as the "CORE of Down-regulated Root Genes" (CDRG).
There are 231 genes shared between the COREs of up regulated genes for root and leaf tissues, coded here as "CORE of Osmotic Stress Up-regulated Genes" (COSUG; Fig 7B).On the other hand, we found 95 genes shared between the COREs of down regulated groups for both leaf and root tissues, and was defined as the "CORE of Osmotic Stress Down-regulated Genes" (COSDG; Fig 7E).
There were 1,956 assembled unigenes corresponding to transcription factors (TF), of which 148 were regulated in T1, 174 TF were regulated in T2, and 127 TF were regulated in T3 treatment for leaf tissues (Table 2).In root tissues, 80, 70, and 84 TF were regulated in T1, T2, and T3, respectively.For the up-regulated TFs in leaves, 15 of them were shared among treatments (CULG) and 62 down-regulated TFs were shared among treatments (CDLG).In roots, 19 up-regulated TFs were shared among treatments (CURG) and only 15 down-regulated TFs were shared in the three treatments (CDRG).Interestingly, only two TFs were shared among all the up-regulated genes in leaf and root tissues (COSUG), the gene NAC25-lke and an HSF-C2a.For the down-regulated genes, only two TFs were shared among all the stress treatments for tissues (COSDG), an ERF27-like and a DREB1G-like.

Genes involved in abiotic stress responses
We analysed individual datasets for the up-and down-regulated DEGs per each stress treatment of DEGs (CULG, CURG, CDLG, CDRG, COSUG, and COSDG; see Fig 7B and 7E).Several genes that play an essential role in plant maintenance and development were transcriptionally affected due to osmotic stress conditions.In CULG and CDLG sets (cores of Up-and down-regulated genes in leaf tissues respectively), we found several TFs families associated to stress responses like bZIP, ZF, HSF, NAC, MYB, PHR1, and AP2/ERF.All the TFs were regulated during the T1, T2, and T3 treatments (Table 2).
Only a few DEGs were found to be associated with secondary metabolism.In the current study, isoflavone 2-hydroxylase and chalcone synthase 1 markedly increased with osmotic stress imposition and reached the highest levels in T3 treatment.Several other genes were found to be importantly related to stress regulation (see Table 2).
Among the CURG and CDRG datasets (cores of Up-and down-regulated genes in root tissues respectively), we detected stress-related TFs representatives as NAC, HSF, MYB, bHLH, and AP2/ERF (Table 3).All of these TFs were significantly regulated during T1, T2, and T3 treatments.Others several stress related genes were regulated in root tissues after osmotic stress imposition like dehydrin DHN4 and dehydrin COR410, a cationic peroxidase SPC4, desiccation-related PCC13-62 protein, LEA 3 protein, BAX-Inhibitor protein, anthocyanidin reductase, among others.
Selected genes (NAC-25 TF, bidirectional sugar transporter SWEET6a, PIP2-1, and Abscisic Acid 8-hidroxilase 3) were selected for RT-qPCR analysis to confirm the reproducibility of gene profiles detected by bioinformatic analysis (S2 Fig) .NAC-25 TF, bidirectional sugar transporter SWEET6a, PIP2-1, and Abscisic Acid 8-hidroxilase 3 (this gene was used for qRT-PCR analysis in leaf and root tissues).Although the exact expression pattern of selected genes varied among the RT-qPCR and RNA-Seq analysis, their expression was quite similar in both cases.

GO functional analysis for DEGs
All the DEGs were annotated against the GO database (http://www.geneontology.org/) by Blast2GO software [21], to reveal significantly enriched GO terms in DEGs, besides, we also divided these terms into up-regulated and down-regulated groups.The GO enrichment analysis was conducted for each up and a down-regulated set of genes for leaf and root tissues using the Fisher's Exact Test.The results were subsequently adjusted using a FDR< 0.05 correction.We further analysed the functional enrichment for the COSUG and COSDG groups (Fig 8).
As shown in Fig 8, response to heat, water deprivation, abscisic acid, secondary metabolite biosynthetic process (as flavonoids), regulation of transcription, glutathione metabolic processes, among others, were the main enriched GO terms for up-regulated DEGs.On the other hand, response to light stimulus, anion transport, sucrose and starch metabolic processes, cellulose biosynthetic processes, among others, were the primary enriched GO terms for downregulated core of DEGs.

Physiological responses of sugarcane under osmotic stress treatments
Tissue culture is a great tool for the assessment of many hypotheses because a homogenous population of plants can be easily controlled, as well as the conditions of the assay [22,23].In this context, PEG has been used for several groups to mimic drought conditions in vitro and ex vitro [22,[24][25][26][27][28].PEG is a non-penetrating, inert, and non-phytotoxic agent that reduces the osmotic potential of culture medium inducing an osmotic stress in plant cells.Nevertheless, an outdoor evaluation of the in vitro tested condition is needed for a better understanding of the real and more complex processes that are triggered by drought conditions [29].
In this study, we describe the response of 90 days-old sugarcane plantlets after severe osmotic stress.The decrease in RWC is dependent to the time of stress imposition and reached its lowest value in T3 treatments with 31.17% less than T0.Comparatively, [30] reported a decrease of 25% in the RWC during mild stress (3 DWW; days by withholding water) and 33% during severe stress (9 DWW).
Sugarcane plantlets were significantly affected in CO 2 assimilation and transpiration rate after 48 h of osmotic stress.This observation can be explained by taking into account that drought conditions prevent CO 2 from entering the leaves, influence absorption of CO 2 by carboxylation centre result in a decrease of net photosynthetic rate [31].Over production of proline results in augmentation level tolerance to osmotic stress in plants.Here, proline content increased in both leaves and roots along the stress.Similar physiological parameters were reported in [30] for sugarcane variety GT21, showing significant alterations in response to drought stress.Our physiological data suggest a significant change in the regulation of molecular mechanisms that trigger the physiological responses caused by sensing the osmotic stress imposed by PEG 8000.

The sugarcane transcriptome
The osmotic stress is a principal component of many abiotic stresses, which interferes with the capability of plants to sense cell wall damage due to reduced turgor originated in an inappropriate activation of defence mechanisms [32].Recently, important efforts have been made towards the identification of stress-regulated genes aiming to improve the resilience capability against drought in economically important crops [33][34][35].
Rocha et al. reported 93 differentially expressed genes in sugarcane plants exposed to drought using macro-array technology which contained 1,545 probes [10]; among the DEG were hsp70, expansines, tyrosine phosphatases, catalases, Fructose-bisphosphate aldolase, trehalose metabolism-related proteins, most of which were found in our analysis.Similar to [10], Rodrigues and collaborators found 165 DEG associated to water deficit by using macroarray membranes containing 3,575 probes [11]; among these DEG were found DNAj proteins, PP2C, AP2 transcription factors, calmodulins, peroxidases, hsp70, Fructose-bisphosphate aldolase, among many others, similar to our results.Li and collaborators reported that 1,501 genes were differentially regulated (9.6% of the tested genes) in leaves under water-deficit by using a microarray technique [30], many of them involved in metabolic pathways such as biosynthesis of secondary metabolites, ribosomes, carbon metabolism, among others.
In the present study, a total of 11,796 and 8,611 genes were identified as DEG in leaf and root tissues, respectively.The PCA and the hierarchical clustering showed a tissue-specific expression pattern of the DEGs.The PCA and HCA revealed a significant change in gene expression after the first 24 h of stress imposition for both tissues, and these changes were specific according to the tissues.These results are visually evident on our heatmap (Fig 5A).The non-stressed tissue samples had very similar expression patterns in comparison to the stressed samples.The clusters E, F, and G (with 452, 161, and 31 genes; Fig 5B ) showed a rapid change in the expression patterns starting at T1 treatment, and with little changes for T2 and T3 stressed samples.

DEGs and pathways regulated by osmotic stress conditions
The GO enrichment analysis shows that several of the up-regulated sequences were grouped into many stress-related pathways as response to abscisic acid, water deprivation, heat, secondary metabolite biosynthesis, oxidation-reduction process, cellular chemical homeostasis, regulation of transcription, among others.On the other hand, several down-regulated sequences were grouped into photosynthesis related pathways as photosynthesis-dark reaction, chlorophyll biosynthetic process, and photosynthesis-light harvesting in photosystem I. Furthermore, some important carbohydrate-related pathways were enriched such as starch metabolic process, sucrose metabolic process, cellulose biosynthetic process, cell wall biogenesis, regulation of transcription, among others (Fig 8).
These pathways are of great interest since the most economically important characteristic of sugarcane is its biomass production and accumulation of sucrose [36].The results show negative effects of osmotic stress on biosynthesis of sucrose [37].Being a high water-demanding crop, the major limiting factors in sugarcane cultivation are abiotic stresses such as drought, salinity, and extreme temperatures [38].
Comparatively, data mining analysis for sugarcane done by using the Sugarcane Expressed Sequence Tags (SUCEST) database, revealed an abundant expression of genes encoding chaperones, co-chaperones and other proteins linked to protection against stress [46].
A well represented group of coding genes found in our transcriptome is the group of LEA proteins, commonly modulated by ABA [47].Some LEA proteins have been associated with diverse specific functions as anti-aggregation activity, membrane stabilization, ion scavenging, metal binding and water sequestration, among others [48].Along the kinetic, some members of this group of proteins were detected, for example LEA14A, Dehydrin DHN3, Dehydrin DHN4 and COR410 (see Fig 9).
Transcription factors are another interesting group of genes involved in osmotic stress that were enriched in our analysis.Among the up-regulated TFs in sugarcane, we detected three principal representative groups belonging to the HSF, NAC, and AP2/ERF families, previously reported as stress-responsive TFs [49].The transcript level of P5CS, an important gene of the proline biosynthesis pathway, was significantly altered upon exposure to PEG stress, in contrast with the reported by [50].However, induction of transcript expression of OsP5CS was reported in rice in response to high salt, dehydration, ABA and cold treatments [51].
Osmotic stress conditions also enhance the production of reactive oxygen species (ROS) that propel the cellular damage inducing oxidative stress.All pathways that activate mechanisms leading to ROS scavenging have been shown to play an important role in protecting plants against different abiotic stresses [52].Several genes coding detoxifying enzymes related to the response against oxidative stress are triggered by osmotic stress to mitigate the potential damage induced by ROS and other reactive.Other important groups of proteins involved in the redox metabolism are the thioredoxins, thiol-disulfide oxidoreductase enzymes participating in the cellular signalling pathways mediated by redox mechanism.
All these groups of proteins were subjected to regulation at transcriptional level in the experimental conditions of the present study, suggesting an induction of detoxifying systems in response to a parallel oxidative stress.
Additionally to the characterized and identified transcripts, a total of 20.2% of the 20,407 DEGs showed an unknown function (no match + uncharacterized protein).Li et al., (2016) determined by using microarray analysis that ~38.7% of DEG were assigned as proteins with unknown functions [30].In accordance, this finding remarks the importance of the study of proteins with assigned unknown function, as a benchmark for the discovery of new genetic sources for the understanding of molecular signalling during drought phenomenon.

Conclusions
The lack of genome and transcriptome data has restricted molecular studies in S. officinarum.Here, we provide a benchmark set of transcript profiles responding to osmotic stress.The knowledge gained from these studies play an important role on the grasp of the molecular processes involved in the osmotic stress response in sugarcane.The next challenge is to test those drought responsive genes in monocot plant models as rice as well as sugarcane itself, to evaluate the functionality of the selected genes.Finally, the development of stress tolerant sugarcane cultivars will be substantially beneficial for the production of this crop in drought predicted scenarios.
Osmotic stress was induced in in vitro-plantlets (90-days) which were transferred to magenta dishes containing liquid MS media and enriched with 40% of PEG 8000 (Sigma-Aldrich1) This PEG 8000 concentration corresponds to an osmotic potential of -4.5 MPa, measured with a digital hygrometer (Bapro1).

Tissue collection and total RNA extraction
RNA extraction from S. officinarum aerial tissues (stem and leaves) and roots was carried out using TRIzol reagent (Invitrogen1) according to the manufacturer's instructions.The integrity and quantity of the purified RNA were assessed both on 1% agarose gel and 2100 Bioanalyzer (Agilent Technologies1).Equal quantities of high-quality RNA from each sample were pooled for the library preparation.

Library construction, RNA-seq raw data cleaning, and assembly
Libraries were constructed and sequenced at the Genomic Services Laboratory, (LANGEBIO-CINVESTAV, Mexico).Libraries were prepared with TruSeq RNA Sample Prep Kit v2 (Illumina), and library quantification was performed with a 2100 Bioanalyzer (Agilent Technolo-gies1).For high-throughput sequencing, paired-end sequence samples were multiplexed with 16 libraries in one line, yielding ~410 million clean reads and with ~22 million reads per library for S. officinarum transcriptome.All libraries were sequenced in a 2x100 bp configuration on the HiSeq2500 Illumina platform.Raw reads were checked for quality with FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) [54], processed by clipping adapters, cleaned and filtered for quality scores greater than 30 and read length greater than 80 nucleotides using the Trimmomatic v0.36 toolkit [55].

Generating a reference transcriptome of sugarcane
The transcriptome reference was obtained from a RNA mix from leaves, stems and roots of plantlets of sugarcane var.Mex 69-290.RNA-seq was carried out on the complementary DNA libraries derived from non-stressed plantlets and groups of plantlets subjected to 24, 48, and 72 h of osmotic stress.Four libraries for leaves+stems (L0, L1, L2, L3; L = leaf) and four libraries for roots (R0, R1, R2, and R3; R = root) were sequenced with their respective biological replicates.A total, 412,440,102 raw reads were generated, and 410,382,788 clean reads with a Q>30 were retained.De novo assembly of all QC cleared reads was conducted with Trinity software v2.2.0 [56] and following the protocol reported in [57].The number of unigenes, N50 statistics, and the average coverage were used to evaluate the sequence assembly quality (see Table 1).
All unigenes were annotated by assigning their associated generic Gene Ontology (GO) terms and Enzymes codes (EC) using Blast2GO software [21].The annotation was based on the homology to proteins from other plant species determined by BLASTX.Plant ref-seq homology searches were complemented with information from INTERPROSCAN [58], the Kyoto Encyclopedia of Genes and Genomes (KEGG; [19] biochemical pathways, and expanded using ANNEX [59].BLAST searches of the unigenes were conducted against Swiss-Prot reviewed database (www.uniprot.org)containing 551,705 manual annotated proteins and TrEMBL database (www.uniprot.org)with 65,378,749 un-reviewed proteins.TransDecoder v3.0.0 software was used to identify candidate coding regions within unigene sequences generated by our de novo RNA-Seq transcript assembly, yielding a total of 94,416 open reading frames (ORFs) of at least 100 amino acids long.Detected ORFs were used for homology searches against Pfam-A v9.0 database [60] containing 5,724 families and 705,698 proteins, and another search was conducted against the Conserved Domain Database (CDD) [61] by the Batch CD-Search tool [62].The original KEGG map used for depict the starch and sucrose metabolism genes in Fig 3 is found in S5 File.

Differential expression analysis and GO annotation
Expression abundances were determined by mapping the pair-end read libraries (each replicate for each tissue) independently against our S. officinarum de novo transcriptome, using RSEM software v1.2.27 (RNA-Seq by Expectation-Maximization [20], implemented with bowtie aligner [63]).Expression was normalized to a transcripts per million transcripts (TPM), representing gene expression level.Nevertheless, the implemented script calculates both TPM and FPKM (Fragments per kilobase of exon per million fragments mapped) values.For a detailed explanation of all the steps and scripts implemented in this work (assembly, filtering, comparisons, gene clustering methods, and differential expression analysis) please see [57].The implemented pipeline is available at the github repository (https://github.com/trinityrnaseq/trinityrnaseq.wiki.git).
DEGs between treatments and tissues was determined by EdgeR [64] in R [65] following the pipeline described in [57].For each gene, a significant threshold of 0.001 was applied after the p-value was adjusted with false discovery rate (FDR) via Bonferroni-Holms correction [66].For the prediction analysis of gene homology, sugarcane DEGs were aligned to the Plant ref-seq proteins from NCBI with a 1e -3 E-value threshold by using BLASTX to annotate gene functions.Functional annotations were also created by Blast2GO software [21] with the GO and EC terms, and KEGG EC biochemical pathways.

Data and statistical analysis
Data analysis and graphics were performed with the R statistical package [65] unless stated otherwise.All the graphs were done by ggplot2 [67] and Vennerable R packages (https:// github.com/js229/Vennerable).To provide statistical support for physiological and biochemical analyses, we use the Tuckey's HSD (Honest Significance Difference; p<0.05) test for each treatment by using STATGRAPHICS Centurion XVI.I software and SigmaPlot 11.0 software.Figures editing were done by Inkscape v0.91 software (www.inkscape.org).

Gene validation by RT-qPCR
Five selected candidate genes were evaluated by quantitative real-time PCR.The designed primers used on the analysis are on S1 Table .The relative quantification of each transcript was calculated by the 2 -ΔΔC T method using ACT as internal control [68].

Determination of physiological parameters
Photosynthesis rate was determined with a portable Li-6400 photosynthetic system (LI-COR Biosciences1), using three samples per treatment.Plant water status was determined by measuring the osmotic potential on fully expanded leaves with a thermocouple psychrometer (HR-33T Dew point microvoltmeter, sample chambers type C-52; USA).RWC was computed using the formula in [69] RWC = [(FW-DW)/(TW-DW)] x 100, using the flag leaf of four plantlets.Proline content determination was performed for root and leaves according to [70].

Fig 1 .
Fig 1. Physiological and biochemical characterization of in vitro sugarcane plantlets under osmotic stress induced by PEG 8000.Four physiological parameters and one biochemical parameter were measured in sugarcane plantlets to characterize the effect of the osmotic stress imposition.A) Plantlets exposed to osmotic stress by PEG 8000 (40%).The time of stress exposure (from top to bottom) is indicated by an arrow showing the severity of the stress treatment increased according to time (0, 24, 48, and 72 hours respectively).Physiological measurement of B) Relative water content, C) Osmotic potential, D) Assimilation rate, E) Transpiration rate, and F) Biochemical measurement of proline content in roots and leaves of sugarcane.https://doi.org/10.1371/journal.pone.0189271.g001 95 in all cases (S1 Fig).Dynamic changes of gene expression patterns were retrieved by using a gathering cut-off value of fold change (FC) >2 and an adjusted p-value (FDR) of <0.001.On the basis of similar expression kinetics, 12,662 DEGs were depicted by a heatmap analysis to show major clusters of DEGs according to tissue specificity and time-stress treatment (Fig 5A).Similar expression profiles were observed for T0 samples and specific clusters were denoted for stressed samples (T1, T2, and T3) depending on plant tissues.Temporal expression patterns of genes following osmotic stress treatments are shown in Fig 5A where it is

Fig 2 .Fig 3 .
Fig 2. Summary of the annotated unigenes of S. officinarum var.Mex 69-290 transcriptome.A) Gene Ontology (GO) classification of the unigenes annotated.B) Distribution lengths of the assembled unigenes.C) Distribution of BLASTX-hits for annotated unigenes among the 10 principal plant species represented in the PlantRef protein database.D) Distribution of the closest homology matches (Top BLASTX-hits) for mapped unigenes in the ten principal species represented in the PlantRef protein database.E) Venn diagram is showing BLAST results of unique and shared unigenes among CDD, GO, Pfam, Sprot, TrEMBL, and PlantRef databases.F) Venn diagram showing shared unigenes among GO, Pfam, Sprot, and PlantRef databases.G) Venn diagram showing the BLAST results for the three principal protein databases, Pfam, Sprot, and CDD.Plant species names were coded for C) and D) as follows: Bd = Brachypodium distachyon, Eg = Elaeis guineensis, Ma = Musa acuminata subsp.Malaccensis, Ob = Oryza brachyantha, Os = Oryza sativa Japonica, Pd = Phoenix dactylifera, Sb = Sorghum bicolor, Si = Setaria italica, and Zm = Zea mays.https://doi.org/10.1371/journal.pone.0189271.g002

Fig 4 .
Fig 4. Correlation of sample replicates and hierarchical clustering analysis of gene expression levels in each sample (0, 24, 48, and 72 h) under osmotic stress treatment.A) Sample relatedness is plotted using the first two principal component (PCs) showing the variability between replicates and different treatments.Non-stressed leaf and root samples are in blue shade background.In dark-green and in red background colors are the samples corresponding to stressed leaves and roots (24, 48, 72 h ASI), respectively.C = control, D = Drought.B) Correlation dendrogram of gene expression in all samples.Correlation r values are coded into key colors where green indicates lower r values; reds indicates higher r above the mean.Blue squares in the nodes represent non-stressed samples, green stars represent leaf stressed samples, and the red stars represent root stressed samples.L = leaf, R = root, L0 and R0 = non-stress treatment, L1 and R1 = 24 h treatment ASI, L2 and R2 = 48 h treatment ASI, L3 and R3 = 72 h treatment ASI.* = biological replicate.https://doi.org/10.1371/journal.pone.0189271.g004 Fig 7A) and 2,319 down-regulated genes (Fig 7D).For T2 treatment (48 h ASI), 5,033 genes were differentially expressed, including 1,753 up-regulated transcripts (Fig 7A) and 3,280 down-regulated transcripts (Fig 7D), and finally, for T3 treatment (72 h ASI), 3,224 transcripts were found differentially expressed, which included 1,382 up-regulated transcripts (Fig 7A) and 1,842 down-regulated transcripts (Fig 7D).

Fig 5 .
Fig 5. Correlation and clustering analysis of the sugarcane DEGs subjected to osmotic stress treatments.A) Tissue-specific heatmap clustering of 12,662 significant DEGs (FDR 0.001 and FC !2) for the 0, 24, 48, and 72 h treatments ASI.B-G) Osmotic stress DEG profiles for the organ specific behavior clusters.Clusters were built according to their expression patterns during the stress treatments.Gray lines indicate the behavior of each gene into the treatment and the mean expression profile for each cluster is plotted as blue lines.Bottom number indicates the number of time series genes.L = leaf, R = root.L0 and R0 = non-stress treatment, L1 and R1 = 24 h treatment ASI, L2 and R2 = 48 h treatment ASI, L3 and R3 = 72 h treatment ASI.* = biological replicate.https://doi.org/10.1371/journal.pone.0189271.g005

Fig 6 .
Fig 6.Differentially expressed gene analysis along the osmotic stress treatments in sugarcane.A) Volcano plots show DEGs in each of the stress treatments for leaf and root tissues.The x-axis represents the log2FoldChange of genes, and the y-axis represents the statistical significance degree by-log10 (PValue).DEGs are shown in red dots and were detected using their statistical significance differences (FDR<0.001and a log2FoldChange>2).B) Comparison of the up-and down-regulated DEGs for leaf and root tissues in all three stress treatments and comparison of shared genes among all the treatments (CORE gene sets).https://doi.org/10.1371/journal.pone.0189271.g006

Fig 7 .
Fig 7. Comparative analysis of the unique and shared DEGs among the osmotic stress treatments for sugarcane var.Mex 69-290.The Venn diagrams show the overlap of the DEGs for leaves and roots submitted to a 3 time-series of osmotic stress treatments.A and C) Comparison among the genes differentially up-regulated for the stress treatments for leaf and root tissues respectively.D and F) Comparison among the genes differentially downregulated for the stress treatments for leaf and root tissues respectively.B and E) Comparison among the core set of up-and down-regulated genes respectively for both tissues.CULG = CORE of Up-regulated Leaf Genes, CURG = CORE of Up-regulated Root genes, COSUG = CORE of Osmotic Stress Up-regulated Genes, CDLG = CORE of Down-regulated Leaf Genes, CDRG = CORE of Down-regulated Root Genes, and COSDG = CORE of Osmotic Stress Down-regulated Genes.The Venn diagrams were done by using Vennerable R package.https://doi.org/10.1371/journal.pone.0189271.g007

Fig 8 .
Fig 8. Functional enrichment analysis for the CORE set of DEGs regulated in the osmotic stress treatments for sugarcane var.Mex 69-290.Significant groups (FDR<0.05)for COSUG and COSDG groups are plotted according to their enriched adjusted p-value (x axis).The size is proportional to the number of enriched genes among the total number of their GO term associated unigenes.COSUG = CORE of Osmotic Stress Up-regulated Genes for both leaf and root tissues, COSDG = CORE of Osmotic Stress Down-regulated Genes for both leaf and root tissues.https://doi.org/10.1371/journal.pone.0189271.g008

Fig 9 .
Fig 9. Principal findings of the analysis of sugarcane transcriptome under osmotic stress.A) Genes regulated under osmotic stress treatments.Inside blue square are the DEGs regulated in the three treatments.The green square represents the DEG regulated only in the T1 treatment, within the yellow square are the DEGs expressed only in T2 treatment, and within the red square are the DEGs regulated only in T3 treatment.Inside gray circles are DEGs shared among treatments.Letters in black color represent the up-regulated DEGs and letters in orange color represent the down-regulated DEGs per treatment.B) Schematic depiction of the physiological results obtained in sugarcane var.Mex 69-290, imposed to an osmotic stress treatment by PEG 8000 for 24, 48, and 72 h.The green arrow under the depiction indicates the directions of the time-stress treatments.RCW = relative content water, OsPot = osmotic potential, AR = assimilation rate.Complete name of the genes in A) are listed in the S4 File.https://doi.org/10.1371/journal.pone.0189271.g009