Leaf transcriptome profiling of contrasting sugarcane genotypes for drought tolerance under field conditions

Drought is the most detrimental abiotic stress to sugarcane production. Nevertheless, transcriptomic analyses remain scarce for field-grown plants. Here we performed comparative transcriptional profiling of two contrasting sugarcane genotypes, ‘IACSP97-7065’ (drought-sensitive) and ‘IACSP94-2094’ (drought-tolerant) grown in a drought-prone environment. Physiological parameters and expression profiles were analyzed at 42 (May) and 117 (August) days after the last rainfall. The first sampling was done under mild drought (soil water potential of −60 kPa), while the second one was under severe drought (soil water potential of −75 kPa). Microarray analysis revealed a total of 622 differentially expressed genes in both sugarcane genotypes under mild and severe drought stress, uncovering about 250 exclusive transcripts to ‘IACSP94-2094’ involved in oxidoreductase activity, transcriptional regulation, metabolism of amino acids, and translation. Interestingly, the enhanced antioxidant system of ‘IACSP94-2094’ may protect photosystem II from oxidative damage, which partially ensures stable photochemical activity even after 117 days of water shortage. Moreover, the tolerant genotype shows a more extensive set of responsive transcription factors, promoting the fine-tuning of drought-related molecular pathways. These results help elucidate the intrinsic molecular mechanisms of a drought-tolerant sugarcane genotype to cope with ever-changing environments, including prolonged water deficit, and may be useful for plant breeding programs.


Results
Physiological responses to drought. The impacts of drought stress on sugarcane plants were investigated by analyzing physiological parameters. As 'IACSP97-7065' presented decreases in CO 2 assimilation rate (41%) and photochemical efficiency (8%) after 42 days of the last rainfall (Fig. 1c,e), it can be argued that this genotype was facing mild drought stress. On the other hand, 'IACSP94-2094' did not show any changes in the same parameters analyzed at the same sampling point. After 117 days after the last rainfall, both 'IACSP97-7065' and 'IACSP94-2094' genotypes showed reductions in total chlorophyll content (16% and 30%, respectively) and CO 2 assimilation (41% and 67%, respectively) ( Fig. 1a-c), indicating a severe drought stress condition. However, only 'IACSP97-7065' exhibited an increase in non-photochemical quenching processes (260%), and decreases in potential (8%) and effective (50%) quantum efficiencies of PSII ( Fig. 1d-f).
Analysis of deregulated pathways under drought. Differentially expressed genes in 'IACSP97-7065' and 'IACSP94-2094' were classified according to KEGG pathway categories (p ≤ 0.05) for mild and severe drought stresses (Table 1). Under mild drought stress, 'IACSP97-7065' displayed DEGs associated with metabolism of energy, amino acids, carbohydrates, lipids, and environmental adaptation; conversely, the drought-tolerant genotype 'IACSP94-2094' showed no enriched pathways. Moreover, KEGG pathway categories were substantially enriched for both genotypes under severe stress. Besides the induced pathways due to mild drought stress, both genotypes also showed activities associated with the metabolism of other amino acids, and translation and transcription pathways. Regarding the common plant slim categories between both genotypes, enriched DEGs in 'IACSP94-2094' genotype often outnumbers those of 'IACSP97-7065' . Noteworthy, genotype-specific plant slim categories were also revealed. DEGs of the drought-sensitive genotype 'IACSP97-7065' were classified into exclusive categories, such as citrate cycle, fructose and mannose metabolism, pyruvate metabolism, glyoxylate and dicarboxylate metabolism, amino sugar and nucleotide sugar metabolism, glutathione metabolism, cyanoamino acid metabolism, selenocompound metabolism, carbon fixation in photosynthetic organisms, fatty acid elongation, and arachidonic acid metabolism. On the other hand, the drought-tolerant genotype 'IACSP94-2094' comprises DEGs that are represented by photosynthesis (antenna proteins), plant-pathogen interaction, spliceosome and metabolism of amino acids (i.e., alanine, aspartate, glutamate, tyrosine), glycerophospholipid, and nitrogen. The identification of genes involved in each enriched KEGG pathway category is shown in Supplementary Table S3 online.
Validation of microarray results. The correlation analysis of DEGs detected by microarray and validated by RT-qPCR indicated a strong similarity (Pearson's r = 0.849, p-value = 0.001; Spearman's rank = 0.834; p-value < 0.0001) for the gene expression data (see Supplementary Fig. S1 online), which supports the reproducibility and accuracy of the transcriptional profiles presented in this study. As shown in Fig. 4, RT-qPCR analysis confirmed the differential expression of 18 genes in at least one of the sugarcane genotypes. In common, only the FLA16 gene was found differentially expressed for both genotypes under mild stress. Moreover, six genes were significantly differentially expressed in both genotypes under severe stress, among which LOX2, DHN1, and NAM were up-regulated, whereas MYB1, CesA, and PHO1 were down-regulated.

Discussion
Water is essential for plant development, and the sugarcane formative phase (i.e., grand growth stage) is critically dependent on water supply 19 since the shortage can reduce cane yield by up to 70% 20 . In a drought scenario, physiological, biochemical, and molecular responses are triggered specifically to the organ type 21 . In general, low water availability is detrimental to photosynthesis, leading to low functioning of photosystem II (PSII) 22,23 , which turns the maximum quantum efficiency of PSII (F v /F m ) into an indicator of drought stress 24 . Moreover, in the water-limiting context, the absorbed light energy exceeds its use in photosynthetic reactions, which results in energy dissipation by non-photochemical processes 25 . Accordingly, physiological analyses in leaves of droughtstressed sugarcane cultivars ('IACSP97-7065' and 'IACSP94-2094') revealed several changes in photosynthetic variables (Fig. 1), with long term exposure to low water availability intensifying plant responses. However, only 'IACSP97-7065' was promptly affected at 42 days after the last rainfall and showed impacts on carbon assimilation and F v /F m , indicating mild drought stress. In addition, this sensitive genotype showed increases in the nonphotochemical quenching and reduced potential and effective quantum efficiencies of PSII at 117 days after the last rain (Fig. 1d-f), which revealed plants were facing severe drought stress. These results suggest that excessive light energy was partially dissipated by non-photochemical processes (NPQ) and caused damage to PSII. On the other hand, 'IACSP94-2094' showed no decreases in F v /F m and Φ PSII , indicating that the maximum conversion of light energy into chemical energy was not affected by drought stress. Interestingly, such photosynthetic response of 'IACSP94-2094' to water deficit was previously highlighted under greenhouse conditions 26 www.nature.com/scientificreports/ its drought tolerance when cultivated in a drought-prone field. Curiously, the 'photosynthesis-antenna proteins' pathway was enriched only for 'IACSP94-2094' under severe drought (Table 1), representing the repressions of light-harvesting chlorophyll-binding proteins encoding genes. We would argue that this may be a molecular strategy to avoid excessive energetic pressure in PSII 28 and partially justify the maintenance of F v /F m and Φ PSII even after 117 days of water deficit. Transcriptional responses comprise multiple sets of drought-responsive genes that are spatiotemporally repressed or induced in an intricate manner 9 . Global gene expression regulation has already been investigated in sugarcane plants under drought conditions [10][11][12][13][14][15][16][17][18] . However, few studies have addressed the contrasting transcriptional profile between drought-tolerant and drought-sensitive sugarcane genotypes 11,13,15,17,18 ; even rarer are those  www.nature.com/scientificreports/ conducted within natural settings, exploring the dynamic influence of the environment 15 . In this study, we used microarray technology to decipher the distinct molecular adaptations of contrasting sugarcane (Saccharum spp.) genotypes, 'IACSP97-7065' (drought-sensitive) and 'IACSP94-2094' (drought-tolerant), to long-term water-deficit stress under natural environmental conditions. From the analysis of 14,552 SAS, we identified a set of 662 (4.6%) DEGs in sugarcane leaves under mild and severe stress (Fig. 2a). Among those, most are not shared between the sugarcane genotypes, which reinforces their divergent responses to drought at the transcriptional level (Fig. 2b,c). Yet, several DEGs found in both sugarcane genotypes are similar to the major drought stress-related hub genes (DSRhub genes) from sorghum (Sorghum bicolor) regulatory networks 29 , including transcription factors and other regulatory proteins. In particular, the lower number of DEGs in the 'IACSP94-2094' genotype under mild stress (fivefold less than 'IACSP97-7065') may be attributed to its drought-tolerant profile, in which a mild drought is unlikely to cause cell damage and/or the activation of drought-responsive signaling networks. Moreover, the functional enrichment analyses revealed that the DEGs of the sensitive and tolerant sugarcane genotypes are associated with different cellular processes (Fig. 3). Although hormone signaling pathways were not significantly enriched, some related genes already reported in drought-stressed sugarcane 30 and sorghum 29 were genotype-specific regulated. For instance, whereas only 'IACSP97-7065' shows the ethylene biosynthesis pathway induction, our RT-qPCR analysis revealed the auxin-responsive protein (IAA12-SCJFRT1007H07.g) as up-regulated in only 'IACSP94-2094' under severe drought stress (Fig. 4b). Collectively, the microarray data from this study may indicate distinct strategies of those plants in coping with drought. www.nature.com/scientificreports/ Environmental stresses such as drought are conservatively sensed by plants, leading to the overproduction of reactive oxygen species (ROS) as key components of signal transduction pathways; however, with detrimental effects to cells 31 . The oxidative stress in sugarcane causes significant damage to PSII and can virtually abolish photosynthetic efficiency and electron transport rate 32 . Under mild stress, despite the lower number of DEGs in the 'IACSP94-2094' genotype, some were uniquely associated with redox enriched pathways. Those included polyphenol oxidase (SCSGLR1045A10.g), cinnamyl alcohol dehydrogenase (SCSBAD1054A07.g), chalconeflavanone isomerase 1 (SCCCLR1C03G09.g), and cytochrome P450-related genes (SCCCCL3001H12.g), which have been reported to participate in acclimation to water and oxidative stresses in sugarcane [33][34][35] , sorghum 29 and model species 36 . In special, induction of CYP450 gene expression may represent an early mechanism for the tolerant genotype to avoid oxidative injury induced by high levels of hydrogen peroxide, as reported in droughtstressed CYP450-RNAi knockdown cotton plants 37 . Remarkably at severe drought stress, the tolerant genotype showed a greater number of deregulated genes in oxidoreductase activities than the sensitive genotype, with 20 exclusive DEGs. Additionally, some common oxidoreductase-related DEGs may exhibit relevant opposite regulations on gene expression of sugarcane genotypes. For example, glutathione s-transferase (GSTU) is strongly associated with preventing oxidative burst in plant cells exposed to abiotic stress 38 . Interestingly, our RT-qPCR results confirmed that GSTU was down-regulated in the sensitive genotype (log 2 FC = −1.12) and up-regulated in the tolerant genotype (log 2 FC = 1.62) under severe drought stress. In a greenhouse study, 'IACSP94-2094' exhibited the antioxidant system as a key feature for the recovery of the photosynthetic system after exposure to drought stress, thus enabling normal plant growth 7 . Overall, these results reinforce that the tolerant genotype 'IACSP94-2094' harbors an enhanced antioxidant system that might alleviate drought-induced damage.
Drought tolerance in plants is ensured by signaling pathways that elicit critical players in gene expression modulation-the transcription factors (TFs) 39 . Despite the lack of enriched transcriptional regulators in the 'IACSP94-2094' genotype under mild stress, we identified that R2R3-MYB was down-regulated by RT-qPCR analysis. We previously reported that overexpression of the most abundant alternative form of ScMYBAS1 transcript (ScMYBAS1-2) resulted in reduced biomass in transgenic rice plants subjected to drought 40 . Accordingly, the repression of R2R3-MYB was also detected in another tolerant sugarcane genotype (KPS01-12) under drought stress 17 , suggesting this down-regulation as a beneficial trait when plants face drought. Under severe drought stress, the enriched TF genes in 'IACSP94-2094' outnumbered those of the sensitive genotype, unveiling particular up-(MADS-box 1, SCCCCL4015C03.g; auxin-responsive protein IAA12, SCEQRT2093D08.g; and RNA polymerase II transcriptional coactivator KELP, SCCCLR2001H06.g) and down-regulated genes (a bZIP family member TGAL3, SCSGFL4194F09.g and nuclear transcription factor Y subunit A-10, SCRLLR1109E12.g). Most of these genes have already been reported in plant responses to abiotic stresses 29,[41][42][43] . For example, A. thaliana transgenic lines overexpressing TaNF-Y10 were more sensitive to water and salt stress, as compared to wild-type plants 41 . Furthermore, overexpression of TaMADS51 improves growth, biomass yield, and phosphate accumulation in tobacco plants, as well as the antioxidant enzymatic activities, specifically under phosphorusstarvation conditions 42 . Interestingly, a network-based study in sorghum revealed the responsiveness of MADSbox and bZIP transcription factors to co-occurring stresses, including drought 29 , which may cross-talk with some transcriptional responses of sugarcane found herein. Thus, our results indicate that the drought tolerance of the 'IACSP94-2094' genotype may be partially explained by an enhanced transcriptional regulation activity, which has conserved elements with sorghum drought-related pathways.
Carbohydrate and amino acid metabolism act synergistically in cellular acclimation to environmental stresses [44][45][46] . The enriched carbohydrate metabolism pathway comprises genes that are most significantly repressed in both genotypes under drought stress; except for trehalose-6-phosphate synthase (T6P, SCCCCL2001H04.b), which plays a role in osmoprotection and dramatically improves drought tolerance when overexpressed in model plants 47,48 . Interestingly, the tolerant genotype 'IACSP94-2094' exhibits a particular altered amino acid metabolism under severe drought stress, relying on up-regulated genes, such as NADH-dependent glutamate synthase 1 (GLT1, SCCCCL3120F01.g), glutamate decarboxylase 1 (GAD1, SCCCRZ2C03D05.g), and aspartate kinase-homoserine dehydrogenase (AK-HseDH, SCACCL6010C02.g). GLT1 and GAD1 genes are involved in glutamine synthetase/glutamate synthase (GOGAT/GS) cycle and the γ-Aminobutyric acid (GABA) metabolism, which are associated with the maintenance of the tricarboxylic acid (TCA) cycle and nitrogen assimilation in wild tomato (Solanum pennellii), serving as key players in drought tolerance 49 . These results suggest the metabolism of amino acid and nitrogen as an additional layer of the adaptive molecular machinery of 'IACSP94-2094' genotype.
Natural field conditions harbor multiple biotic and abiotic stresses that trigger overlapping molecular networks in plants to withstand limiting environments 52 . For instance, some of the DEGs of both sugarcane genotypes are strictly associated with pathogen resistance, such as receptor-like serine/threonine-protein kinase (ALE2, SCCCLR1079G06.g). Therefore, not every DEG found in this study can be readily anticipated to be www.nature.com/scientificreports/ related to drought. Moreover, we also detected around 5% of hypothetical drought-induced proteins in leaves of 'IACSP94-2094' under severe drought stress. Both findings require further investigations, using gene editing/ knockdown and/or gene overexpression strategies, combined with structural bioinformatics (structural modeling and molecular docking) to identify precisely the functions of those proteins. To illustrate, a putative 32.7 kDa jasmonate-induced gene (SCJLLR1103A10.g)-currently known as dirigent-jacalin (ShDJ)-was exclusively up-regulated in 'IACSP94-2094' under severe drought. Further, it was overexpressed in transgenic rice lines and improved drought tolerance, plant biomass accumulation, and saccharification efficiency 53 . Finally, this study revealed the gene expression profiles underlying the differential responses of contrasting sugarcane genotypes for drought tolerance under field conditions. After 117 days of the last rainfall, the tolerant genotype displayed significant changes in biological pathways associated with transcriptional regulation, oxidoreductase activity, amino acid metabolism, and translation process. Collectively, these molecular traits portray the molecular basis of 'IACSP94-2094' to cope with severe drought stress. For example, (i) the high oxidoreductase activity might protect PSII against oxidative stress, thus favoring photosynthetic activity even under water shortage 32 ; (ii) the enhanced amino acid metabolism may regulate the balance of compatible osmolytes, signaling molecules, and secondary metabolites precursors, which are critical in abiotic stress responses 54 ; and (iii) the remarkable transcriptional regulation activity may orchestrate the expression of multiple genes, including those associated with drought tolerance 55 . Thus, our study contributes substantially to the understanding of drought-induced response networks in sugarcane leaves at a trustworthy water-limiting scenario, as well as the discovery of new putative genes, which together may provide novel strategies for plant breeding programs. Physiological analyses. The physiological responses of 'IACSP97-7065' and 'IACSP94-2984' were assessed under irrigated (control) and non-irrigated (treatment) conditions at mild and severe drought stress stages. Chlorophyll fluorescence was measured with a fluorometer coupled to an infrared gas analyzer (model LI-6400XT, LI-COR, Lincoln, NE, USA) and some photochemical parameters were estimated: the maximum (F v /F m ) and effective (Φ PSII ) quantum efficiency of PSII, and the non-photochemical quenching (NPQ). While F v /F m is based on the minimum (F o ), maximum (F m ) and variable (F v = F m − F o ) fluorescence signals in dark-adapted leaf tissues, Φ PSII is based on the steady-state (F s ), maximum (F m ′) and variable (ΔF = F m ′ − F s ) fluorescence signals in light-adapted leaves, following the saturation pulse method (λ < 710 nm, Q ~ 12,000 µmol m −2 s −1 , 0.8 s) 56 . Chlorophyll a and b contents were also indirectly measured with a chlorophyll-meter (clorofiLOG, Falker, Brazil), by calculating the Falker Chlorophyll Index (FCI), which represents the chlorophyll mass/leaf mass ratio. In addition, the CO 2 assimilation (P n ) was evaluated, as a gas exchange-associated parameter, by using an infrared gas analyzer (model LI-6400XT, LI-COR, Lincoln, NE, USA). These measurements were performed on sugarcane leaves (the first fully expanded leaf with visible ligule, known as leaf + 1) under photosynthetic photon flux density of 2000 µmol m −2 s −1 ) and air CO 2 concentration of 385 ± 2 ppm, considering the natural variation of air temperature (34.7 ± 1.0 °C) and relative humidity (41 ± 2%). Measurements were taken between 09:00 and 11:00 am, when leaf temperature was 34.6 ± 1.0 °C. Technical and biological triplicates were subjected to the three-way analysis of variance (ANOVA), followed by Fisher's least significant difference test, which was performed using GraphPad Prism version 8.0.0 for Windows, GraphPad Software, San Diego, California USA (www. graph pad. com).

Methods
Microarray procedures and data analysis. Total RNA was extracted from 1 g of leaf tissue by the lithium chloride (LiCl) method and purified using the RNeasy kit (QIAGEN, Hilden, Germany) according to the manufacturer's instructions. Isolated RNA was quantified using a Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) spectrophotometer, and its integrity was evaluated via 1.2% agarose electrophoresis. Then, RNA samples were treated with the DNase I RNase-free (Fermentas, Vilnius, Lithuania).
The oligo microarray (4 × 44 K) contained 14,552 Sugarcane Assembled Sequences (SAS) released from Sugarcane EST Project (SUCEST) database 57  www.nature.com/scientificreports/ hybridization of the samples were performed according to the Two-Color Microarray-Based Gene Expression Analysis protocol. Hybridized microarrays were scanned using the Agilent Microarray Scanner (Agilent Technologies, Santa Clara, CA, USA). Thus, two biological and technical replicates were used for each microarray experiment. Microarray data sets can be found at the Gene Expression Omnibus (GEO) (series, record GSE187416) and at the SUCEST database (http:// sucest-fun. org/). The statistical analyses were carried out using the HTself method 58 to determine the differentially expressed genes (DEGs) in each experiment (confidence level set to 0.9). For each experimental condition, DEGs were determined when at least 70% of the spots showed similar and positively correlated patterns in the two technical replicates. Finally, the gene expression was calculated and normalized in log 2 Fold-change (FC). To confirm the reliability of the microarray data, the Pearson's and Spearman's rank correlation coefficients (p-value ≤ 0.05) were determined using the log 2 FC values obtained from the microarray and RT-qPCR analyses. Pearson's and Spearman's correlation analyses were performed using Graph-Pad Prism version 8.0.0 for Windows, GraphPad Software, San Diego, California USA (www. graph pad. com). We randomly selected 20 DEGs for the RT-qPCR assay to confirm the gene expression patterns. Each primer pair is described in Supplementary Table S4 online. Melting curve analysis between 72 and 95 °C was performed to confirm the specificity of the reaction. For each sample, three biological and technical replicates were used (n = 3). The reference genes rpl35-4, rRNA1, and ubq1 (see Supplementary Table S4 online) were used to normalization of the relative expression values. Finally, relative expression levels (drought treatment against irrigated control) were calculated using the comparative C t method 65 , followed by log 2 normalization.

Data availability
The datasets generated and analyzed during the current study are available in the Gene Expression Omnibus repository (Accession code: GEO series record GSE187416), https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE18 7416.