The bermudagrass (Cynodon dactylon) materials were collected from 16 different wild sites in different regions of China (Supplementary Table 1). All genotypes were grown in the experimental field of Yangzhou University. In order to eliminate the effect of cytotype on gene expression, three replicated tetraploids with similar growth status were respectively, selected from different latitudes in this study, based on our previous chromosomal observations of the bermudagrass root tip. The 48 samples from different latitudes were transplanted into 48 plastic pots (20 cm depth × 16 cm diameter). The pots were filled with a mixture of sand and organic soil (1:1, v/v). Plants were cultivated in a greenhouse with a photoperiod of 12 h, a temperature of 25°C, and a humidity of 60%. After growth for 1 month, all the samples were equally divided into three groups. Each group was respectively, transferred into their respective phytotron under standardized growth conditions (12 h photoperiod, 20°C). After 7 days, the samples were subjected to 35°C (T35), 20°C (T20), and 5°C (T5) for 24 h. The fully expanded leaves were harvested immediately after the completion of treatment. All of the leaves were frozen in liquid nitrogen and stored at −80°C for RNA extraction.

RNA-sequencing was performed at Beijing Biomarker Technologies Co. Ltd. (Beijing, China). In brief, total RNA was isolated from mature leaf tissues using TRIzol according to the manufacturer’s protocol. A total of 137 cDNA libraries were constructed as previously described and sequenced by Illumina NovaSeq. After removing the adaptor sequences and low-quality reads, trimmed reads were de novo assembled into reference contigs with Trinity software using default settings (Grabherr et al., 2011). Finally, the contigs were clustered into unigenes using the TIGR Gene Indices clustering (TGICL) software system (Pertea et al., 2003). Raw data were submitted to the National Center for Biotechnology Information (NCBI) under the BioProject PRJNA646313 and PRJNA649353.

The transcriptomic resources of C. dactylon were submitted as BLASTx queries to functionally annotate expressed genes according to the NCBI non-redundant (Nr) protein sequences, clusters of orthologous groups (COG), gene ontology (GO) database, euKaryotic orthologous groups (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) public databases (Altschul et al., 1997). After predicting the amino acid sequences of unigenes, HMMER v3.1 (hidden Markov models) software was used to compare sequences to the Protein family (Pfam) database to obtain annotation information of unigenes (Finn et al., 2016).

The fragments per kilobase of transcript per million mapped reads (FPKM) were used to represent the expression abundance of each assembled unigene. Differential expression analysis of two samples was performed using edgeR software (Robinson et al., 2010). Fold change (FC ≥ 2) and a false discovery rate (FDR) < 0.01 (P-value < 0.05) were used to analyze the significant differentially expressed genes (DEGs) by pairwise comparisons to compare three treatments.

The plant heights (HTs) and root lengths (RLs) were determined by direct measurement and five replications were set in each pot.

Electrical conductivity (EL) assay: a fragment of 0.1 g leaves was soaked in water for 24 h at room temperature, and the initial conductivity (Ci) and final conductivity (C_max) after autoclaving at 121°C for 30 min were measured with a conductance meter (DDS-11A, Aolilong, China). The relative EL% = (Ci)/C_max × 100% (Hu et al., 2016).

The photosynthetic index, including net photosynthetic rate (Pn), stomatal conductance (Gs), and transpiration rate (Ts), were measured by inserting 4–5 individual leaves (second fully expanded leaves from the top) with a portable photosynthesis system (Li-6400, Li-COR Inc., Lincoln, NE, United States). The measurement conditions were: flow rate set at 500 μmol s^–1 and a fluorescent light intensity of 1,000 μmol m^–2 s^–1 (Chater et al., 2018).

Antioxidant assay and malondialdehyde (MDA) content: 0.2 g fresh leaves were ground in a mortar and a crude enzyme solution was extracted with 4 mL of phosphate buffered saline (PBS) buffer. The crude enzyme was applied for the determination of MDA content. After the thiobarbituric acid (TBA) reaction, the absorbance of the supernatant was examined at 532 and 600 nm, and the MDA content was calculated as described in a previous study (Dhindsa, 1981). The peroxidas (POD) mixture contained 50 μL of crude enzyme solution, 1.85 mL sodium acetate-acetic acid buffer (0.1 M, pH 5.0), 0.25 mL 0.25% guaiacol (dissolved in 50% ethanol solution), and 0.1 mL 0.75% (v/v) H₂O₂. The reaction system of superoxide dismutase (SOD) mixture contained 1 mL of crude enzyme solution, 2.2 mL sodium phosphate buffer (50 mM, pH 7.8), 0.039 mM methionine, 0.3 nM ethylene diamine tetraacetic acid (EDTA), 0.012 μM riboflavin, and 0.225 μM nitro blue tetrazolium (NBT). The activities of POD and SOD were measured with a spectrophotometer (UH5300, Hitachi, Japan) based on Hu et al. (2012).

Statistical analyses were performed using SPSS21.0 software. For correlation analysis, the Pearson correlation coefficient (r) was calculated and a two-tailed test was carried out. For weight evaluation analysis, the technique for order preference by similarity to ideal solution (TOPSIS) method was used in R software.

We used weighted gene co-expression network analysis (WGCNA) to build plant traits and bermudagrass gene modules (clusters of genes displaying similarly correlated patterns of transcription) (Langfelder and Horvath, 2008). As it is interesting to understand the relationship between gene expression and plant traits of different bermudagrass germplasms treated at different temperatures, the DEGs of the T5/T20 group and the T35/T20 group of one site were gathered to form a new group using R software, and so on. Finally, all DEGs of the 16 sample sites were gathered to form a differential expression pool. We used the pool of DEGs and quantitative traits to build bermudagrass gene modules. The default WGCNA “step-by-step network construction” analysis was used to build the modules. We calculated the adjacency relationship between genes and constructed a topological matrix. A hierarchical clustering tree with the dissimilarity of the topological matrix was generated, and the modules of dynamic trees were cut by a standard method. Then, similar modules were merged by calculating the eigengenes module. Finally, a cluster dendrogram was formed with a soft-power of 6, a minimum module size of 50 genes, and a distance threshold cut of 0.1 (Supplementary Figure 1).

Additionally, the modules-trait membership was determined by measuring the Pearson correlation between bermudagrass module eigengenes and plant quantitative traits using the default WGCNA relating modules to external information analysis. Then, the correlation and connectivity for each trait to the corresponding correlated module were calculated and shown on the modules-trait membership. According to the significance of the correlation between the trait and gene expression of each module, the mean value of each module represents the significance of the trait within the module, so the module with the largest absolute significance value has the greatest correlation with the given trait. Based on this information, a data set of “hub genes” in the module with the greatest significance of each trait was selected, which was most correlated to the modules of the bermudagrass trait. Finally, the sequences and related information of the “hub genes” dataset were extracted from the transcriptome sequencing data, and then we analyzed and interpreted these genes by Mapman software (Usadel et al., 2009).

Total RNAs were extracted as described above for the samples used for RNA-seq. The quantitative real-time (qRT)-PCR analysis was performed to validate the reliability of the RNA-seq data. Four DEGs were chosen for qRT-PCR analysis. Each sample at different latitudes had three replicates. The qRT-PCR assays were performed using the UltraSYBR Mixture (Kwbiotech, Beijing, China) and were conducted on a Roche LightCycler 96 Sequence Detection System. The reference gene, CdActin, was used to normalize the expression levels of target genes. All primers are shown in Supplementary Table 2. The expression levels were calculated using the 2^–ΔΔCt method (Livak and Schmittgen, 2001).

In total, 20–42 million clean reads per library were obtained from wild bermudagrass leaves after removing adapter-related and low-quality reads. The GC content ranged from 41.65 to 57.51%, and the respective Q30 percentage (sequencing error rates < 0.1%) was 89.84–94.91%, respectively. After assembly, a total of 553,110 unigenes were obtained with an N50 of 663 bp and a mean length of 510 bp, respectively. Sequence alignment was performed with the clean data of each sample and the assembled unigene library, and mapped reads were used for further analysis. A total of 313,701 unigenes were annotated against the NCBI Nr protein database, Swiss-Prot, KOG, COG, GO databases, and KEGG, respectively.

As shown in Table 1 and Figure 1, the DEGs in the T5/T20 group of most samples along latitudinal gradients are higher than those in the T5/T35 and T35/T20 groups. The results showed that compared with higher temperatures, most genes were expressed following low-temperature treatment. Second, compared to the samples at low and high latitudes, the number of DEGs in the samples at mid-latitude L7 and L8 of each temperature comparison group was less, which suggested that the wild bermudagrass of the two sites may be not sensitive to low or/and high temperature. In addition, among all the same latitudes, the number of DEGs shared by the three temperature comparison groups is less, which indicates that these wild bermudagrasses have altered responses to different temperatures at the transcription level (Figure 1).

In order to better compare the responses of bermudagrass at low and high temperatures, the expression patterns in the T20 group were assigned as the reference samples. A total of 14,654 DEGs from the T5/T20 and T35/T20 comparisons were identified and combined across all latitudes. To understand the function of DEGs in two comparative treatments, these DEGs were searched against the KEGG database using the TBtools software (Chen et al., 2020). Among these pathways, the DEGs of the T5/T20 group were enriched in functions about the ubiquitin system, protein kinases, and glycosyltransferases. Unlike the low-temperature group, chaperone catalysts and metabolism of other amino acids were only found in the T35/T20 comparison (Supplementary Figure 1).

The WGCNA was used to analyze the association between the 14,654 DEGs identified from T5/T20 and T35/T20 comparisons and physiological indicators. A total of 39 modules related to physiological indicators were observed (Figure 2). Analysis of the module-trait relationships revealed that the “green” module of 256 genes was mostly correlated with SOD (r = −0.86, P = 1 × 10^–14), POD activities (r = −0.86, P = 3 × 10^–14), and EL (r = −0.5, P = 3 × 10^–4). Therefore, the genes in the “green” module may be related to the antioxidant processes of bermudagrass in response to cold/heat stress. The “red” module of 212 genes, “dark turquoise” of 34 genes, and “magenta” of 222 genes appeared to be associated with HT, RL, and MDA content. The genes related to these traits are closely related to the plant growth of bermudagrass under cold/heat stress conditions.

First of all, based on the MapMan functional annotation, 724 genes related to different physiological phenotypes in the above modules were identified. As shown in Figure 3, these genes were enriched in the regulation of transcription, protein metabolism, photosystem, minor carbohydrates, and hormone metabolism.

The differentially expressed gene expression profiles in plant growth and metabolism of bermudagrass under different temperature treatment were covered in various pathways, including regulation of transcription, hormone and terpenoid metabolism, tetrapyrrole synthesis, and other processes (Figure 4A). The pathway of transcriptional regulation comprises some TFs, including one NAC (NAM, ATAF1/ATAF2, and CUC2 homologous proteins) domain-containing protein 1 (NAC1), five basic helix-loop-helix DNA-binding superfamily proteins (bHLH), and one pseudo-response regulator 37 (PRR37). The homologous genes PRR37 and PRR7 were reported to be related to the components of a temperature-sensitive circadian system. The gene expression of PRR37 was also upregulated in T20 and T35 compared with the T5 group. In addition, some TFs, such as the TCP (Teosinte Branched/Cycloidea/Proliferating cell factors) and DUF (domain of unknown function) family TFs, also participate in the metabolic pathway of plant hormones (Figure 4A).

The primary and secondary metabolites of terpene metabolism, such as phytosterols, chlorophyll, and phytohormones, are involved in regulating plant growth and development. The hormone pathway covers 10 genes, which encode proteins related to biosynthesis and signal transduction of auxin, cytokinin, and ethylene, such as auxin signaling F-box 2 (AFB2), auxin-responsive GH3-like protein 5 (GH3.5), histidine kinase 3 (HK3) of cytokinin receptor, and adenine nucleotide alpha hydrolases-like superfamily protein (ANAH) of ethylene signal transduction. In terms of terpene and chlorophyll biosynthesis, except for CHLG encoding chlorophyll synthase, seven genes were covered and upregulated in the T35 treatment group (Figure 4A). Chlorophyll synthase has been shown to perform the esterification of chlorophyllide (a and b), the last step of chlorophyll biosynthesis. Compared to T35 treatment groups, the CHLG gene was upregulated in the T5 treatment, while the gene encoding magnesium-protoporphyrin IX monomethyl ester cyclase (MPEC) was downregulated in the T5 groups.

Generally, HT is a characteristic of turfgrass growth and a measure of the development and resistance under abiotic stress. Through our analysis of bermudagrass HT at different latitudes after different temperature treatments, compared with the T20 treatment, it is clear that T5 treatment was significantly lower without considering the latitude, while the T35 treatment had no significant difference (Figure 4B). These results were consistent with the expression of the plant growth-related genes as mentioned above.

In order to explore whether latitude affects HT of wild bermudagrass from different regions, the correlation analysis between latitude and HT was performed on the samples of the T20 treatment group (Figure 4C). First of all, the coefficient between latitude and height of bermudagrass was only −0.07 in all regions. Then, according to the latitude data (Supplementary Table 1), the sampling position could be divided into subtropical (L1–L10) and mesothermal zones (L11–L16). So, an exciting discovery is that the height of subtropical bermudagrass was negatively correlated with latitude (coefficient is −0.9120), while the mesotherm zone was positively correlated with latitude (correlation coefficient is 0.8004). Thus, the HT of wild bermudagrass was related to latitude and climate, which also provides a reference range for exploring wild germplasm with cold resistance in these materials.

As mentioned above, the related genes of chlorophyll synthesis in bermudagrass exhibited decreased expression to some extent after temperature stress. We next studied whether the related photosynthetic system was similarly affected. Taking the DEGs in WGCNA modules as experimental objects, the transcriptional expression pattern of genes related to photosynthesis was explored. Analysis of 32 mapped DEGs showed that the most important genes involved in photosynthesis were downregulated under low-temperature treatment, including the genes encoding the subunit of photosystem I (PSA), photosystem I/II light-harvesting complex (LHCA, LHCB), chlorophyll a/b-binding protein (CAB), etc. (Figures 5A,B). Furthermore, the net photosynthetic rate (Pn), stomatal conductance (Gs), and transpiration rate (Tr) of plants with different temperature treatments were determined. The Figures 5C–E showed the effect of the different temperature treatments on the estimated photosynthetic indexes (Pn, Gs, and Tr). The Pn of bermudagrass in the T5 treatment was lower than that of the T35 and T20 groups. Overall, these results indicate that the photosynthesis was blocked at low temperature, while the CHLG gene (the enzyme encoding gene of the last step of chlorophyll synthesis) was upregulated in T5 treatment (Figures 4A, 5B).

To shed light on the tolerance mechanism of bermudagrass under temperature stress, 11 resistance-related genes obtained from WGCNA association modules were annotated and analyzed. In addition, the MDA content, EL, and antioxidant enzyme activity (POD, SOD, and ascorbate peroxidase (APX)) of the plants with different treatments were also determined. First, Figures 6A,B shows that MDA content and EL of the T5 group were higher than those of the T20 and T35 groups, which indicated that bermudagrass under 5°C treatment suffered damage, but there was no significant difference between the T20 and T35 treatment groups. For the antioxidant system, the activities of POD, SOD, and APX of bermudagrass at most latitudes in T5 treatment were higher than that in the T20 (Figures 6C,D and Supplementary Figure 2). Similarly, the expression levels of some genes encoding POD and SOD were higher in the T5 groups, such as the peroxidase 3 gene (Supplementary Figure 3). With the decrease of temperature, more antioxidant enzymes accumulated, thus delaying the damage caused by temperature stress. In addition, the temperature treatment time in this study was 24 h, which indicated that these wild bermudagrasses responded quickly to low-temperature stress, which laid a foundation for screening cold-resistant germplasm.

Furthermore, plant HSPs and Hsfs are multifunctional, as they regulate protein folding in processes ranging from abiotic stress to plant development. As shown in Figure 6E, the transcript abundance of HSP90.1, HSP70-1, and heat shock cognate protein 70-1 (HSC70-1) were higher in the T5 and T20 treatments compared with T35 groups. On the contrary, heat-stress-associated 32-kD (HSA32), two chaperone DnaJ-domain superfamily protein genes, and a member of Bcl-2-associated athanogene (BAG) family involved in heat stress were upregulated in the T35 treatment.

As a critical TF for regulating abiotic stress, 28 HSF transcripts were found from the pool of bermudagrass DEGs and their FPKM values were analyzed. As shown in Figure 7, compared to the expression level of the T20 treatment, HsfA3, A8a, and C1b were significantly upregulated at 5°C, while HsfA2d, A2c, and its homologs were upregulated in the T35 treatment. Besides, HsfB2c and B1a were both downregulated at 5 and 35°C, indicating that the expression of both genes was sensitive to temperature. Taken together, these results indicate that a complex regulatory network of the HSP70, HSP90, HsfA3, and HsfA8a genes could be more involved in response to low-temperature stress.

To explore the cold-resistant germplasms from wild bermudagrass resources collected from 16 regions at different latitudes, the differences of the related physiological indexes between the samples treated at 5 and 20°C were determined, and the weight evaluation of variation was evaluated (we provided Supplementary Table 4 for the detailed data and weight ratio). The physical parameters, including Pn, Gs, Tr, MDA content, EL, POD, SOD, and APX, were integrated by the TOPSIS method for comprehensive resistant evaluation. By calculating the distance ratio between the positive and negative ideal samples, the bermudagrass germplasm resources collected from 16 different areas were sorted. As shown in Table 2, the bermudagrass of the L12 might have strong cold resistance.

In order to further verify the reliability of our transcriptomic data, four DEGs were selected for qRT-PCR analysis, including HSP90.1, HsfA3, and HsfA8, which participate in the response of temperature stress and LHCA3 involved in the regulation of plant photosynthesis. As shown in Figure 8, the FCs determined by the expression qRT-PCR is consistent with the difference of FPKM in most samples, and their expression patterns were similar. For example, HSP90.1 was downregulated in the T35 treatment, while HsfA3 and HsfA8 were significantly upregulated in the low-temperature group. Taken together, our qRT-PCR analysis indicated that the sequencing results were reliable and provided a basis for subsequent analyses.

Although the gene expression and physiological performance associated with cold stress in turfgrass have been extensively studied in several wild and commercial species (Shi et al., 2015; Wei et al., 2015; Bi et al., 2016), combined studies of global gene expression, plant growth, photosynthesis, and enzymatic activities in wild bermudagrass from a latitudinal gradient are scant.

Overall, 14,654 DEGs were identified from the combination of T5/T20 and T35/T20 comparisons in all samples from all latitudes, which accounted for the transcriptional regulation of bermudagrass subjected to temperature stress. The 14,654 DEGs enabled us to identify six WGCNA modules of 724 genes, respectively, correlated with HT, RL, MDA content, EL, POD, and SOD activity (Figure 2). WGCNA has been used widely for similar analyses in other studies (El-Sharkawy et al., 2015; Farcuh et al., 2017). Therefore, the identification of some stress-associated “hub genes” was important and meaningful in this study. According to annotation, these hub genes were linked to plant growth and stress response, involved in terpene metabolism, tetrapyrrole synthesis, photosystem, regulation of TFs, and other pathways (Figures 4, 5A, 6B, 7). First, the low temperature has decreased the HT of bermudagrass (Figure 4B), and there are various and complex regulatory networks of related gene expression, such as well-known terpenoid synthesis and hormone metabolism and signal transduction. For the photosynthesis system, it is speculated that cold stress has a negative effect on photosynthesis of bermudagrass from the term of decreased Pn and lower expression of LHCA gene and others in electron transport chain under low-temperature treatment. However, the internal relationship between Pn and the efficiency of electron transport chain is worth further exploration. While, for the terms of chlorophyll, the CHLG gene expression was upregulated in the low-temperature treatment, which may be due to the adaptation strategy of bermudagrass under cold stress or the negative feedback regulation of the high expression of MPEC gene in the chlorophyll biosynthesis pathway.

It is known that HSPs are the major induced proteins in response to stress, including cold and oxidative stress (Dumas et al., 2019). The inducible HSP70s (70-kDa HSP), as well as the constitutive HSC70, act as molecular chaperones, play crucial roles in protecting denatured proteins and prevent substantial damage to plants exposed to cold stress (Franzellitti and Fabbri, 2005). For example, the high expression of HSP70 transcripts had a strong relationship with resistance to cold stress in Leguminivora glycinivorella (Wang et al., 2014), and the mRNA levels of HSP70s (including HSP70A, B, C) were increased in Chlamydomonas reinhardtii cells exposed to cold stress (Maikova et al., 2016). Similarly, in this study, the expression levels of HSP70 and Hsc70-1 genes were upregulated in both 5 and 20°C grown bermudagrass compared with the 35°C treatment (Figure 6E). Although the FC of its expression in the T5 treatment group was lower than that in the T20 treatment group, there was no significant difference. Consistently, there was no significant change in HSP70 expression of cucumber after 1 day of 4°C treatment compared to a control of 20°C treatment as Ru et al. (2020) reported. These results support that the regulation of HSP70 and HSC70-1 is most likely a strategy that bermudagrass uses to react to cold stress. Then, HSP90 is one of the critical proteins to decrease cellular damage via cross talk with other mechanisms in the face of environmental stress (Wang et al., 2004; Hossain et al., 2012). Overexpression of GmHSP90A reduced chlorophyll loss and lipid peroxidation levels, so HSP90A likely participates in decreasing oxidative stress injury under heat stress (Xu et al., 2013). In the present study, the transcript level of HSP90.1 and CHLG was higher in the T5 and T20 treatments than those in the T35 treatment, and there was no obvious difference in MDA content between the T5/T20 and T35 treatments (Figures 4A, 6A,E).

Genes encoding HSP70 and HSP90A, which are the components of HSP chaperone machine, are both upregulated in Chlamydomonas reinhardtii response to cold stress (Maikova et al., 2016). Likewise, HSP70 and HSP90 transcripts were also significantly increased in Ageratina adenophora after cold treatment, which together could play a fundamental role in protecting cells from stress damage (Gong et al., 2010). In this study, the upregulation of the HSP70/90.1 gene under low temperature is closely related to the regulation mechanism of bermudagrass adapting to cold stress. And then, the higher expression of the HSPs gene in T5 treatment could promote the CHLG gene expression and inhibit the MDA content in T5 treatment. Although the net photosynthetic rate is still decreased in T5 treatment, its value was close to that of the T20 groups, which may be due to the mechanism of bermudagrass adapting to cold stress, including HSPs regulation. It is suggested that HSPs could direct the protein folding, as well proper assembly of the protein complex such as Rubisco (Park and Seo, 2015). For example, a tomato chloroplast-targeted DnaJ protein could target Rubisco to maintain photosynthesis under temperature stress (Wang et al., 2015). Overall, the accumulation of HSP70 and HSP90 genes in bermudagrass to a low temperature of 5°C may reduce the damage under stress by enhancing chlorophyll synthesis and inhibiting peroxidation.

RNA-sequencing data analysis detected 28 differentially expressed Hsfs among the 14,654 DEGs, but none of them were members of the WGCNA highly correlated modules (Figure 7). However, Hsfs, as the large gene family, are highly redundant and participate in a flexible gene network that regulates abiotic stress responses (Miller and Mittler, 2006). Gad and Ron summarized the findings from a considerable number of studies and described that Hsfs function as hydrogen peroxide sensors in mammals and yeast, and then presented the hypothesis that HsfA4a and HsfA8 might function as potential H₂O₂ sensors (Miller and Mittler, 2006). Furthermore, Song et al. (2016) also reported that the expression of HsfA3 was induced by oxidative damage in Arabidopsis. Although Hsf genes in this study are not the “hub genes” of WGCNA, they showed differential expression under different treatments (Figure 7), which indicated that Hsfs might not be directly responsible for the response to cold/heat stress and coincide with flexible and redundant roles. ROS are generally considered as the early signals of various stresses in plants (Mittler, 2004; Karuppanapandian et al., 2011). As mentioned earlier, HsfA3 and HsfA8a of bermudagrass were only upregulated at 5°C (except for HsfA3 in T35/T20_13 samples, Figure 7). Could HsfA3 and/or HsfA8a act as hydrogen peroxide sensors in bermudagrass? Besides, HsfB2c and B1a were responsive to both 5 and 35°C conditions compared with 20°C, indicating that both of these two genes are sensitive to temperature. Of course, it is necessary to determine whether Hsfs such as HsfA3, HsfA8a, HsfB1a, and HsfB2c function as redox sensors in future research.

Under heat/cold stress, hydrogen peroxide can be eliminated by the antioxidative enzymes system, including POD and SOD (Yang et al., 2019). In the current study, the POD, SOD, and APX activities were significantly induced by 5°C stress, and some encoding peroxidase genes were upregulated in T5 and T20 treatments (Figure 6 and Supplementary Figures 1, 2). This begs the question, is there a direct relationship between the levels of antioxidant enzymes, Hsfs, and HSP70 induced by low temperature? Previous studies have shown that regulation of HSP70 and SOD-1 were positively correlated, and some induced Hsfs directly regulated antioxidant enzymes in Arabidopsis, such as APX (Shin et al., 2013; Song et al., 2016). As well as others, our results demonstrate that POD, SOD, and APX activities and the expression level of HsfA3/A8a, Hsp70/Hsp90 coding genes are both increased in 5°C treated bermudagrass at all latitudes.

Based on differential analysis of physiological performance and gene expression related to plant growth and response to cold stress between 5 and 20°C treated bermudagrass, the above discussion about the regulatory network of heat shock modules-photosynthesis-redox system in bermudagrass is illustrated in a schematic display (Figure 9). Although our results are preliminary, we speculate that there is a direct or indirect interaction among HsfAs, antioxidant enzyme activities, and HSPs to coordinate a systematic response to cold stress.

Finally, based on the weight evaluation analysis of a dataset of the physiological indexes (Pn, Gs, Tr, EL, and MDA content as well as antioxidase activities), there may be an assumption that the germplasm located in Zhumadian (L12) of mesothermal zones with particular climate had good performance under cold stress. Of course, whether this accession could be cultivated as a cold-resistant bermudagrass variety needs further study. At the same time, through comparative transcriptomic analysis of different samples under different temperature treatments, we found that a set of genes encoding HSPs and Hsfs were not only regulated by temperature stress, but also directly or indirectly interplayed with peroxide scavenging and chlorophyll synthesis under cold stress. Of course, the specific functions and regulation mechanisms of the HSP70/90 and HsfA3/A8 genes respond to cold stress. This work provides not only essential clues to further clarify how bermudagrass responds to temperature stress, but also provides genetic resources for the cultivation of excellent bermudagrass germplasm.