In this study, plant materials from two sugarcane Saccharum L. cultivars, F172 and Guitang 31 (GT31), with strong and weak drought tolerance, developed by the Taiwan Sugar Research Institute (Taiwan, China) and Sugarcane Research Institute, Guangxi Academy of Agricultural Sciences (Nanning China), respectively, were used (Yang and Li, 1992; Li et al., 2011).

The sugarcane seedlings were observed and photographed at 4, 5, and 7 d after water withdrawal and were categorized under different drought treatment conditions: mild drought stress (B), moderate drought stress (C), and severe drought stress (D) ( Figure 1 ). Both varieties showed significant changes across the three time points, and cultivar F172 showed stronger drought tolerance than Guitang31 ( Figure 1 ). Sugarcane leaf tissues from cultivars F172 and GT31 were obtained under B, C, and K drought conditions for subsequent transcriptomic sequencing and photosynthetic rate detection. In parallel, leaf tissues from the two cultivars under normal watering conditions were collected at the same time points as controls for subsequent transcriptomic sequencing and detection of the photosynthetic rate, which are abbreviated as BCK (on day 4), CCK (on day 5), and DCK (on day 7).

The experiment was performed in triplicates for each condition and cultivar. All the leaf samples were flash-frozen with liquid nitrogen and stored at −80°C until further use.

RNA was extracted using the RNeasy Mini Kit (Qiagen, Beijing, China), followed by purification, fragmentation, and quality control using a NanoDrop 2000 (Thermo Scientific) and an Agilent 2100 Bioanalyzer (Agilent Technologies). Strand-specific libraries were obtained using dUTP for second-strand synthesis and subsequently sequenced on a BGISEQ-500 instrument (BGI, Shenzhen, China). The experiments were conducted in triplicates for each cultivar at each time point under drought (B, C, and K) and control (BCK, CCK, and DCK) conditions.

Preprocessing of the paired-end reads was performed using FASTP (v 0.23.1), and mapped to the sugarcane genome (NCBI accession: ASM2245720v1) using HISAT2 (v 2.20) (Kim et al., 2019). Raw read counts were quantified using the featureCounts software (SUBREAD v2.0.0) (Liao et al., 2014). Differential expression analysis was performed using the DESeq2 package (v1.38.3) (Love et al., 2014) for transcriptome comparisons between cultivars under the same conditions (GT31-B-vs-F172-B, GT31-C-vs-F172-C, GT31-D-vs-F172-D, GT31-BCK-vs-F172-BCK, GT31-CCK-vs-F172-CCK, and GT31-DCK-vs-F172-DCK) and between treatments and the corresponding controls for each cultivar (GT31-BCK-vs-GT31-B, F172-BCK-vs-F172-B, GT31-CCK-vs-GT31-C, F172-CCK-vs-F172-C, GT31-DCK-vs-GT31-D, and F172-DCK-vs-F172-D) to identify the differentially expressed genes (DEGs) with an absolute value of log2 FC > 1.0 and false discovery rate < 0.05. KEGG pathway analyses of the identified DEGs were conducted using ClusterProfiler (v4.3.1) (Yu et al., 2012).

Clusters of genes with the same expression profile over different time points were identified using the short time-series expression miner (v1.3.13) (Ernst and Bar-Joseph, 2006) for cultivar F172 under drought stress (F172), cultivar F172 controls (F172CK), cultivar GT31 under drought stress (GT31), and cultivar GT31 controls (GT31CK). The statistical significance of the number of genes for each profile compared to the expected number was computed using a permutation-based test. Unique and common trend genes for significant profile clusters (p < 0.05) from the above-mentioned four groups were selected and KEGG enrichment analysis was performed as described above.

Weighted correlation network analysis (v 1.69) (Langfelder and Horvath, 2008) was used to infer the network modules (parameters: softPower = 20, mergeCutHeight = 0.7, minModuleSize = 30)for 58,873 genes after filtering those with low expression levels (Fragments per kilobase of transcripts per million fragments mapped < 0.5). Module-trait associations were estimated using the correlation between module eigengenes and traits, including cultivars, drought conditions, and the strength of drought stress. Module–trait associations were considered statistically significant at p < 0.05. Trait-related genes with significant correlations were extracted from the module and subjected to KEGG enrichment analysis as described above. Hub genes in the modules were identified using the CytoHubba module in Cytoscape software.

A portable photosynthesissystem (Li-6800, Li-COR Biosciences, Lincoln, NE, USA) was used to observe the net photosynthetic rate for the functional topvisible dewlap leaf (leaf + 1) of sugarcane. as previously described (Verma et al., 2020). The photosynthesis rate was measured with three biological replicates for each cultivar at each time point under both drought treatment and control conditions.

The pheatmap package (v1.0.12) (https://CRAN.R-project.org/package=pheatmap) was used to plot the heatmaps. Differences were calculated using the t-test and were considered statistically significant at p < 0.05. Photosynthetic rate data were statistically analyzed using one-way ANOVA followed by Duncan’s multiple range test.

RNA sequencing generated a total of 379.6 G of raw data for all 36 samples (each in the range of 7.8 G–13.9 G) ( Table S1 ). Approximately 373.8 G clean reads were obtained and passed through quality control; all samples were of high quality (Q20 ≥ 96.68% and Q30 ≥ 91.86%) ( Table S1 ). The mapping rate of clean reads to the sugarcane reference genome ranged from 76.16% to 79.84% ( Table S1 ). Principal component analysis showed low inter-replicate variability, and the samples in each group clustered together ( Figure 2A ). The principal component analysis clearly distinguished between the water deficit and control conditions, and the cultivars were also well separated ( Figure 2A ), indicating a large variability between the F172 and GT31 cultivars.

Differential gene expression analysis was performed between both cultivars grown under the same conditions, and between those grown under different conditions along with the corresponding controls ( Figure 2B ; Table S2 ). A total of 8,546, 10,036, 13,517, 11,776, 14,564, and 5,506 DEGs were identified between the GT31 and F172 cultivars under the same drought conditions, including B, BCK, C, CCK, D, and DCK, respectively ( Figure 2B ). On comparing intra-cultivar drought treatments and controls, we identified 33,058, 34,197, and 39,741 DEGs from the comparisons of three conditions (B vs. B-vs-BCK, C vs. C-vs-CCK, and D-vs-DCK) for strain GT31, and 30,618, 35,060, and 38,763 DEGs for strain F172, respectively ( Figure 2B ). In general, the number of DEGs between the intra-cultivar drought treatments and controls was much greater than that between the inter-cultivar differences under the same conditions. These DEGs between cultivars grown under the same conditions were mainly enriched in starch and sucrose metabolism, linoleic acid metabolism, glycolysis/gluconeogenesis, glyoxylate and dicarboxylate metabolism, nitrogen metabolism, carotenoid biosynthesis, and ribosome pathways among others ( Figure 2C ). Compared with the DEGs between cultivars under normal watering conditions, those under drought conditions enriched in pyruvate metabolism, beta-alanine metabolism, glutathione metabolism specifically under moderate and severe drough stress ( Figure 2C ). In the comparison groups BCK-vs-GT31-B, CCK-vs-GT31-C, and DCK-vs-GT31-D, we identified 33,058 DEGs (12,697 up and 20,361 down), 34,197 DEGs (14,328 up and 19,869 down), and 39,741 DEGs (15,821 up and 23,920 down) for strain GT31, and 30,618 DEGs (12,141 up and 18,477 down), 35,060 DEGs (14,444 up and 20,616 down), and 38,763 DEGs (15,876 up and 22,887 down) for stain F172, respectively ( Figure 2B ). The DEGs between cultivars grown under different conditions along with the corresponding controls were mainly enriched in starch and sucrose metabolism, linoleic acid metabolism, glycolysis/gluconeogenesis, glyoxylate and dicarboxylate metabolism, valine, leucine and isoleucine degradation, phenylpropanoid biosynthesis, flavone and flavonol biosynthesis, photosynthesis-antenna proteins, glycosphingolipid biosynthesis (globo- and isoglobo- series), and phagosome pathways among others ( Figure 2C ). Many DEGs were identified between the drought stress and control groups for the same time points for both cultivars ( Figure 2B ; Table S2 ), indicating that drought stress has a significant impact on gene expression.

Time-series expression analysis was performed to examine the dynamic transcriptomic differences in drought tolerance between cultivars. The gene expression patterns for cultivars grown under different conditions and their corresponding controls were analyzed for both cultivars across different drought time points, and the profiles with p < 0.05 in the permutation test were considered as significant ( Figure 3 ). Profiles 0, 1, and 6 were considered to be significant for cultivar F172 in drought stress ( Figure 3A ; Table S3 ); profiles 4 and 3 were considered to be significant for cultivar F172 controls ( Figure 3B ; Table S4 ); profiles 1, 6, and 0 were considered to be significant for cultivar GT31 under drought stress ( Figure 3C ; Table S5 ); and profiles 1 and 6 were considered to be significant for cultivar GT31 controls ( Figure 3D ; Table S6 ).

To identify the trend genes specific to F172 in response to drought stress, an intersection analysis was conducted for the trend genes in different groups, including the cultivar F172 under drought stress (F172), cultivar F172 controls (F172CK), GT31 under drought stress (GT31), and GT31 controls (GT31CK) ( Figure 4A ). Among them, 5,103 genes exhibited a specific trend in the F172 group that was inconsistent with those in the GT31 and the F172 controls, potentially related to the drought tolerance of F172. The KEGG pathway enrichment results showed that these genes were significantly enriched in photosynthesis, mitogen-activated protein kinases (MAPK) signaling pathway, biosynthesis of various plant secondary metabolites, and cyanoamino acid metabolism pathways ( Figure 4B ).

We employed WGCNA to identify potential co-expression modules and key regulatory networks, and further elucidate their roles in the response of the F172 cultivar to drought stress. WGCNA categorized all genes into 15 modules ( Figure 5A ). Module-trait relationships were explored to extract significant associations between cultivars, drought conditions, strength of drought stress, and modules ( Figure 5B ). The blue4 and plum1 modules showed significant negative (–0.86, p < 0.01) and positive (0.96, p < 0.01) correlations, respectively, with drought conditions, ( Figure 5B ). The tan and salmon4 modules showed significant negative (-0.93, p < 0.01) and positive (0.8, p < 0.01) correlations, respectively, with the cultivars ( Figure 5B ).

All unique trend genes ( Figure 4A ) were mapped to weighted correlation network analysis modules ( Table S7 ). Among the unique trend genes identified in F172 compared to GT31 and controls, most were in the blue4 (1,301) and plum1 (983) modules ( Table S7 ); these were then used for KEGG enrichment analysis. Genes in the blue4 module were mainly enriched in phagosome, endocytosis, photosynthesis-antenna proteins, amino sugar and nucleotide sugar metabolism, and oxidative phosphorylation, among others ( Figure 5C ; Table S8 ). Those genes in the plum1 module were mainly enriched in spliceosome, RNA transport, aminoacyl-tRNA biosynthesis, basal transcription factors, peroxisome, and alanine, aspartate, and glutamate metabolism, among others ( Figure 5D ). Several enriched pathways, including those of photosynthesis, purine metabolism, starch and sucrose metabolism, beta-alanine metabolism, photosynthesis-antenna proteins, and plant hormone signal transduction, were shared between those genes in the blue4 module and the unique trend genes in F172 ( Figure 5C ; Table S8 ). The hub genes identified in the blue4 module included ARF14 (auxin response factor 14) and Os10g0147400 (similar to auxin influx carrier protein), which are associated with the plant hormone signal transduction pathway, and Os02g0733300 (glycoside hydrolase gene), which is associated with the starch and sucrose metabolism pathway ( Table S9 ). The hub genes of plum1 module included the spliceosome-related genes RS31 (serine/arginine-rich splicing factor RS31) and CLPC1 (chaperone protein ClpC1, chloroplastic), and peroxisome-related genes PEX1 (peroxisomal biogenesis factor 1) and PEX2 (peroxisomal biogenesis factor 2) ( Table S10 ).

Photosynthesis is closely associated with drought stress in plants ( Figure 6 ). In the photosynthesis-antenna protein pathway for F172, the expression levels of Lhca1 (light-harvesting complex I chlorophyll a/b binding protein 1) and Lhca4 (light-harvesting complex I chlorophyll a/b binding protein 4) were high at stage 1 and significantly decreased at later stages, exhibiting a distinct downward trend ( Figure 6 ). Further, PsbO (photosystem II oxygen-evolving enhancer protein 1), PsbP (photosystem II oxygen-evolving enhancer protein 2), PsbQ (photosystem II oxygen-evolving enhancer protein 3), PsbW (photosystem II PsbW protein), PsbY (photosystem II PsbY protein), Psb27 (photosystem II Psb27 protein), PsaD (photosystem I subunit II), PsaE (photosystem I subunit IV), PsaF (photosystem I subunit III), PsaH (photosystem I subunit VI), PsaL (photosystem I subunit XI), PsaN (photosystem I subunit PsaN), PetF (ferredoxin), and petH (ferredoxin–NADP+ reductase), which participate the photosynthesis pathway, also showed a similar expression pattern for F172 ( Figure 6 ).

The photosynthetic rates for two cultivars F172 and GT31 were significantly down-regulated under drought conditions. However, this decrease was more pronounced for the F172 variety, which was consistent with the declining trend observed in the expression of genes related to the photosynthetic pathway as described above ( Figure 7 ). These results indicate that drought stress resulted in a significantly greater reduction in photosynthesis for cultivar F172, which is associated with stronger drought resistance.