Comparative Transcriptome Analysis Reveals Potential Gene Modules Associated with Cold Tolerance in Sugarcane (Saccharum officinarum L.)

Sugarcane is an important crop worldwide, and most sugar is derived directly from sugarcane. Due to its thermophilic nature, the yield of sugarcane is largely influenced by extreme climate conditions, especially cold stress. Therefore, the development of sugarcane with improved cold tolerance is an important goal. However, little is known about the multiple mechanisms underlying cold acclimation at the bud stage in sugarcane. In this study, we emphasized that sensitivity to cold stress was higher for the sugarcane variety ROC22 than for GT42, as determined by physical signs, including bud growth capacity, relative conductivity, malonaldehyde contents, and soluble sugar contents. To understand the factors contributing to the difference in cold tolerance between ROC22 and GT42, comparative transcriptome analyses were performed. We found that genes involved in the regulation of the stability of the membrane system were the relative determinants of difference in cold tolerance. Additionally, genes related to protein kinase activity, starch metabolism, and calcium signal transduction were associated with cold tolerance. Finally, 25 candidate genes, including 23 variety-specific and 2 common genes, and 7 transcription factors were screened out for understanding the possible cold resistance mechanism. The findings of this study provide candidate gene resources for cold resistance and will improve our understanding of the regulation of cold tolerance at the bud stage in sugarcane.


Introduction
The frequency of extreme weather, including low temperatures, flooding, waterlogging, and drought (high temperatures), has recently increased. Among these abiotic stresses, cold stress (i.e., chilling and freezing) is particularly common, and the average minimum temperature of approximately 64% of the world total land is under 0 °C (Rihan et al. 2017). Cold stress mainly influences plant growth and development and reduces crop yields and/or the quantity of other useful materials (Pearce 2001). Some plants have evolved mechanisms to withstand cold stress. Under cold stress, physiological, biochemical, and cellular changes take place in plant cells. At the physiological and cellular levels, the osmotic potential, ice crystal formation, cell membrane stability, and concentration and distribution of reactive oxygen species (ROS) are altered substantially (Dong et al. 2009). Moreover, the soluble sugar concentration, cold resistance proteins, and proline content, which are important indexes  (Kaplan et al. 2007). A reduction in plasma membrane fluidity appears to be a primary event in cold signal transduction (Ding et al. 2019). In molecular aspect, the well-understood cold-signaling transduction pathway is the CBF-COR transcriptional cascade ). This pathway plays important roles in the promotion of several downstream cold-regulated (COR) genes (Thomashow 2010). Another crucial component of the CBF-COR cascade is the C-repeat binding factors (CBFs) gene family, which known as central nodes in the cold acclimation pathway and is conserved in many plant species. The basic helix-loop-helix transcriptional activator inducer of CBF expression 1 (ICE1) has been identified as a transcriptional activator. It can bind to MYC cis-elements within the CBF3 promoter, thereby repressing the expression of CBF3 and its target COR genes during cold acclimation. In addition, overexpression of ICE1 in Arabidopsis results in increased cold resistance, suggesting that ICE1 has a regulatory role in the cold stress response (Chinnusamy et al. 2003).
Sugarcane is an important crop with high economic and energy-related value (Li et al. 2017). In Southeast Asia, sugarcane is the main sugar crop (Zhang et al. 2011). In China, 60% of sugarcane is produced in Guangxi province. As a tropical crop, sugarcane is sowed in the spring (from late January to mid-march) in most production areas. However, freeze injury during the wintering stage can destroy the growth and development of sugarcane, thereby decreasing sugar production. For example, in 2008, freezing injury over a 23-day period (< 4 °C) was recorded in Guangxi province during the wintering stage, leading to major losses in some sugarcane varieties, especially ROC22. In January of 2018, a 5-day chilling period resulted in a reduction in sugar production in Guangxi province . During the last few decades, climate change, especially in southern China, the frequency of abnormal conditions is increasing. Chilling and freezing disasters are growing more frequent and unpredictable, resulting in high variation in production and economic losses. Thus, a comprehensive understanding of the mechanism underlying cold tolerance is necessary for the prompt development of solutions.
Several regulators of cold tolerance in sugarcane have been identified. For example, a NAC domain protein family member, SsNAC23, plays roles in the response to cold stress in sugarcane (Nogueira et al., 2005). The overexpression of the sugarcane alcohol dehydrogenase gene family member ScADH3 enhances the cold tolerance of transgenic Nicotiana benthamiana (Su et al. 2020). Furthermore, an uncharacterized gene encoding a coldregulated protein (Cor413) in Saccharum spontaneum is upregulated in leaf and root tissues under low temperatures (Dharshini et al. 2020). Generally, sugarcane is considered sensitive to low temperatures; however, field observations have shown that the degree of sensitivity varies among varieties (Nogueira et al. 2003a, b). However, these limited information is still insufficient for us to establish a systematic and complete recognition of cold tolerance mechanism in sugarcane. Thus, genetic resources related to cold tolerance can be detected by comparative analyses of different sugarcane varieties.
Although several genes involved in the regulation of cold responses in sugarcane have been identified at different growth stages, few studies have evaluated coldresponse genes during the bud stage, despite the important implications in the context of climate change. To address this gap in knowledge, we performed a comparative transcriptome assay using different sugarcane bud materials. Our results provide useful candidate genes for improving the cold resistance of sugarcane at the bud stage.

Plant Materials, Treatments, and Sample Collections
Sugarcane varieties ROC22 and GT42 were used as plant materials in this study. ROC22 is the main cultivar in Guangxi province, China. GT42 is introduced from Taiwan, China. The production and saccharinity in GT42 overmatch those in ROC22 (Deng et al. 2018). The culture method was mainly referenced the previous study (Yang et al. 2016;Qiu et al. 2019) with slight modifications. The individual stems of ROC22 and GT42 were cut to be onebud node and immersed with flowing water for 24 h to promote the sprout of buds and remove impurities, and then cultured in sand soil in the growth chamber (MLR-352H-PC, Panasonic, Japan) at 28 °C with the continuous illumination (30,000 lx) and 65% relative humidity. After the buds grown out and the length reached 4 cm, both of ROC22 and GT42 buds were transported into a lowtemperature (4 °C) growth chamber under the continuous illumination (30,000 lx), and the cold treatment maintained 48 h. The control plants were cultured under same condition except the culture temperature was 28 °C. The young buds were collected for the subsequent experiments. For the RNA sequencing sample collections, five buds of control group and cold treated group at 0 h, 24 h, and 48 h were collected using sterilized forceps and knifes, and frozen in liquid nitrogen immediately. Each of RNA-seq sample contained eight buds and was performed three biological replications. For the measurement of physiological indexes, five buds of ROC22 and GT42 at 0 h, 24 h, and 48 h post treatments were collected and stored in − 80 °C refrigerator or used for metrical experiments immediately.

Measurements of Morphological and Physiological Indexes
To monitor the growth status of ROC22 and GT42 under cold treatment, the lengths of buds were measured. The measurements of weights and lengths of buds were performed by measuring 20 buds of ROC22 and GT42 under cold treatment at 0 h, 24 h, and 48 h, respectively. To understand the physiological changes of ROC22 and GT42 under cold treatment, the relative conductivity (a physiological index that reflects cell membrane permeability under stress condition) and malondialdehyde (MDA) content were measured using the relative methods described previously (Li et al. 2008;Ranieri et al. 1993). For relative conductivity measurement, five buds of each time point were cleaned using distilled water, and then dry using absorbing paper. The buds were cut as 1 cm segments and immersed in 20 ml distilled water, and then were vacuumed for 8 min. The conductivity (recorded as S1) was measured using conduct meter at 25 °C. Finally, the bud tissues with distilled water were boiled for 15 min, and then measured the conductivity again (recorded as S2). The relative conductivity was the value of S1/S2. For MDA measurement, five buds of each time point were collected and ground in mortar with liquid nitrogen. Then, 1 g of tissue powder was suspended in 8 ml trichloroacetic acid (10%, v/v), and centrifuged for 10 min under 4000 rpm. The absorbency of supernate was measured using spectrophotometer. The calculation method was referenced as the description in the citing article. Meanwhile, the soluble sugar contents in bud tissues were also measured using anthrone method and the metrical method was described previously (Zou 1995), and five buds of each time point were used. All experiments were replicated three times.

Library Construction and Sequencing
Total RNA was extracted from the bud tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA concentration and purity were determined using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). RNA integrity was also checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). All samples had RNA integrity number (RIN) above 7.0. A total amount of 1.0 μg RNA per sample was used to prepare sequencing libraries using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, USA). An Illumina HiSeq 2000 platform was used for sequencing, and 150 bp paired-end reads were generated.

The Quality Control of Transcriptomic Data and Mapping on the Reference Genome
The total raw sequencing data were filtered by Phred quality score (Q ≥ 30) and read length (≥ 25 bp) using the Solex-aQA package. The reads containing adapters and poly-N, and low-quality reads were removed from the raw reads. The clean reads were mapped to the reference sugarcane genome (http:// sugar cane. zhang jisen lab. cn/ sgd/ downl oad/ Sspon. HiC_ chr_ asm. fasta. gz) (Zhang et al. 2018) using TopHat2 (Kim et al. 2013). The sugarcane genome assembly and gene annotation are deposited in the NCBI database under accession number QVOL00000000, BioProject number PRJNA483885 and BioSample number SAMN09753102 (Zhang et al. 2018).

Isolation and Functional Analyses of DEGs
The expression level of each gene were quantified based on the mapped reads counts using featureCounts software (v 2.0.1) with default parameters. For biologically duplicated samples, we performed an analysis using DESeq2 to screen differential expression genes (DEGs) (Anders and Huber 2010). The screening threshold were expression fold change ≥ 1 and a false discovery rate (FDR) < 0.01. The published sugarcane gene annotation was referenced to annotate all DEGs (Zhang et al. 2018). Gene Ontology (GO) enrichment analysis was performed using topGO (v 2.42.0) (corrected p-value < 0.05) (Alexa and Rahnenfuhrer 2020). The DEGs that were gained from control groups of ROC22 and GT42 were thought as developmental background and removed from DEG sets isolated from cold treated groups.

Gene Validation Using Quantitative Real-Time PCR
The First-strand cDNA synthesis was performed using M-MLV Reverse Transcriptase (Takara, Japan) according the protocol. The CFX96 Touch™ Real-time PCR detection system (Bio-Rad Laboratories, Hercules, CA, USA) was used for DEGs expression analysis using the SYBR Premix Ex Taq Kit (Takara, Japan). The thermal cycling program was used for PCR: 45 cycles of 95 °C for 20 s, 60 °C for 20 s, and 72 °C for 20 s after an initial denaturation step for 5 min at 95 °C. The 2 −ΔΔCT method was used to calculate gene expression. Each experiments performed three biological replications. The sugarcane glyceraldehyde-phosphate dehydrogenase (GAPDH) gene was used as the inner reference (Thiebaut et al. 2012), and all primers used in this study were listed in Table S1.

Statistical Calculation and Data Availability
SPSS 18.0 software was used to perform statistical calculation. One-way ANOVA with post hoc Tukey HSD (Honestly Significant Difference) test (p < 0.05) was carried out to detect the differences among different data sets. The RNA sequencing data sets were uploaded in Sequence Read Archive (SRA) database with the BioProject Accession Number PPJNA667204 (https:// www. ncbi. nlm. nih. gov/ sra/? term= PRJNA 667204).

ROC22 and GT42 Show Different Cold Sensitivity at the Bud Stage
Previous studies have shown that the sugarcane variety ROC22 is highly sensitive to low temperatures (Sun et al. 2012); however, little is known about this trait in GT42. To better understand the differences between ROC22 and GT42 under low-temperature treatment during the bud stage, we monitored the bud growth status at different time points. Firstly, we observed the morphological changes of buds of both ROC22 and GT42. Under normal condition, the growth of ROC22 and GT42 buds grew normally, and both of them started leaves at 48 h. After cold treatment, the growth of ROC22 was inhibited, while the growth of GT42 buds were similar as buds grown under normal condition (Fig. 1A).

Fig. 1
Changes of growth and physiological indexes of ROC22 and GT42 buds. A Morphological observations of ROC22 and GT42 buds under normal and cold conditions. Scale bars: 1 cm; B fresh weights of ROC22 and GT42 buds under normal and cold conditions. C Bud lengths of ROC22 and GT42 under normal and cold conditions; D relative conductivity of ROC22 and GT42 buds under normal and cold conditions; E MDA contents of ROC22 and GT42 buds under normal and cold conditions. F Soluble sugar contents in ROC22 and GT42 buds under normal and cold conditions. Data was based on 20 individual samples and One-way ANOVA with post hoc Tukey HSD (Honestly Significant Difference) test (**p < 0.05) was carried out to detect the differences among different data sets. The alphabets indicate differences among different data sets To quantify these growth differences, we next measured the fresh weights and lengths of buds. Comparing with the buds grown under normal condition, the fresh weights and lengths of ROS22 and GT42 buds were both decreased, indicating that the growth of ROC22 and GT42 was restrained under cold treatment (Fig. 1B, C). However, GT42 displayed higher relative growth rate due to GT42 showed obviously higher fresh weight and longer length after cold treatment (Fig. 1B, C).
Plant cell membrane permeability is an important determinant of plant stress resistance, including cold resistance (Liu et al. 2013). To evaluate the physiological changes of cell membrane in ROC22 and GT42 under low-temperature conditions, we evaluated relative conductivities and malonaldehyde (MDA) contents in bud samples at different time points. As shown in Fig. 1D, both ROC22 and GT42 displayed increases in relative conductivity after cold treatment, and ROC22 showed significantly higher relative conductivities than those in GT42. Under normal condition, both ROC22 and GT42 revealed low relative conductivity, and they displayed the similar increased rate of relative conductivity (Fig. 1D). Moreover, we found higher levels of MDA in ROC22 than in GT42 at each time point under cold treatment, while there were no clear differences between different time points under normal condition. (Fig. 1E). As a crucial cellular component and metabolite, soluble sugar also play key roles in plants to against cold stress. Thus, we next measured the soluble sugar contents in ROC22 and GT42 buds at each time point under normal and cold treatments. The soluble sugar contents in ROC22 and GT42 obviously increased after cold treatment but lower than that at each time point under normal condition. However, the increase in soluble sugar in GT42 was significantly higher than that in ROC22 at 24 and 48 h after cold treatment (Fig. 1F). Taken together, the results shown above indicated that ROC22 and GT42 maintained clear differences in the response to cold stress, and cold resistance was higher in GT42 than in ROC22 at the bud stage.

Summary of Transcriptome Data
The responses of ROC22 and GT42 to cold treatment differed significantly at 24 and 48 h, prompting us to evaluate the genetic factors underlying the response to cold treatment. To clarify the gene expression differences between ROC22 and GT42, we collected bud tissues at 0 h, 24 h, and 48 h after cold treatment and performed transcriptome sequencing. After filtering out adaptors and low-quality reads, transcriptome sequencing of 18 samples generated,  (Table S2). At least 94.62% of the clean reads showed eligible Q30 (%) values in all libraries, and the GC ratio of clean reads was higher than 54% (Table S2), indicating that our sequencing data were of sufficiently high quality for subsequent analyses. On average, 73.05% of clean reads were mapped to the sugarcane reference genome, and over 94% of mapped clean reads were uniquely mapped (Table S3).

Functional Analyses of Variety-Specific Differentially Expressed Genes (DEGs)
To investigate the effects of cold stress on ROC22 and GT42 from a genetic perspective, we isolated differentially expressed genes (DEGs) at 0 h, 24 h, and 48 h. In total, 703 DEGs, including 44 up-and 659 downregulated genes, were screened out in the 24 h vs. 0 h comparison in the cold-sensitive sugarcane variety ROC22 ( Fig. 2A), while 887 DEGs, including 318 up-and 569 downregulated genes, were identified in the 48 h vs. 24 h comparison ( Fig. 2A). In the cold-tolerant sugarcane GT42, 53 (35 up-and 19 downregulated) and 356 (10 up-and 346 downregulated) DEGs were identified in the 24 h vs. 0 h and 48 h vs. 24 h comparisons, respectively (Fig. 2B). Since continuous changes in physical signs were observed, we expected DEGs identified in both the 24 h vs. 0 h and 48 h vs. 24 h comparisons to play particularly important roles in the regulation of cold tolerance. Hence, we further isolated 66 and 13 DEGs in the intersection between the 24 h vs. 0 h and 48 h vs. 24 h comparisons in ROC22 and GT42 ( Fig. 2C and D). We found that 66 DEGs displayed three different expression patterns in ROC22 during cold treatment (Fig. 2E), while 13 DEGs in GT42 showed two different expression profiles (Fig. 2F).
To understand the functions of these variety-specific DEGs, a GO enrichment analysis was carried out to filter the genes involved in cold stress-related processes. We found that ROC22-specific DEGs were enriched for the terms starch catabolic process (4 DEGs) and membrane systemrelated terms, including Golgi membrane (1 DEG), membrane (1 DEG), plasma membrane (3 DEGs), vacuolar membrane (3 DEGs), cytoplasmic membrane-bounded vesicle (4 DEGs), and integral component of membrane (15 DEGs). In addition, 4 and 4 DGEs were related to calcium ion binding and beta-amylase activity, respectively (Fig. 3A). We also found that several GT42-specific DEGs were related to membrane system-related terms, including Golgi membrane (1 DEG), integral component of plasma membrane (1 DEG), cytoplasmic membrane-bounded vesicle (2 DEGs), and integral component of membrane (3 DEGs). Additionally, some GT42-specific DEGs were also involved in the stress-activated protein kinase signaling cascade (1 DEG) and activation of protein kinase activity (1 DEG). One DEG was associated with calcium ion transmembrane transport and one DEG was associated with calcium-transporting ATPase activity (Fig. 3B).

Functional Analyses of Common DEGs in Both ROC22 and GT42
We further evaluated common DEGs with predicted roles in the regulation of cold tolerance in both ROC22 and GT42. We compared variety-special DEGs and screened out six common DEGs (Fig. 4A). Each of these six common DEGs participated in several biological, cellular, and molecular processes. Notably, we found that Sspon.004C0011332 was related to the term 'Golgi membrane' and Sspon.004D0012690 was related to the terms 'integral component of membrane and cytoplasmic membrane-bounded vesicle' (Fig. 4B, Table 1). These results indicated that cell membrane system-related genes are collectively involved in the regulation of cold resistance in ROC22 and GT42, and the stability of the membrane system might explain the difference in cold resistance between the varieties.

Expression Profiles of Candidate Genes
In addition to broad analyses of similarities and differences in gene modules involved in the regulation of cold tolerance in ROC22 and GT42, detailed analyses of the spatiotemporal expression patterns of candidate genes can provide insight into gene function. Based on functional analyses, we identified 25 candidate genes, including 18 ROC22-specific, 5 GT42-specific, and 2 common DEGs (Table 1). To further investigate the expression patterns of these genes in ROC22 and GT42, we performed qRT-PCR assays. All 18 ROC22specific genes were upregulated in ROC22, consistent with expression trends determined by RNA-seq data. However, Venn diagrams of common DEGs in ROC22 and GT42, respectively; E, F Heatmap showed the expression patterns of common DEGs in ROC22 and GT42, respectively these genes did not display significant changes in expression in GT42 after cold treatment (Fig. 5). Inversely, five GT42-specific genes were upregulated under cold stress but were unchanged in ROC22 (Fig. 5). The two common genes were upregulated; however, levels of both genes were clearly higher in GT42 than in ROC22 (Fig. 5). These results demonstrated that these genes might be involved in cold tolerance in sugarcane via different expression patterns in ROC22 and GT42, and two common genes encoding a B-box zinc finger protein (Sspon.004C0011332) and CMP-sialic acid transporter (Sspon.004D0012690) might be important candidate markers for cold resistance in sugarcane.

Isolation of Transcription Factors that Responded to Cold Stress
Transcription factors (TFs) have important functions in the regulation of target genes at the transcriptional level. Several kinds of TFs, such as WRKY, NAC, and ethylene response factors (ERF), have been reported to play essential roles in the regulation of cold resistance. We next surveyed the number and type of TFs among candidate genes. Among ROC22-specific DEGs, five TFs, including two WRKY TFs (Sspon.003C0035890 and Sspon.007A0013201), two NAC TFs (Sspon.002B0025720 and Sspon.003B0009181), and one ERF (Sspon.004C0008870), were screened out ( Table 2). Two TFs were identified from GT42-specific DEGs, including one ERF (Sspon.005C0001440) and one C2H2 zinc finger TF (Sspon.001A0036160) ( Table 2). There were no TFs in the common DEG set. As determined by qRT-PCR analyses, all five ROC22-specific TFs were upregulated in ROC22 during cold treatment but were not affected in GT42 (Fig. 6). In contrast, GT42-specific ERF and C2H2 TFs were upregulated in GT42 under cold stress and showed no expressed differences in ROC22 (Fig. 6).

Validation of DEGs
To confirm the reliability of the RNA-seq results, 20 randomly selected DEGs in the 24 h vs. 0 h and 48 h vs. 24 h comparisons for ROC22 and GT42 were evaluated by qRT-PCR. As shown in Table 3, all DEGs showed consistent trends, apart from some differences in fold change values. These results indicated that our RNA-seq data were highly reliable.

Discussion
Despite an increasing global temperature, abnormal weather and cold extremes are still observed worldwide (Gupta and Deswal 2014;Horton et al. 2015). Especially in southern China, the frequency of low temperatures is increasing. This is a challenge from the perspective of tropical agriculture. To expand the growth region to high latitudes and altitudes, it is necessary to developed improved cold-tolerant crop varieties (Ding et al. 2019), including important economic and energy-related crops, such as sugarcane. However, our understanding of cold tolerance in sugarcane is still limited, especially at the bud stage. Therefore, uncovering the mechanism underlying cold tolerance and the identification

Cell Membrane Stability is Essential for Cold Tolerance in Sugarcane Buds
When plants are exposed to cold stress, the cell membrane is thought to be the primary site of injury due to its central role in the regulation of various cellular processes (Takahashi et al. 2013). Therefore, cell membrane stability is a key trait for analyses of abiotic stresses and the selection of tolerant genotypes in plants (ElBasyoni et al. 2017). Under stresses, membrane lipid peroxidation often occurs in the cell membrane, and MDA is a representative product of membrane lipid peroxidation. MDA contents can reflect the extent of injury in the cell membrane (Robert and Bewlery 1980). Membrane lipid peroxidation can also alter membrane conductivity. Under cold conditions, membrane conductivity is often increased in sugarcane (Zhu et al. 2018). In this study, we found that ROC22 was more sensitive to cold stimulus than GT42, as evidenced by the better growth of GT42 buds under cold stress (Fig. 1A). In addition, GT42 showed lower increases in relative conductivity and MDA contents than those of ROC22 (Fig. 1B, C). A previous study has reported that cell membrane stability is greater in plant varieties with high cold resistance (Zhu et al. 2018). Our results are consistent with this finding, indicating that cell membrane stability is essential for cold tolerance in sugarcane buds.
We also found similarities in the molecular mechanisms underlying cold resistance in ROC22 and GT42. Several DEGs, including phylloplanin, AAA-ATPase ASD, and sodium/hydrogen exchanger 2 encoding genes, were involved in processes related to the cell membrane in both  Table 1). Furthermore, we found that the calcium-related GO terms calcium ion binding, calcium ion transmembrane transport, and calciumtransporting ATPase activity are important in ROC22 and GT42 (Fig. 3), indicating that the calcium ion transmembrane transport capacity and calcium-dependent processes are altered under cold stress. These results provide insights at the genetic level and suggest that cell membrane stability is the core determinant of differences of cold resistance between ROC22 and GT42. Two common genes encoding a B-box zinc finger protein 24 (Sspon.004C0011332) and a CMP-sialic acid transporter 5 (Sspon.004D0012690) may play key roles in the regulation of cell membrane stability in sugarcane. The possible functions of these relative DEGs are in line with our observations of physical signs.
We also observed differences between the strains. Several variety-specific DEGs contributed to the response to cold stress. For example, we isolated a gene encoding mitogen-activated protein kinase kinase A (MAPKK; Sspon.003C0017630) in GT42 but not in ROC22 (Fig. 3, Table 1). Its expression levels are increased under cold stress in GT42 but are not affected by cold in ROC22 (Fig. 5). Given the essential roles of the MAPK cascade in several stresses (Chen et al. 2021), MAPK-mediated pathways may contribute to the increase in cold resistance in GT42, and these pathways are highly correlated with cell membrane stability. Taken together, our data reveal that cell membrane stability is an important indicator of cold resistance, and the candidate genes can be used to improving cold tolerance in sugarcane at the seeding stage.

Soluble Sugar and Starch Metabolism May Influence Cold Tolerance in Sugarcane Buds
Soluble sugars exert protective effects on plant cells against damage caused by cold stress by many mechanisms (e.g., by acting as osmoprotectants or nutrients as well as by interacting with the lipid bilayer) (Ma et al. 2009). Several soluble sugars, such as sucrose, glucose, fructose, ribose, and trehalose, often accumulate under low-temperature stress (Wei et al. 2019). These kinds of sugars function in free radical clearance and indirectly induce protein synthesis, thereby altering the cold tolerance of plants (Selvarajan et al. 2018).
In our study, we found that the soluble sugar contents in  (Fig. 1D). Consistent with this result, we found that several DEGs were enriched in beta-amylase activity in ROC22 (Fig. 3, Table 1), indicating that starch metabolism is a variety-specific process involved in regulating cold tolerance in sugarcane. As a storage polysaccharide, starch is found in sugarcane stalks, leaves and growing point regions (de Figueira et al. 2011). The accumulation of soluble sugars may be at least partially triggered by an increased rate of starch degradation under cold treatment (Sitnicka and Orzechowski 2014); thereby, starch metabolism is thought to be closely related to tolerance against multiple stresses (Krasensky and Jonak 2012). Four genes encoding betaamylase 3 were isolated, and Sspon.001A0018210 and Sspon.001B0023140 showed clear differences in expression levels under cold stress (Fig. 5), indicating that these two genes may be possible factors involved in regulating starch metabolism.

Transcription Factors May Play Roles in Regulating Cold Resistance in Sugarcane Buds
TFs are regulated by different signal transduction pathways and can directly or indirectly bind to cis-acting elements to alter the transcription efficiency of target genes related to plant growth, development, and stress responses. In this study, we found that some WRKY TFs were involved in the response to cold stress in ROC22 (Table 2). These two putative WRKY genes (Sspon.003C0035890 and Sspon.007A0013201) were upregulated in ROC22 but did not change in GT42 after cold treatment (Fig. 6). WRKY TFs, which are among the well-characterized class of TFs in plants, have reported roles in cold tolerance, and 39 WRKY TFs have been isolated in sugarcane (Javed et al. 2020). Proteins in the NAC TF family are specific to plants and play important roles in plant growth, development, and senescence (Kikuchi et al. 2000). In a pioneering study, three putative NAC TFs were identified as cold-response genes in sugarcane (Nogueira et al. 2003a, b). Another NAC  6 Expression patterns of transcription factors (TFs). qRT-PCR assay was performed to detect the relative expression levels of each TF. The expression level at 0 h in ROC22 was set as 1. Each of experiment was replicated three times. Bars shown means ± SD TF, SsNAC23, is induced by low temperatures in sugarcane (Ditt et al. 2011). Here, we found that the expression levels of two NAC TFs increased by cold treatment in ROC22 but remained unchanged in GT42 (Table 2, Fig. 6). ERFs are a plant-specific TF family that was initially identified in tobacco as binding proteins for reduced sensitivity to disease (Wang et al. 2014). Increasing studies have revealed that ERFs play key roles in plant stress responses (Wang et al. 2014). SodERF3, a sugarcane ERF, enhances drought and osmotic stress tolerance in transgenic tobacco (Trujillo et al. 2008). However, the functions of ERFs in cold resistance in sugarcane are largely unknown. In this study, we isolated one ERF from ROC22 and one ERF from GT42 (Table 2) and observed distinct expression profiles in the two varieties (Fig. 6), suggesting that these two ERFs function in a variety-specific manner. More important, cold stress may result in water shortage, therefor give rise to the expression differences of these ERFs. Meanwhile, we also screened out a C2H2 zinc finger TF from GT42 ( Table 2) and found that it was a cold-response protein in GT42 but insensitive to cold in ROC22 (Fig. 6). So far, there are very limited information of the biological roles of C2H2 in regulating cold tolerance in sugarcane. Thus, our data provide novel insights into the molecular mechanisms underlying cold tolerance in sugarcane during bud stage.

Conclusion
In this study, we confirmed that the sugarcane variety ROC22 is more sensitive than GT42 by monitoring physical signs under cold treatment during the bud stage. Under cold stress, GT42 buds show higher fresh weights and lengths. The results of relative conductivity and MDA contents reveal that GT42 buds have stronger stabilities of cell membrane under cold treatment. Through performing a comparative transcriptome analysis using the bud tissues, we identified candidate genes involved in the regulation of cold acclimation at the bud stage. The candidate genes, which involved in cell membrane stability, soluble sugar metabolism, and starch metabolism, support us new insights into cell membrane-based regulatory mechanism of cold tolerance in sugarcane. Two common genes, which encode B-box zinc finger protein 24 and CMP-sialic acid transporter 5, may be potential candidates involved in cold resistance regulation. More importantly, TFs isolated in this study, including WRKY, NAC, ERF, and C2H2, may regulate these two gene clusters at the transcriptional level. Interactions among these candidate genes may have synergistic effects on cold resistance during sugarcane bud stage, and the possible functions of these candidates should be studied from genetic aspect in future. Based on these findings, we speculate that the TFs including WRKY, NAC, ERF, and C2H2 may synchronously regulate the cell membrane stability and soluble sugar and starch metabolisms to response cold stress, in turn to control cold resistance in sugarcane from transcriptional control aspect during bud stage, and more verified works need to be done in the future (Fig. 7). Overall, our results provide useful information on differentially expressed gene modules and metabolic and signaling pathways related to the response to cold stress in sugarcane during bud stage. Our data also provide a genetic resource for studies of cold tolerance regulatory genes and for sugarcane breeding.
need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.