Next Article in Journal
Insights into Metabolic Reactions of Semi-Dwarf, Barley Brassinosteroid Mutants to Drought
Previous Article in Journal
The COVID-19 Pandemic: Does Our Early Life Environment, Life Trajectory and Socioeconomic Status Determine Disease Susceptibility and Severity?
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptomic Profiling of Young Cotyledons Response to Chilling Stress in Two Contrasting Cotton (Gossypium hirsutum L.) Genotypes at the Seedling Stage

1
College of Agronomy, Northwest Agriculture and Forestry University, Yangling 712100, China
2
State Key Laboratory of Cotton Biology, Institute of Cotton Research, Chinese Academy of Agricultural Sciences, Anyang 455000, China
3
College of Agronomy, Hebei Agricultural University, Baoding 071001, China
*
Authors to whom correspondence should be addressed.
Int. J. Mol. Sci. 2020, 21(14), 5095; https://doi.org/10.3390/ijms21145095
Submission received: 17 June 2020 / Revised: 14 July 2020 / Accepted: 17 July 2020 / Published: 19 July 2020
(This article belongs to the Section Molecular Biology)

Abstract

:
Young cotyledons of cotton seedlings are most susceptible to chilling stress. To gain insight into the potential mechanism of cold tolerance of young cotton cotyledons, we conducted physiological and comparative transcriptome analysis of two varieties with contrasting phenotypes. The evaluation of chilling injury of young cotyledons among 74 cotton varieties revealed that H559 was the most tolerant and YM21 was the most sensitive. The physiological analysis found that the ROS scavenging ability was lower, and cell membrane damage was more severe in the cotyledons of YM21 than that of H559 under chilling stress. RNA-seq analysis identified a total of 44,998 expressed genes and 19,982 differentially expressed genes (DEGs) in young cotyledons of the two varieties under chilling stress. Weighted gene coexpression network analysis (WGCNA) of all DEGs revealed four significant modules with close correlation with specific samples. The GO-term enrichment analysis found that lots of genes in H559-specific modules were involved in plant resistance to abiotic stress. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that pathways such as plant hormone signal transduction, MAPK signaling, and plant–pathogen interaction were related to chilling stress response. A total of 574 transcription factors and 936 hub genes in these modules were identified. Twenty hub genes were selected for qRT-PCR verification, revealing the reliability and accuracy of transcriptome data. These findings will lay a foundation for future research on the molecular mechanism of cold tolerance in cotyledons of cotton.

1. Introduction

Temperature is an essential environmental factor limiting the geographical distribution and growing season of plants [1]. Cold stress is classified as chilling (0–20 °C) and freezing (<0 °C), which are regulated by different mechanisms [2]. Low temperature is detrimental to plant growth and survival, as well as the yield and quality of thermophilic crops. As a commercial crop originating in low-latitude regions, the optimum growth temperature of cotton is 20–30 °C [3]. However, different varieties of the same plant species are generally with different tolerance to cold stress [4,5]. Tropical plants have been reported to be sensitive to cold stress, especially during the transition from heterotrophic to the autotrophic stage, and have different tolerance to cold stress at different growth stages [6,7]. Therefore, identification of tolerant and sensitive varieties at a specific growth stage is of great help to reveal the cold tolerance mechanism.
Physiological and biochemical changes in cold stress are often used as indexes to assess cold tolerance of plants. Antioxidant enzymes such as superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) help eliminate ROS accumulated excessively under abiotic stress to prevent plants from oxidative damage [8]. Activities of these antioxidant enzymes are often used as physiological indices for cold tolerance [9]. However, the severity of membrane damage and ROS burst levels reflect the degree of oxidative damage under stress and are usually used as negative indicators for evaluating cold tolerance [10,11]. Malondialdehyde (MDA) is the ultimate product of lipid peroxidation, and its dynamic accumulation in plant cells indicates the degree of membrane damage. ROS including hydrogen peroxide, hydroxyl radicals, and peroxides not only act as second messengers of cold stress, but their excessive accumulation in cells can lead to oxidative damage to plants. Additionally, soluble sugars and proline act as osmotic substances to protect plants from cold stress [12]. These biochemical indexes can help to understand the cold tolerance mechanism of plants.
In recent years, some cold tolerance mechanisms and genes have been discovered. For example, the ICE-CBF-COR transcriptional cascade is thought to be induced by cold stress and plays a vital role in cold response and adaptation of many plant species [13,14]. It is well known that CBFs, as AP2/ERF transcription factors, can bind to the CRT/DREB elements on promoters of cold-regulated (COR) genes to activate their expression. Overexpression of AtCBF1 and AtCBF3 in tomato, maize, tobacco, and rice enhanced cold tolerance [15,16,17,18]. Many cold-regulated genes have been reported to be related to plant cold tolerance. For instance, overexpression of COR15a in Arabidopsis could enhance the freezing resistance of chloroplasts and protoplasts [19]. Besides, some members of the ERF, bHLH, MYB, and C2H2 transcription factor families, such as ERF105, ICE1, MYB4, MYB15, ZAT10, and ZAT12, are also involved in regulating the expression of genes responding to cold stress [20,21,22]. Previous studies have shown that some phytohormones and their signal transduction-related genes are also crucial in response to cold stress [23]. However, the cold tolerance mechanism of cotton cotyledon remains unknown.
In plants, some candidate genes or loci have been identified by quantitative trait loci (QTL) mapping and genome-wide association analysis [24,25,26,27]. However, it is still challenging to mine cold resistance genes according to the genetic mapping of QTLs for some plant species. Fortunately, with the development of molecular biology techniques, RNA-seq is becoming a popular strategy to identify a large number of cold-responsive genes quickly and to elucidate the molecular mechanisms of cold tolerance in different plant species and organs at different development stages. In wheat, RNA-seq was used to study the cold response mechanisms of meristematic crown and leaves of spring and winter wheat seedlings, and flag leaves of durum wheat at the reproductive stage [28,29]. In rice, transcriptomic analysis has been used to study the cold response mechanisms of leaves, roots, and shoots at the seedling stage, as well as anthers at the reproductive stage [30,31,32]. Recently, scientists have studied the cold response mechanisms of true cotton leaves at the seedling stage and seed embryos at the germination stage by using RNA-seq technology [3,33]. However, reports on the transcriptome profiling of cotton seedling cotyledons under cold stress had not been disclosed.
Young cotyledons, as the first two opposite leaves, are an indispensable part of early cotton seedlings. They are not only the essential autotrophic organs before the first true leaf appeared but also the most vulnerable parts to cold stress. Previous studies on the response of cotton to cold stress mainly focused on the germination stage of seeds and the true leaf stage of seedlings [3,33]. However, the seedling cotyledons are more susceptible to cold stress than true leaves. Therefore, we focused our research on the cold tolerance of young cotyledons at the early seedling stage. In this study, to better understand the diversity in cotyledon cold tolerance of natural populations and to select valuable varieties, we first used the chilling injury (CI) index to evaluate the cold sensitivity of young cotyledons of 74 upland cotton varieties and selected two varieties (H559 and YM21) with contrasting phenotypes for subsequent studies. To comprehensively understand the molecular mechanism of cold tolerance in early seedlings, we performed physiological experiments and comparative transcriptome sequencing on the young cotyledons of both varieties. Through comparative transcriptome analysis, important expression patterns, significant pathways, and hub genes were identified, which would provide help for future research on the molecular mechanism of cotton cotyledon cold tolerance and breeding.

2. Results

2.1. Evaluation of Cold Tolerance Based on CI Index

We evaluated the chilling sensitivity of 74 upland cotton varieties from different ecological regions in China. The CI index of young cotyledons was adopted to measure the cold sensitivity of cotton seedlings. Seedlings on the seventh day after sowing were subjected to a chilling treatment regime of 36 h at 4 °C plus a week at 25 °C. After surveying the chilling injury phenotype and calculating the CI index of each variety, we found that there was a wide variety of cold tolerance among these varieties (Figure 1A). Of these varieties, YM21 appeared to be the most chilling-sensitive (CI index = 43.17%), and extensive cotyledon damage of most young seedlings was observed after seven-day recovery. In contrast, H559 was the most cold-tolerant (CI index = 4.47%), and its cotyledons showed minor damage after recovery (Figure 1A,B). Furthermore, the seedling survival rate of the two varieties after 48 h of cold stress (4 °C) was also investigated (Figure 1C). Two days of continuous chilling stress resulted in the death of a large number of YM21 seedlings (survival rate = 15.24%) while H559 was less affected (survival rate = 90.20%). Therefore, YM21 (chilling-sensitive) and H559 (chilling tolerant), two varieties with distinctively different chilling sensitivity, were selected for further analysis.

2.2. Physiological Responses to Chilling Stress in Cotton Cotyledons

To reveal the physiological response of young cotyledons of both varieties at different time points under cold stress, we determined the activities of SOD, POD, and CAT (Figure 2). When seedlings exposed to 4 °C for 24 h, the activities of the three enzymes in H559 increased at 3 h and then decreased, while the changing trend in YM21 was the opposite. In general, the enzyme activities in YM21 were higher than that of H559 during the cold stress process, except for CAT activity at 3 h. Obviously, there was a conflict between the results of enzyme activity tests and the CI evaluation.
To eliminate this misunderstanding, we measured the ROS (O2− and H2O2) levels and ROS scavenging capability in cotyledons of the two varieties in the cold (Figure 2). With the extension of cold stress time, the level of O2− increased slightly, but it was significantly lower in H559 than that of YM21. In the cotyledons of H559, the content of H2O2 decreased continuously and was always substantially lower than that of YM21. Compared with the control, H2O2 content in YM21 significantly increased when exposed to 4 °C for 24 h. In response to chilling stress, the superoxide anion scavenging (SAC) rate in YM21 was consistently lower than that in H559. These results suggest that H559 is more cold-tolerant in terms of ROS accumulation levels. MDA is the final product of membrane lipid peroxidation and is often used as a physiological index to measure the degree of cold damage. We found that the MDA level in H559 was significantly higher than that in YM21 at 3 h after cold stress, but it was the opposite at 24 h, indicating that the cold tolerance of H559 gradually increased along with the extension of cold stress (Figure 2).
Osmoregulation substances, including soluble sugar (SS), soluble protein (SP), and proline, protect plants from cold stress, and their concentrations were measured here (Figure 2). Under continuous cold stress, the levels of all osmotic adjustments in both varieties showed a rising trend. The content of SS and SP in YM21 was always higher than H559, while the proline content was lower than H559 after cold stress. These results suggest that the cold tolerance mechanism of young cotton cotyledon can be complicated, and changes of many physiological parameters in H559 are superior to YM21 in response to chilling stress.

2.3. Transcriptome Sequencing and Differential Expression Gene Analysis

To reveal the molecular mechanism of the difference in cold tolerance between YM21 and H559 at the early seedling stage, we performed transcriptome sequencing using total RNA extracted from the cotyledons at 0, 3, and 6 h of chilling stress. Eighteen cDNA libraries were constructed and sequenced using Illumina HiSeq2500 by Gene Denovo Biotechnology Co. (Guangzhou, China). Each cDNA library produced about 6.5 Gb to 8.9 Gb of raw data. All RNA-seq raw datasets were deposited in the NCBI database with an SRA accession number SUB7558830. After removing reads containing adapters and low-quality reads, the Q30 of each library ranged from 91.08% to 92.33%, and the GC content ranged from 45.83% to 46.78%. Approximately 92.52% to 97.47% clean reads in each library were mapped to the cotton genome (http://ibi.zju.edu.cn/cotton/), of which 87.6% to 92.27% were uniquely mapped. Besides, 4668 novel genes were identified (Table S1). Details of transcriptome sequencing and alignment with the reference genome were shown in Table 1. The mapped reads of each sample were assembled by using StringTie v. 1.3.1 in a reference-based approach. For each transcription region, an FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify its expression abundance and variations, using StringTie software. A total of 64,461 transcripts were obtained with FPKM values > 0 in at least one sample, of which 44,999 were expressed with FPKM values of ≥1. Therefore, we found that 44,998 genes were expressed in T0, and 344 of these genes were preferentially expressed in T0 (Figure 3A,C). Interestingly, the expressed genes in the T0 sample included all the expressed genes in the others, suggesting that differences in the genetic background might play an essential role in chilling tolerance of the two varieties.
Pearson correlation coefficient analysis based on the expression levels of all expressed genes showed that correlation coefficients between biological repeats ranged from 0.98 to 1.00, indicating that the sample repeatability was reliable (Figure S1). Cluster analysis showed that there was a close correlation between the precold stress and early response stages, and they were different from the samples at 24 h after cold stress. This indicates that the gene expression was less changed in the initial response stage of chilling stress.
Gene differential expression analysis was performed by using DESeq2 software between two different groups with the threshold of false discovery rate (FDR) < 0.05 and an absolute fold change of ≥2 [34]. In H559, 1626 DEGs were found in T0-vs-T3 (954 upregulated, 672 downregulated), 11,902 DEGs in T3-vs-T24 (4224 upregulated, 7678 downregulated), and 13,558 DEGs in T0-vs-T24 (4591 upregulated, 8967 downregulated; Table 2). In YM21, 2030 DEGs were found in S0-vs-S3 (1056 upregulated, 972 downregulated), 8984 DEGs in S3-vs-S24 (4152 upregulated, 4832 downregulated), and 11,795 DEGs in S0-vs-S24 (4522 upregulated, 7273 downregulated; Table 2).

2.4. Genetic Difference between Varieties at the Transcriptomic Level

Previous studies have reported that transcription differences among varieties before cold stress may be related to their genetic differences in response to cold stress [35]. Before cold stress, 44,999 genes expressed in H559 included all the genes (40,067) expressed in YM21. In addition, 1504 DEGs were detected between S0 and T0, of which 964 and 540 DEGs were upregulated in H559 and YM21, respectively (Figure 3B,C and Table 2). These findings suggested that H559 had a more complex background of transcriptional regulation to cope with cod stress than YM21. To further investigate the functional differences of these DEGs, GO-term enrichment analysis was performed (Table S2). For the 964 DEGs upregulated in H559, 84 GO-terms were significantly enriched (p < 0.05) and the top three abundant terms in the biological process category were “single-organism metabolic process”, “response to stimulus”, and “carbohydrate metabolic process”. For the 540 DEGs upregulated in YM21, 77 GO-terms were significantly enriched (p < 0.05) and the terms “single-organism process”, “single-organism metabolic process”, and “single-organism cellular process” were the most abundant in the biological process. However, the GO-term “single-organism metabolic process” was the only one shared by both varieties. To remove the functional redundancy among GO-terms, we obtained overrepresented GO-terms by using the REVIGO program (http://revigo.irb.hr/). In H559, the top three significant clusters related to biological processes included “protein-DNA complex subunit organization”, “regulation of intracellular signal transduction”, and “malate transport” (Figure 4A). However, they were “pigment biosynthetic process”,“pigment metabolic process”, and “single-organism biosynthetic process” in YM21 (Figure 4B). These results suggest that the genes preferentially expressed in the two varieties have significant functional differences under suitable temperature, and H559 has a more sophisticated strategy than YM21 to deal with the possible sudden cold stress.

2.5. Weighted Gene Coexpression Network Analysis

To reveal the differences in gene regulation of cold response in cotton varieties with contrasting cold tolerance, weighted gene coexpression network analysis (WGCNA) was performed using 19,982 DEGs, and 11 merged coexpression gene modules were identified (Figure 5A). Module–sample relationship analysis found that four of the 11 modules had closer and significant correlations with particular samples (Figure 5C). Among them, the saddlebrown module contained 103 genes, which had a significant and positive correlation with T3 (r = 0.95, p = 2 × 10−9). Blue module containing 3743 genes had a closer correlation with T24 (r = 0.87, p = 2 × 10−6). This indicated that the genes of these two modules were related to the cold response of H559. In contrast, the paleturquoise and green modules had significant positive correlations with S3 (r = 0.86, p = 5 × 10−6) and S24 (r = 0.87, p = 2 × 10−6), respectively, suggesting that genes of these two modules were related to cold response of YM21. Additionally, a module–trait relationship analysis was performed using module eigengene and physiological data. As shown in Figure 5D, the saddlebrown and blue modules were highly positively correlated with MDA and proline, respectively. Previous studies have found that MDA mainly reflects the damage degree of cell membrane caused by abiotic stress, and an increase in proline content means enhanced cold resistance. The paleturquoise module was positively correlated with O2−, SS, and SP, but negatively correlated with SAC. It is known that level changes of O2− (a kind of ROS) reflect the ability of cells to resist stress, while SS, SP, and SAC are positively correlated with cold tolerance of plants. There was a significant positive correlation between the green module and SP. Furthermore, we found that the darkmagenta module was highly positively correlated with SAC and that it was distinctively correlated with samples of contrasting genotypes. These results revealed the physiological and transcriptional differences between the cotyledons of H559 and YM21 in response to cold stress.

2.6. GO Enrichment Analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Annotation

To reveal the functional differences of DEGs in the four significant modules, we performed GO-term enrichment analysis and genes in each module were recruited into three categories, including biological process, molecular functions, and cellular component (Table S3). To get nonredundant GO-terms, de-redundancy analysis on the significant GO-terms (p < 0.01) was performed by using the REVIGO program [36]. Overall, 181 nonredundant and significant GO-terms were represented (Figure 6). In the blue module, “response to organonitrogen compound”, “defense response”, and “regulation of biological process” were the top three significant GO-terms for the biological process; “protein kinase activity”, “phosphotransferase activity”, and “kinase activity” were the top three significant GO-terms for the molecular function; whereas only “nuclear ubiquitin ligase complex” was significant for the cellular component. In the saddlebrown module, “glutamine family amino acid metabolic process”, “glutamine metabolic process”, and “response to stress” were the top three significant GO-terms for the biological process; “tetrapyrrole binding” was the most significant GO-term for the molecular function; whereas only “cell periphery” was significant for the cellular component. For the biological process, “response to stress” and “response to stimulus” were the common GO-terms enriched in both modules (blue and saddlebrown). Most terms in the two modules were related to plants resisting abiotic stress. In the green module, “regulation of gene expression”, “regulation of metabolic process”, and “regulation of biological process” were the top three significant GO-terms for the biological process; “nucleic acid binding transcription factor activity”, “DNA binding”, and “oxidoreductase activity” were the top three significant GO-terms for the molecular function; “membrane-bounded organelle” and “intracellular membrane-bounded organelle” were the most significant GO-terms for the cellular component. This suggested that transcription factors (TFs) played an essential role in response to chilling stress in cold-sensitive variety. In the paleturquoise module, “alpha-amino acid metabolic process”, “single-organism catabolic process”, and “CTP biosynthetic process” were the top three significant GO-terms for the biological process; “ligase activity”, “hydrolase activity”, and “anion binding” were the top three significant terms for the molecular function; whereas “cell wall” was the most significant for the cellular component.
To understand the function of genes in different modules more comprehensively, we also conducted KEGG pathway enrichment analysis. As shown in Table 3, a total of 24 significant pathways were obtained from the four significant modules. Among them, the plant hormone signal transduction pathway was shared by three modules of blue, green, and paleturquoise, and MAPK signaling pathway was shared by blue and green modules, indicating that the two pathways played a central role in the process of cotton cotyledon response and resistance to cold stress. In the blue module, the most abundant pathway was “plant–pathogen interaction”. In the saddlebrown module, “phenylpropanoid biosynthesis” and “glutathione metabolism” pathways were specially and significantly enriched. These results indicated that there were both common and unique pathways in cotyledons of the two cotton varieties in response to cold stress.

2.7. Construction of Gene Coexpression Network and Identification of Centrally Connected Genes

Hub genes usually have high connectivity in coexpression networks as indicated by high Kme value. According to the gene expression heatmaps and eigengene histograms of the four significant modules, we discovered that genes in the same module generally had the same preponderant expression stage and that genes in different modules often had different preponderant expression stages (Figure 7A, Figure 8A, Figure S2A and Figure S3A). The predominantly expressed genes at every time point can be assessed by stage specificity score (SSC) [37]. To mine highly expressed genes in preponderant expression stages, we normalized the expression levels of all genes within the same stage by using the Z-score method. The genes with higher Z-score values were considered to have higher expression levels at this specific stage. In our study, genes with SSC > 0.5 and Kme > 0.8 were defined as hub genes (Figure 7B, Figure 8B, Figure S2B and Figure S3B). A total of 662, 186, 39, and 49 hub genes were identified in blue, green, paleturquoise, and saddlebrown modules, respectively. The top five hub genes with high Z-score values in the preponderant expression stage of each module were shown in Table 4.
In the blue module, the top 100 hub genes with higher Z-score values at T24 and weight values >0.45 were used to construct the coexpression network. As shown in Figure 7C, 13 transcription factors were identified, of which ERF4, ERF026, ERF017, and SCL13 were the critical transcription factor genes. Among them, ERF4 (ethylene-responsive element-binding factor 4) was one of the top three hub genes with a higher degree. In the coexpression network of the saddlebrown module, the gene with the highest degree was BH0283 (46 edges), followed by BETV1F (45 edges; Figure 8C). Key genes in the coexpression network of the green module included GH_A10G0211 (hypothetical protein CCACVL1_15597), GH_A05G4197 (sigma factor binding protein 2) and GH_D08G1892 (probable serine/threonine-protein kinase PBL19; Figure S2C). In the coexpression network of paleturquoise module, hub genes were almost correlated with each other, of which TIP1-1 (aquaporin TIP1-1-like) and BZIP1 (basic leucine zipper 1-like) were identified as key genes (Figure S3C).

2.8. Transcription Factors in Response to Chilling Stress

Transcription factors (TFs) play an essential role in plant response to abiotic stress [38]. GO-term analysis of the genes in the blue and green modules revealed that quite a number of them were enriched to “regulation of biological process”, suggesting that TFs participated in regulating cold response. To better understand the transcription differences of TF genes in different genotypes and response stages, we performed TF analysis on the significant modules. A total of 594 TFs were identified, including 401 in the blue module, 155 in the green module, 35 in the paleturquoise module, and two in the saddlebrown module (Table S4). This indicated that more TFs were involved in regulating cold response after exposure to cold stress for 24 h compared with 3 h. In the blue module, the top three abundant TF families were ERF (70 genes, 17.46%), WRKY (41 genes, 10.22%), and MYB (39 genes, 9.73%), while the key TF genes include ZAT10 (GH_D05G2128), RAV1 (GH_A03G1130), and ERF4 (GH_D01G0569; Figure 9 and Table S4). In the green module, the three major TF families were ERF (41 genes, 26.45%), WRKY (16 genes, 10.32%), and NAC (15 genes, 9.68%), and the key TF genes included ERF5 (GH_A12G2350) and ERF1B (GH_D02G0425; Figure 8 and Table S4). There was only one hub TF gene BZIP1 (GH_D09G0794) in the paleturquoise module, but none had been identified in the saddlebrown module.

2.9. Genes Related to Plant Hormone Signal Transduction

Phytohormones play an important role in the response of plants to various abiotic stresses [39]. In our study, 133 genes from blue, green, and paleturquoise modules were significantly enriched in the “plant hormone signal transduction pathway”, and three genes from the saddlebrown module were also detected in this pathway (Figure 10). The differences in expression patterns and expression levels of phytohormone-related genes were shown in Figure 10B. In the auxin (IAA) signal transduction pathway, a total of 46 genes were detected, including two AUX1 genes, 12 AUX/IAA genes, one ARF gene, 13 GH3 genes, and 18 SAUR genes, of which 30 genes were from the blue module. After 24 h of cold stress, most IAA genes were significantly upregulated, especially those in the blue module showed higher expression levels, indicating that the IAA genes might have an important role in regulating the cold response of cotyledons. A total of 20 genes were found in the abscisic acid (ABA) pathway, including two ABF genes, six PP2C genes, eight PYL genes, and four SnPK2 genes. Moreover, a total of 22 genes were found in the brassinosteroid (BR) pathway, 20 of which were from the blue module, including three BAK1 genes, one BRI1 gene, one BSK gene, one BIN2 gene, two CYCD3 genes, and 12 TCH4 genes. These TCH4 genes were all from blue and green modules and had higher expression levels, indicating that they mainly played a role in cotyledon cold adaptation and tolerance. Furthermore, 22 genes were also found in the ethylene (ET) pathway, including three ETR genes, one SIMKK gene (MKK4/5), one EIN2 gene, three EIN3 genes, five EBF1/2 genes, and nine ERF1 genes. In addition, 12 genes were identified in the jasmonic acid (JA) pathway, including one JAR1 gene, four JAZ genes, and seven MYC2 genes. We also found two GID1 genes, three DELLA genes and three PIF genes (TF) in the gibberellin (GA) pathway, and one NPR1 gene, three TGA genes, and one PR1 gene in salicylic acid (SA) pathway. Remarkably, all three genes in the cytokinin (CK) pathway were from the blue module, including one AHP gene and two ARR-A genes.

2.10. Verification of RNA-Seq Data by qRT-PCR

To validate the accuracy and reliability of RNA-seq data and obtain candidate genes for further functional studies, we performed quantitative real-time PCR (qRT-PCR) for 20 interesting hub genes selected from four significant WGCNA modules. The 20 genes and their qRT-PCR primers were listed in Table S5. We found that the relative expression levels of genes obtained by qRT-PCR were highly consistent with the RNA-seq data (Figure 11 and Figure S4). Linear regression analysis showed that the correlation between RNA-seq and qRT-PCR data was significantly positive (Figure 11B,C). These results suggested that the RNA-seq data was credible and accurate.

3. Discussion

Cotton (Gossypium hirsutum L.) is sensitive to chilling stress, especially at seedling emergence and early cotyledon stage. We had previously found that under low natural temperature, the germination and cotyledon stages of cotton would last for a long time. Relevant evidence also supported that low temperatures could reduce germination speed and delay plant growth [40]. Once the ambient temperature is suitable, cotton seedlings will quickly resume growth and development. Therefore, improving cold tolerance of young cotyledons is more meaningful for cotton seedlings. Some researchers have noticed this and assessed the cold sensitivity of young cotyledons of different cotton varieties with the help of CI index [41]. We found that cold stress (4 °C) for 24–36 h could not cause significant damage to roots and stems of most cotton seedlings, but would cause varying injury degrees of cotyledons. We used the revised CI index to evaluate the cold sensitivity of 74 cotton varieties, and found that H559 was the most tolerant and YM21 was the most sensitive. To elucidate the molecular mechanism of cotton cotyledons response and resistant to cold, physiological changes and transcriptomic analysis was performed in our study.

3.1. Physiological Response of Cotton Cotyledons to Chilling Stress

A plant will face various temperature stresses during its lifetime, including freezing, chilling, and high temperature. However, molecular mechanisms of plant responses and resistance to different temperature stresses are usually different [42]. To cope with cold stress, plants have evolved many strategies to help them survive and maintain growth, such as keeping ROS balance and accumulating osmotic substances [43,44]. In plants, ROS include hydrogen peroxide (H2O2), singlet oxygen, and superoxide anions, and their levels in cells are usually regulated by various antioxidant enzymes, such as CAT, POD, and SOD [45]. MDA is the ultimate product of phospholipid peroxidation, and its content reflects the integrity of cell membranes [46]. Under cold stress, MDA content increased rapidly and then decreased to normal levels in H559, indicating that cold-tolerant variety could quickly start the cold defense system to maintain and enhance the integrity of the cell membrane. Besides, the contents of H2O2 and O2− reflect the ROS levels in cells, and their excessive accumulation will disrupt intracellular ROS balance and cause oxidative damage [47]. In our study, the contents of O2− and H2O2 in H559 were significantly lower than YM21 under both normal and low temperatures, and the H2O2 content in H559 showed a reduction trend, indicating that the tolerant variety could restrict the accumulation of ROS in the cold. Under abiotic stress, various osmoregulatory substances such as soluble sugars (SS) and soluble proteins (SP) accumulate rapidly to enhance plant stress resistance [44,48]. Before and within 24 h of cold stress, the intracellular SS and SP concentrations of cold-sensitive variety were consistently and significantly higher than that of tolerant variety, which was similar to previous studies [41,49]. This may indicate that cotyledons of tolerant variety are not highly dependent on osmotic substances to improve cold tolerance in the early stage of cold response. It is well known that the activity of antioxidant enzymes represents ROS removal capability [8]. Compared with YM21, antioxidant enzyme activities increased rapidly in H559 in cold stress to remove over accumulated ROS. Therefore, we believe that these results of enzyme (SOD and POD) activities, which are different from previous research [41], are mainly due to the different mechanisms of cold resistance between cotyledons and true leaves. The cold tolerance of cotton cotyledons may be more dependent on genetic background and transcriptional regulation of the varieties.

3.2. Transcriptional Differences in Cotyledons of Varieties with Different Cold Sensitivity

Transcriptome technology has been used to reveal differences in the molecular mechanisms of cold tolerance among different varieties [33,50,51,52,53]. Previous studies have shown that there are more constitutive expression genes in the tolerant variety than the sensitive one before cold stress treatment [54,55]. It is believed that there should be relatively more genes involved in forming cold-tolerant genotypes. Under normal temperature, 40,067 and 44,999 expressed genes were detected in young cotyledons of YM21 and H559, respectively. It was noteworthy that expressed genes in T0 cotyledons included all the expressed genes of S0 and all the DEGs obtained by comparative analysis. GO enrichment analysis of the up and downregulated expressed genes acquired by comparing S0 with T0 revealed that their biological functions were significantly different, which was consistent with previous findings [5]. Thus, the constitutive diversities between varieties can result in different cold tolerance. In this study, WGCNA based on all DEGs helped us get four significant modules: early response module of YM21 (paleturquoise), late response module of YM21 (green), early response module of H559 (saddlebrown) and late response module of H559 (blue). In terms of the number of genes contained in each module, the paleturquoise module has more genes than the saddlebrown module in the early response stage, while the green module has fewer genes than the blue module in late response stage. Additionally, we found that the saddlebrown module was significantly correlated with MDA, while the paleturquoise module was significantly correlated with SS and SP. These results indicate that the early cold response of H559 does not perform as well as that of YM21. Many previous studies on plant cold tolerance have found that genes related to plant hormone signal transduction and protein kinases are induced by cold stress [56,57,58]. KEGG pathway analysis found that plant hormone signal transduction and MAPK signal pathways were significantly enriched in both tolerant and sensitive varieties at 24 h of cold stress, suggesting that they played important roles in regulating cold tolerance of cotton cotyledons.

3.3. Genes Related to Phytohormone Signal Transduction Involved in Cotyledon Cold Response

Phytohormones mainly include IAA, GA, CK, ABA, ET, BR, SA, JA, and strigolactone, which are small endogenous signaling molecules. Among them, BR, ABA, IAA, ET, GA, and JA have been reported to respond to cold stress and regulate cold resistance [59,60,61]. Previous studies have demonstrated that JA positively regulated plant cold tolerance [62]. In Arabidopsis thaliana, the low temperature can elevate the level of endogenous JA by inducing the expression of JA biosynthetic genes [63]. Moreover, the external application of JA improves cold tolerance in Arabidopsis. In this study, we found that some genes related to JA biosynthesis, including LOXs, AOS, and AOC, and genes associated with JA signal transduction including JAR1, JAZ, and MYC2 were predominantly expressed in T24. This means that the JA signal plays an important role in regulating cotton cotyledon cold tolerance. It is known that BR treatment can improve the cold resistance of plants [64,65]. It was reported that overexpression of BRI1 improved cold tolerance in Arabidopsis, whereas the bri1 mutants showed cold hypersensitivity [66]. Some studies have suggested that TCH4, as a member of the BR signal transduction pathway, can respond to auxin and BR, as well as various environmental stimuli, including cold and heat shock [67,68]. In the present study, we detected one BRI1 gene and 12 TCH4 genes predominantly expressed in T24 and found that most genes in the BR signaling pathway were exclusively present in the blue module, indicating that BR signaling was necessary for cotyledon cold tolerance. Additionally, some studies have demonstrated that genes involved in the auxin signaling pathway are differentially expressed when plants are subjected to cold stress [30,69]. In these remarkable modules, we identified 136 genes involved in plant hormone signaling pathway, of which 46 were related to auxin signal transduction, most of which were from the blue module, indicating that auxin played a significant role in regulating cotyledon cold tolerance.

3.4. Protein Kinase-Related Genes Involved in Cotyledon Cold Response

Plant MAPK cascade participates in signal transduction induced by abiotic stress [70,71,72]. In Arabidopsis, ROS activate MAPK cascades to mediate the transmission of cold signals [73]. This indicates that the MAPK signal transduction pathway is important for plant response to cold stress. Our results demonstrated that the MAPK signaling pathway was significantly enriched in blue and green modules. The MAPK cascade transmits and amplifies signals by phosphorylating target proteins through three phosphorylation kinases (MPKs, MEKs, and MEKKs), and those phosphorylated proteins will alter the expression of downstream genes, including transcription factors and COR genes [26,71,72]. Some studies have shown that the MAPK pathway mediated by MEKK1-MKK2-MPK4/6 positively responds to cold stress [73,74]. In Arabidopsis, MKK2 can be induced by cold stress, and plants overexpressing MKK2 exhibit strong cold tolerance [73]. It is reported that overexpression of OsMKK6 activates OsMPK3 and enhances the cold endurance of rice [75]. However, the single mpk3 mutant of rice showed freezing resistance [76]. These two diametrically opposite results are presumed to be caused by different degrees of low temperature. In our study, both MKK2 and MPK6 were upregulated under chilling stress. Additionally, we also identified that MPK4, MKK4/5, MKK3, MKK9, and MEKK17/18 were highly expressed when H559 was subjected to cold stress, suggesting that these genes might be involved in regulating cold tolerance of H559 cotyledons. It is known that CDPKs and CIPKs are two protein kinases involved in the calcium signal transduction pathway and participate in plant response to cold stress [77,78]. Under low temperature, the transcription of many CDPK and CIPK genes in plants were induced [52,79]. It has been reported that some CIPK and CDPK genes play important roles in improving the cold resistance of plants [80,81,82,83,84,85]. In this study, 18 CIPK genes and 14 CPK genes were detected in blue (15 CIPK and 13 CPK genes) and green (three CIPK genes and one CPK gene) modules, indicating that most CPK and CIPK genes were preferentially expressed in the cold-tolerant variety. These genes will provide important implications for further research on the cold tolerance of cotton cotyledons.

3.5. Major TFs Involved in the Cold Response of Cotton Cotyledons

TFs play an essential role in regulating plant growth and development, as well as during abiotic stress response. When exposed to a low-temperature environment, TFs induced by cold stress are activated by various signaling pathways. The activated TFs can specifically bind to corresponding cis-acting elements of their target gene promoters to activate cold-responsive genes, thereby regulating the cold tolerance of plants [86]. In the present study, the AP2/EFR family contained most members, followed by WRKY and MYB families, suggesting that they played an important role in the response of cotton cotyledons to chilling stress. DREB is a subfamily of ERF TF family, and its members DREB1B/CBF1, DREB1C/CBF2, and DREB1A/CBF3 are induced by cold stress [87]. Overexpression of DREBs/CBFs in Arabidopsis improved cold tolerance of plants, while suppression of DREB1B or DREB1A reduced plant cold tolerance [88,89,90]. In our study, two DREB1B/CBF1 genes (GH_A12G2453 and GH_D12G2465) were identified, of which GH_A12G2453 had higher expression levels induced by chilling stress in H559. Similarly, some members of WRKY and MYB families, such as VaWRKY12, VbWRKY32, AtMYB14, and MdMYB23 have been reported to participate in regulating cold tolerance of plants [91,92,93,94]. In addition, some other TFs families were upregulated in the cold-tolerant variety, such as C2H2, NAC, GRAS, bHLH, and bZIP. It is reported that ZAT6, ZAT10, and ZAT12, as members of the C2H2 family, participate in cold response [95]. In our study, a ZAT10 gene was identified as a hub gene in the blue module, with higher expression in H559. Excitingly, most of the differentially expressed TF genes in the tolerant variety were involved in response to cold stress, and they would provide promising prospects for the improvement of cold tolerance in cotton cotyledons.

4. Materials and Methods

4.1. Plant Materials and Cold Treatments

Seventy-four upland cotton varieties (provided by the Institute of Cotton Research of CAAS) were used to evaluate cotyledon cold tolerance. When the seedlings grew to the 6th day after sowing, plants with deformed cotyledons or weak growth were pulled out. Three biological replicates were set for each variety, and 50 to 60 seedlings were retained for each replicate. On the seventh day, plug seedlings were transferred to cold light plant growth chamber (15 h light/9 h dark), and after 36 h of cold stress (4 °C), they were returned to the cotton cultured room (15 h light/9 h dark, 25 °C) for another seven days. After a week of recovery, the chilling injury phenotype was investigated, and cold sensitivity was evaluated with the CI index. The cold-tolerant variety H559 and sensitive variety YM21 (also known as Yumian21) were selected from the 74 varieties. The seventh-day seedlings of both varieties were subjected to 4 °C for 0 to 24 h. Three replicates were set for each variety at each time point, and six cotyledons (from 6 seedlings) were collected for each replicate. T0, T3, and T24 represent cotyledon samples of H559 under cold stress at 0 h, 3 h, and 24 h, respectively. Similarly, S0, S3, and S24 represent cotyledon samples of YM21 under cold stress at 0 h, 3 h, and 24 h, respectively. The collected samples (T0 to T24 and S0 to S24) were used for physiological indices determination and RNA-seq. At the same time, seedlings of the two varieties were cold-treated for 48 h to investigate the survival rate after 14 days of recovery, and the emergence of true leaves was used as the criterion for seedling survival.

4.2. CI Index Evaluation and Investigation of Seedling Survival Rate

To evaluate the CI index of different varieties, we first graded cold-damage symptoms of single cotyledon according to the necrotic area (X): 0 (healthy), 1 (yellowish), 2 (X < 1/8), 3 (1/8 ≤ X < 1/4), 4 (1/4 ≤ X < 1/2), 5 (1/2 ≤ X < 3/4), and 6 (X ≥ 3/4). Then, the cold-damage symptoms of the cotyledons of individual seedling were graded according to the necrotic area (Y): 0 (Y = 0), 1 (Y < 1/16), 2 (1/16 < Y ≤ 1/8), 3 (1/8 < Y ≤ 1/4), 4 (1/4 < Y ≤ 3/8), 5 (3/8 < Y ≤ 1/2), 6 (1/2 < Y ≤ 3/4), and 7 (Y > 3/4). The CI index for a specific variety is denoted as: CI = 100 × (0 × S0 + 1 × S1 + 2 × S2 + 3 × S3 + 4 × S4 + 5 × S5 + 6 × S6 + 7 × S7)/8N, where Si represents the number of seedlings at grade i, and N represents the total number of seedlings. Therefore, the CI value ranged from 0 to 1, and the lower the CI value, the more cold-resistant the variety. The CI evaluation method used here was an improvement and optimization based on a previous study [41].

4.3. Determination of Physiological Indices

Young cotyledons were sampled at 0, 3, and 24 h after cold stress and frozen quickly in liquid nitrogen. All samples were prepared for three biological replicates. All physiological indices, including superoxide dismutase (SOD) activity, peroxidase (POD) activity, catalase (CAT) activity, malondialdehyde (MDA) content, soluble sugar (SS) content, soluble protein (SP) content, superoxide anion (O2−) content, H2O2 content, proline content, and superoxide anion scavenging (SAC) rate, were determined using assay kits (Comin, Suzhou, China). The SPSS statistical software (SPSS version 19.0) was used for the analysis of variance.

4.4. RNA Extraction, Transcriptome Sequencing, and Data Analysis

Young cotyledons of H559 and YM21 seedlings were collected after 0, 3, and 24 h of cold stress and quickly frozen in liquid nitrogen. Each sample included three biological replicates, and each replicate consisted of six young cotyledons collected from six seedlings. Total RNA for each sample was extracted according to the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) guidelines. RNA quality was examined with an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase free agarose gel electrophoresis. Preparation of 18 cDNA libraries and transcriptome sequencing using the Illumina HiSeq2500 platform were completed by Gene Denovo Biotechnology Co. (Guangzhou, China). After raw reads were filtered, the obtained clean reads were aligned to the latest version of the Gossypium hirsutum genome, which can be available at http://ibi.zju.edu.cn/cotton. The FPKM value was calculated using StringTie v. 1.3.1 software to assess the abundance of each transcript [96]. Pearson correlation coefficient analysis was performed with all FPKM > 0 genes. The DESeq2 software was used to analyze the DEGs, and the genes with the parameter of FDR < 0.05 and an absolute fold change of ≥2 were considered differentially expressed [34]. For the specific comparison of DEGs, GO enrichment, and REVIGO (http://revigo.irb.hr) deredundancy analysis were performed [36].

4.5. Weighted Gene Coexpression Network Analysis

A total of 19,982 DEGs were used for WGCNA analysis. Coexpression modules were discovered using WGCNA (v. 1.47) package in R software with default settings, except that the power was 12, mergeCutHeight was 0.75, and minModuleSize was 50. To find out biologically significant modules, we used module eigengenes to calculate correlation coefficients with samples or physiological traits. Intramodular connectivity and stage specificity score (SSC) of each gene were calculated, and genes with higher SSC and connectivity tended to be hub genes. The SSC score was calculated as previously reported by Zhan et al. [37]. The networks were visualized using Gephi 0.9.2 [97]. For genes in each module, GO-term and KEGG pathway enrichment analysis were conducted to reveal the biological functions of modules. The REVIGO program (http://revigo.irb.hr/) was used to remove redundant GO-terms, and heatmaps were used to visualize nonredundant GO-terms.

4.6. RNA-Seq Data Validation

Young cotyledons were sampled at 0, 3, 6, 12, and 24 h after cold stress and frozen quickly in liquid nitrogen. The total RNA of each sample was extracted according to the instructions of RNAprep Pure Plant Plus Kit (TIANGEN Biotech Co.,Ltd., Beijing, China). All qRT-PCR primers were designed using Oligo7 software and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). The GhACTIN7 (AY305733) gene was used as an internal reference gene [98]. According to the system and procedures recommended by the UltraSYBR Mixture (Beijing ComWin Biotech Co.,Ltd., Beijing, China) instructions, qRT-PCR was performed on a 7500 real-time PCR system. The relative expression of genes was analyzed by the 2−ΔΔCt method. The verification experiment was performed in three biological replicates, with three technical replicates in each replicate.

5. Conclusions

In this study, we successfully identified outstanding cold-tolerant varieties through CI evaluation of young cotyledons after chilling stress and found that ROS levels and ROS scavenging capability were reliable physiological indices that could best reflect the cold resistance of cotyledons. Transcriptome data revealed that tolerant cotyledons had more sophisticated cold-induced transcriptional regulation capabilities than sensitive cotyledons before and after cold stress. WGCNA provided us with gene expression modules (blue, green, saddlebrown, and paleturquoise) and functional pathways (plant hormone signal transduction, MAPK signal, and plant–pathogen interaction) that positively respond to cold stress. These findings could not only help us understand the potential molecular mechanisms of cold tolerance in cotton cotyledons, but also laid a foundation for future research on the roles of identified hub genes in the regulation network of cold tolerance and breeding of cold-resistant varieties.

Supplementary Materials

Supplementary materials can be found at https://www.mdpi.com/1422-0067/21/14/5095/s1. Table S1: New genes, Figure S1: Pearson correlation coefficient analysis of transcriptome data from cotyledons of two cotton varieties at different time points under chilling stress, Table S2: Significant GO-terms in YM21 and H559 before cold stress, Table S3: Significant GO-term enrichment of four significant WGCNA modules, Figure S2: Coexpression network construction of green module, Figure S3: Coexpression network construction of paleturquoise module, Table S4: Transcription factor genes in four significant WGCNA modules, Table S5: Primer sequences used in qRT-PCR, Figure S4: Clustering heatmap of the hub gene expression levels revealed using RNA-seq and qRT-PCR. All authors have read and agreed to the published version of the manuscript.

Author Contributions

Conceptualization, H.W. (Hengling Wei) and S.Y.; data curation, G.C.; formal analysis, G.C. and L.Z.; funding acquisition, S.Y.; investigation, G.C. and L.Z.; mthodology, H.W. (Hengling Wei); validation, G.C. and L.Z.; visualization, G.C.; writing—original draft, G.C.; writing—review and editing, L.Z., H.W. (Hantao Wang) and J.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Shandong TAISHAN Industry leading talent program (LJNY201608).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

ABAAbscisic acid
BRBrassinosteroid
CATCatalase
CKCytokinin
DEGDifferentially expressed gene
ETEthylene
FPKMFragment per kilobase of transcript per million mapped reads
GAGibberellin
GOGene ontology
IAAAuxin
JAJasmonic acid
KEGGKyoto Encyclopedia of Genes and Genomes
MDAMalondialdehyde
PCAPrincipal component analysis
PODPeroxidase
ROSReactive oxygen species
SASalicylic acid
SACSuperoxide anion scavenging
SODSuperoxide dismutase
SPSoluble protein
SSSoluble sugar
SSCStage specificity score
TFTranscription factor
WGCNAWeighted gene coexpression network analysis

References

  1. Wang, Q.; Liu, N.; Yang, X.; Tu, L.; Zhang, X. Small RNA-mediated responses to low- and high-temperature stresses in cotton. Sci. Rep. 2016, 6, 35558. [Google Scholar] [CrossRef] [PubMed]
  2. Guo, X.; Xu, S.; Chong, K. Cold Signal Shuttles from Membrane to Nucleus. Mol. Cell 2017, 66, 7–8. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Li, Z.B.; Zeng, X.Y.; Xu, J.W.; Zhao, R.H.; Wei, Y.N. Transcriptomic profiling of cotton Gossypium hirsutum challenged with low-temperature gradients stress. Sci. Data 2019, 6, 197. [Google Scholar] [CrossRef] [PubMed]
  4. Zhang, T.; Zhao, X.; Wang, W.; Pan, Y.; Huang, L.; Liu, X.; Zong, Y.; Zhu, L.; Yang, D.; Fu, B. Comparative transcriptome profiling of chilling stress responsiveness in two contrasting rice genotypes. PLoS ONE 2012, 7, e43274. [Google Scholar] [CrossRef] [Green Version]
  5. Li, Z.; Hu, G.; Liu, X.; Zhou, Y.; Li, Y.; Zhang, X.; Yuan, X.; Zhang, Q.; Yang, D.; Wang, T.; et al. Transcriptome Sequencing identified genes and gene ontologies associated with early freezing tolerance in maize. Front. Plant Sci. 2016, 7, 1477. [Google Scholar] [CrossRef] [Green Version]
  6. Leipner, J.; Stamp, P. Chilling Stress in Maize Seedlings. In Handbook of Maize: Its Biology; Bennetzen, J.L., Hake, S.C., Eds.; Springer: Berlin/Heidelberg, Germany, 2009; pp. 291–310. [Google Scholar]
  7. Huang, J.; Zhang, J.; Li, W.; Hu, W.; Duan, L.; Feng, Y.; Qiu, F.; Yue, B. Genome-wide association analysis of ten chilling tolerance indices at the germination and seedling stages in maize. J. Integr. Plant Biol. 2013, 55, 735–744. [Google Scholar] [CrossRef]
  8. Mutlu, S.; Karadağoğlu, Ö; Atici, Ö; Nalbantolu, B. Protective role of salicylic acid applied before cold stress on antioxidative system and protein patterns in barley apoplast. Biol. Plant. 2013, 57, 507–513. [Google Scholar] [CrossRef]
  9. Han, B.; Ma, X.; Cui, D.; Wang, Y.; Geng, L.; Cao, G.; Zhang, H.; Han, L. Comprehensive Evaluation and Analysis of the Mechanism of Cold Tolerance Based on the Transcriptome of Weedy Rice Seedlings. Rice 2020, 13, 12. [Google Scholar] [CrossRef]
  10. Anwar, A.; Bai, L.; Miao, L.; Liu, Y.; Li, S.; Yu, X.; Li, Y. 24-Epibrassinolide ameliorates endogenous hormone levels to enhance low-temperature stress tolerance in cucumber seedlings. Int. J. Mol. Sci. 2018, 19, 2497. [Google Scholar] [CrossRef] [Green Version]
  11. Chen, Z.Y.; Wang, Y.T.; Pan, X.B.; Xi, Z.M. Amelioration of cold-induced oxidative stress by exogenous 24-epibrassinolide treatment in grapevine seedlings: Toward regulating the ascorbate–glutathione cycle. Sci. Hortic. 2019, 244, 379–387. [Google Scholar] [CrossRef]
  12. Ruelland, E.; Vaultier, M.N.; Zachowski, A.; Hurry, V. Cold Signalling and Cold Acclimation in Plants. Adv. Bot. Res. 2009, 49, 35–150. [Google Scholar]
  13. Guo, J.; Ren, Y.; Tang, Z.; Shi, W.; Zhou, M. Characterization and expression profiling of the ICE-CBF-COR genes in wheat. PeerJ 2019, 7, e8190. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Wang, D.Z.; Jin, Y.N.; Ding, X.H.; Wang, W.J.; Zhai, S.S.; Bai, L.P.; Guo, Z.F. Gene Regulation and Signal Transduction in the ICE-CBF-COR Signaling Pathway during Cold Stress in Plants. Biochemistry 2017, 82, 1103–1117. [Google Scholar] [CrossRef] [PubMed]
  15. Hsieh, T.H.; Lee, J.T.; Yang, P.T.; Chiu, L.H.; Charng, Y.Y.; Wang, Y.C.; Chan, M.T. Heterology expression of the Arabidopsis C-repeat/dehydration response element binding factor 1 gene confers elevated tolerance to chilling and oxidative stresses in transgenic tomato. Plant Physiol. 2002, 129, 1086–1094. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Qin, F.; Sakuma, Y.; Li, J.; Liu, Q.; Li, Y.Q.; Shinozaki, K.; Yamaguchi-Shinozaki, K. Cloning and functional analysis of a novel DREB1/CBF transcription factor involved in cold-responsive gene expression in Zea mays L. Plant Cell Physiol. 2004, 45, 1042–1052. [Google Scholar] [CrossRef] [Green Version]
  17. Kasuga, M.; Miura, S.; Shinozaki, K.; Yamaguchi-Shinozaki, K. A combination of the Arabidopsis DREB1A gene and stress-inducible rd29A promoter improved drought- and low-temperature stress tolerance in tobacco by gene transfer. Plant Cell Physiol. 2004, 45, 346–350. [Google Scholar] [CrossRef] [Green Version]
  18. Oh, S.J.; Song, S.I.; Kim, Y.S.; Jang, H.J.; Kim, S.Y.; Kim, M.; Kim, Y.K.; Nahm, B.H.; Kim, J.K. Arabidopsis CBF3/DREB1A and ABF3 in transgenic rice increased tolerance to abiotic stress without stunting growth. Plant Physiol. 2005, 138, 341–351. [Google Scholar] [CrossRef] [Green Version]
  19. Artus, N.N.; Uemura, M.; Steponkus, P.L.; Gilmour, S.J.; Lin, C.; Thomashow, M.F. Constitutive expression of the cold-regulated Arabidopsis thaliana COR15a gene affects both chloroplast and protoplast freezing tolerance. Proc. Natl. Acad. Sci. USA 1996, 93, 13404–13409. [Google Scholar] [CrossRef] [Green Version]
  20. Bolt, S.; Zuther, E.; Zintl, S.; Hincha, D.K.; Schmülling, T. ERF105 is a transcription factor gene of Arabidopsis thaliana required for freezing tolerance and cold acclimation. Plant. Cell Environ. 2017, 40, 108–120. [Google Scholar] [CrossRef]
  21. Tang, K.; Zhao, L.; Ren, Y.; Yang, S.; Zhu, J.K.; Zhao, C. The transcription factor ICE1 functions in cold stress response by binding to the promoters of CBF and COR genes. J. Integr. Plant Biol. 2020, 62, 258–263. [Google Scholar] [CrossRef] [Green Version]
  22. Chinnusamy, V.; Zhu, J.; Zhu, J.K. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007, 12, 444–451. [Google Scholar] [CrossRef] [PubMed]
  23. Ku, Y.S.; Sintaha, M.; Cheung, M.Y.; Lam, H.M. Plant Hormone Signaling Crosstalks between Biotic and Abiotic Stress Responses. Int. J. Mol. Sci. 2018, 19, 3206. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Deng, X.; Gan, L.; Liu, Y.; Luo, A.; Jin, L.; Chen, J.; Tang, R.; Lei, L.; Tang, J.; Zhang, J.; et al. Locating QTLs controlling overwintering seedling rate in perennial glutinous rice 89-1 (Oryza sativa L.). Genes Genom. 2018, 40, 1351–1361. [Google Scholar] [CrossRef] [PubMed]
  25. Kruse, E.B.; Carle, S.W.; Wen, N.; Skinner, D.Z.; Murray, T.D.; Garland-Campbell, K.A.; Carter, A.H. Genomic Regions Associated with Tolerance to Freezing Stress and Snow Mold in Winter Wheat. G3 Genes Genomes Genet. 2017, 7, 775–780. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Hu, G.; Li, Z.; Lu, Y.; Li, C.; Gong, S.; Yan, S.; Li, G.; Wang, M.; Ren, H.; Guan, H.; et al. Genome-wide association study Identified multiple Genetic Loci on Chilling Resistance During Germination in Maize. Sci. Rep. 2017, 7, 10840. [Google Scholar] [CrossRef]
  27. Thapa, R.; Tabien, R.E.; Thomson, M.J.; Septiningsih, E.M. Genome-Wide Association Mapping to Identify Genetic Loci for Cold Tolerance and Cold Recovery During Germination in Rice. Front. Genet. 2020, 11, 22. [Google Scholar] [CrossRef] [Green Version]
  28. Gulick, P.J.; Drouin, S.; Yu, Z.; Danyluk, J.; Poisson, G.; Monroy, A.F.; Sarhan, F. Transcriptome comparison of winter and spring wheat responding to low temperature. Genome 2005, 48, 913–923. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Diaz, M.L.; Soresi, D.S.; Basualdo, J.; Cuppari, S.J.; Carrera, A. Transcriptomic response of durum wheat to cold stress at reproductive stage. Mol. Biol. Rep. 2019, 46, 2427–2445. [Google Scholar] [CrossRef]
  30. Guan, S.; Xu, Q.; Ma, D.; Zhang, W.; Xu, Z.; Zhao, M.; Guo, Z. Transcriptomics profiling in response to cold stress in cultivated rice and weedy rice. Gene 2019, 685, 96–105. [Google Scholar] [CrossRef] [PubMed]
  31. Yang, Y.W.; Chen, H.C.; Jen, W.F.; Liu, L.Y.; Chang, M.C. Comparative Transcriptome Analysis of Shoots and Roots of TNG67 and TCN1 Rice Seedlings under Cold Stress and Following Subsequent Recovery: Insights into Metabolic Pathways, Phytohormones, and Transcription Factors. PLoS ONE 2015, 10, e0131391. [Google Scholar] [CrossRef]
  32. Bai, B.; Wu, J.; Sheng, W.T.; Zhou, B.; Zhou, L.J.; Zhuang, W.; Yao, D.P.; Deng, Q.Y. Comparative Analysis of Anther Transcriptome Profiles of Two Different Rice Male Sterile Lines Genotypes under Cold Stress. Int. J. Mol. Sci. 2015, 16, 11398–11416. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Shen, Q.; Zhang, S.; Liu, S.; Chen, J.; Ma, H.; Cui, Z.; Zhang, X.; Ge, C.; Liu, R.; Li, Y.; et al. Comparative Transcriptome Analysis Provides Insights into the Seed Germination in Cotton in Response to Chilling Stress. Int. J. Mol. Sci. 2020, 21, 2067. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Liu, S.; Ge, F.; Huang, W.; Lightfoot, D.A.; Peng, D. Effective identification of soybean candidate genes involved in resistance to soybean cyst nematode via direct whole genome re-sequencing of two segregating mutants. Theor. Appl. Genet. 2019, 132, 2677–2687. [Google Scholar] [CrossRef]
  36. Supek, F.; Bosnjak, M.; Skunca, N.; Smuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE 2011, 6, e21800. [Google Scholar] [CrossRef] [Green Version]
  37. Zhan, J.; Thakare, D.; Ma, C.; Lloyd, A.; Nixon, N.M.; Arakaki, A.M.; Burnett, W.J.; Logan, K.O.; Wang, D.; Wang, X.; et al. RNA sequencing of laser-capture microdissected compartments of the maize kernel identifies regulatory modules associated with endosperm cell differentiation. Plant Cell 2015, 27, 513–531. [Google Scholar] [CrossRef]
  38. Chinnusamy, V.; Zhu, J.K.; Sunkar, R. Gene regulation during cold stress acclimation in plants. Methods Mol. Biol. 2010, 639, 39–55. [Google Scholar]
  39. Verma, V.; Ravindran, P.; Kumar, P.P. Plant hormone-mediated regulation of stress responses. BMC Plant Biol. 2016, 16, 86. [Google Scholar] [CrossRef] [Green Version]
  40. Hussain, S.; Khan, F.; Hussain, H.A.; Nie, L. Physiological and Biochemical Mechanisms of Seed Priming-Induced Chilling Tolerance in Rice Cultivars. Front. Plant Sci. 2016, 7, 116. [Google Scholar] [CrossRef] [Green Version]
  41. Wang, J.J.; Wang, D.L.; Yin, Z.J.; Wang, S.; Fan, W.L.; Lu, X.K.; Mu, M.; Guo, L.X.; Ye, W.W.; Yu, S.X. Identification of the chilling resistance from germination stage to seedling stage in upland cotton. Sci. Agric. Sin. 2016, 49, 3332–3346. [Google Scholar]
  42. Ding, Y.; Shi, Y.; Yang, S. Molecular Regulation of Plant Responses to Environmental Temperatures. Mol. Plant 2020, 13, 544–564. [Google Scholar] [CrossRef] [PubMed]
  43. Foyer, C.H.; Noctor, G. Redox Homeostasis and Antioxidant Signaling: A Metabolic Interface between Stress Perception and Physiological Responses. Plant Cell 2005, 17, 1866–1875. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Guan, C.; Huang, Y.H.; Cui, X.; Liu, S.J.; Zhou, Y.Z.; Zhang, Y.W. Overexpression of gene encoding the key enzyme involved in proline-biosynthesis (PuP5CS) to improve salt tolerance in switchgrass (Panicum virgatum L.). Plant Cell Rep. 2018, 37, 1187–1199. [Google Scholar] [CrossRef] [PubMed]
  45. Gill, S.S.; Tuteja, N. Reactive oxygen species and antioxidant machinery in abiotic stress tolerance in crop plants. Plant Physiol. Biochem. 2010, 48, 909–930. [Google Scholar] [CrossRef]
  46. Zhao, D.Y.; Shen, L.; Fan, B.; Liu, K.L.; Yu, M.M.; Zheng, Y.; Ding, Y.; Sheng, J.P. Physiological and genetic properties of tomato fruits from 2 cultivars differing in chilling tolerance at cold storage. J. Food Sci. 2009, 74, C348–C352. [Google Scholar] [CrossRef]
  47. Monaghan, P.; Metcalfe, N.B.; Torres, R. Oxidative stress as a mediator of life history trade-offs: Mechanisms, measurements and interpretation. Ecol. Lett. 2009, 12, 75–92. [Google Scholar] [CrossRef]
  48. Cao, S.; Cai, Y.; Yang, Z.; Zheng, Y. MeJA induces chilling tolerance in loquat fruit by regulating proline and γ-aminobutyric acid contents. Food Chem. 2012, 133, 1466–1470. [Google Scholar] [CrossRef]
  49. Hui, W.U.; Hou, L.; Zhou, Y.; Fan, Z.; Shi, J.; Aliyan, R. Analysis and Evaluation Indicator Selection of Chilling Tolerance of Different Cotton Genotypes. Agric. Sci. Technol. 2012, 13, 2338–2346. [Google Scholar]
  50. Ji, C.Y.; Kim, H.S.; Lee, C.J.; Kim, S.E.; Lee, H.U.; Nam, S.S.; Li, Q.; Ma, D.F.; Kwak, S.S. Comparative transcriptome profiling of tuberous roots of two sweetpotato lines with contrasting low temperature tolerance during storage. Gene 2020, 727, 144244. [Google Scholar] [CrossRef]
  51. Li, M.; Sui, N.; Lin, L.; Yang, Z.; Zhang, Y. Transcriptomic profiling revealed genes involved in response to cold stress in maize. Funct. Plant Biol. 2019, 46, 830–844. [Google Scholar] [CrossRef]
  52. Buti, M.; Pasquariello, M.; Ronga, D.; Milc, J.A.; Pecchioni, N.; Ho, V.T.; Pucciariello, C.; Perata, P.; Francia, E. Transcriptome profiling of short-term response to chilling stress in tolerant and sensitive Oryza sativa ssp. Japonica seedlings. Funct. Integr. Genom. 2018, 18, 627–644. [Google Scholar] [CrossRef] [PubMed]
  53. Dametto, A.; Sperotto, R.A.; Adamski, J.M.; Blasi, E.A.; Cargnelutti, D.; de Oliveira, L.F.; Ricachenevsky, F.K.; Fregonezi, J.N.; Mariath, J.E.; da Cruz, R.P.; et al. Cold tolerance in rice germinating seeds revealed by deep RNAseq analysis of contrasting indica genotypes. Plant Sci. 2015, 238, 1–12. [Google Scholar] [CrossRef] [PubMed]
  54. Pradhan, S.K.; Pandit, E.; Nayak, D.K.; Behera, L.; Mohapatra, T. Genes, pathways and transcription factors involved in seedling stage chilling stress tolerance in indica rice through RNA-Seq analysis. BMC Plant Biol. 2019, 19, 352. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Bahrman, N.; Hascoet, E.; Jaminon, O.; Depta, F.; Hu, J.F.; Bouchez, O.; Lejeune-Henaut, I.; Delbreil, B.; Legrand, S. Identification of Genes Differentially Expressed in Response to Cold in Pisum sativum Using RNA Sequencing Analyses. Plants 2019, 8, 288. [Google Scholar] [CrossRef] [Green Version]
  56. Ma, L.; Coulter, J.A.; Liu, L.; Zhao, Y.; Chang, Y.; Pu, Y.; Zeng, X.; Xu, Y.; Wu, J.; Fang, Y.; et al. Transcriptome Analysis Reveals Key Cold-Stress-Responsive Genes in Winter Rapeseed (Brassica rapa L.). Int. J. Mol. Sci. 2019, 20, 1071. [Google Scholar] [CrossRef] [Green Version]
  57. Kashash, Y.; Doron-Faigenboim, A.; Bar-Ya’akov, I.; Hatib, K.; Beja, R.; Trainin, T.; Holland, D.; Porat, R. Diversity among Pomegranate Varieties in Chilling Tolerance and Transcriptome Responses to Cold Storage. J. Agric. Food Chem. 2019, 67, 760–771. [Google Scholar] [CrossRef]
  58. Deng, S.; Ma, J.; Zhang, L.; Chen, F.; Sang, Z.; Jia, Z.; Ma, L. De novo transcriptome sequencing and gene expression profiling of Magnolia wufengensis in response to cold stress. BMC Plant Biol. 2019, 19, 321. [Google Scholar] [CrossRef] [Green Version]
  59. Eremina, M.; Rozhon, W.; Poppenberger, B. Hormonal control of cold stress responses in plants. Cell Mol. Life Sci. 2016, 73, 797–810. [Google Scholar] [CrossRef]
  60. Rahman, A. Auxin: A regulator of cold stress response. Physiol. Plant 2013, 147, 28–35. [Google Scholar] [CrossRef]
  61. Ritonga, F.N.; Chen, S. Physiological and Molecular Mechanism Involved in Cold Stress Tolerance in Plants. Plants 2020, 9, E560. [Google Scholar] [CrossRef]
  62. Kazan, K. Diverse roles of jasmonates and ethylene in abiotic stress tolerance. Trends Plant Sci. 2015, 20, 219–229. [Google Scholar] [CrossRef] [PubMed]
  63. Hu, Y.; Jiang, L.; Wang, F.; Yu, D. Jasmonate regulates the inducer of cbf expression-C-repeat binding factor/DRE binding factor1 cascade and freezing tolerance in Arabidopsis. Plant Cell 2013, 25, 2907–2924. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Kagale, S.; Divi, U.K.; Krochko, J.E.; Keller, W.A.; Krishna, P. Brassinosteroid confers tolerance in Arabidopsis thaliana and Brassica napus to a range of abiotic stresses. Planta 2007, 225, 353–364. [Google Scholar] [CrossRef] [PubMed]
  65. Divi, U.K.; Krishna, P. Brassinosteroid: A biotechnological target for enhancing crop yield and stress tolerance. New Biotechnol. 2009, 26, 131–136. [Google Scholar] [CrossRef] [PubMed]
  66. Eremina, M.; Unterholzner, S.J.; Rathnayake, A.I.; Castellanos, M.; Khan, M.; Kugler, K.G.; May, S.T.; Mayer, K.F.; Rozhon, W.; Poppenberger, B. Brassinosteroids participate in the control of basal and acquired freezing tolerance of plants. Proc. Natl. Acad. Sci. USA 2016, 113, E5982–E5991. [Google Scholar] [CrossRef] [Green Version]
  67. Xu, W.; Campbell, P.; Vargheese, A.K.; Braam, J. The Arabidopsis XET-related gene family: Environmental and hormonal regulation of expression. Plant J. 1996, 9, 879–889. [Google Scholar] [PubMed]
  68. Xu, W.; Purugganan, M.M.; Polisensky, D.H.; Antosiewicz, D.M.; Fry, S.C.; Braam, J. Arabidopsis TCH4, regulated by hormones and the environment, encodes a xyloglucan endotransglycosylase. Plant Cell 1995, 7, 1555–1567. [Google Scholar]
  69. Pu, Y.; Liu, L.; Wu, J.; Zhao, Y.; Bai, J.; Ma, L.; Yue, J.; Jin, J.; Niu, Z.; Fang, Y.; et al. Transcriptome Profile Analysis of Winter Rapeseed (Brassica napus L.) in Response to Freezing Stress, Reveal Potentially Connected Events to Freezing Stress. Int. J. Mol. Sci. 2019, 20, 2771. [Google Scholar]
  70. Jonak, C.; Kiegerl, S.; Ligterink, W.; Barker, P.J.; Huskisson, N.S.; Hirt, H. Stress signaling in plants: A mitogen-activated protein kinase pathway is activated by cold and drought. Proc. Natl. Acad. Sci. USA 1996, 93, 11274–11279. [Google Scholar] [CrossRef] [Green Version]
  71. Rodriguez, M.C.; Petersen, M.; Mundy, J. Mitogen-activated protein kinase signaling in plants. Annu. Rev. Plant Biol. 2010, 61, 621–649. [Google Scholar]
  72. Sinha, A.K.; Jaggi, M.; Raghuram, B.; Tuteja, N. Mitogen-activated protein kinase signaling in plants under abiotic stress. Plant Signal. Behav. 2011, 6, 196–203. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  73. Teige, M.; Scheikl, E.; Eulgem, T.; Doczi, R.; Ichimura, K.; Shinozaki, K.; Dangl, J.L.; Hirt, H. The MKK2 pathway mediates cold and salt stress signaling in Arabidopsis. Mol. Cell 2004, 15, 141–152. [Google Scholar] [CrossRef] [PubMed]
  74. Yang, T.; Shad Ali, G.; Yang, L.; Du, L.; Reddy, A.S.; Poovaiah, B.W. Calcium/calmodulin-regulated receptor-like kinase CRLK1 interacts with MEKK1 in plants. Plant Signal. Behav. 2010, 5, 991–994. [Google Scholar] [CrossRef] [PubMed]
  75. Xie, G.; Kato, H.; Imai, R. Biochemical identification of the OsMKK6-OsMPK3 signalling pathway for chilling stress tolerance in rice. Biochem. J. 2012, 443, 95–102. [Google Scholar] [CrossRef] [PubMed]
  76. Li, H.; Ding, Y.; Shi, Y.; Zhang, X.; Zhang, S.; Gong, Z.; Yang, S. MPK3-and MPK6-Mediated ICE1 Phosphorylation Negatively Regulates ICE1 Stability and Freezing Tolerance in Arabidopsis. Dev. Cell 2017, 43, 630–642 e4. [Google Scholar] [CrossRef] [PubMed]
  77. Cheong, Y.H.; Kim, K.N.; Pandey, G.K.; Gupta, R.; Grant, J.J.; Luan, S. CBL1, a calcium sensor that differentially regulates salt, drought, and cold responses in Arabidopsis. Plant Cell 2003, 15, 1833–1845. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  78. Kim, K.N.; Cheong, Y.H.; Grant, J.J.; Pandey, G.K.; Luan, S. CIPK3, a calcium sensor-associated protein kinase that regulates abscisic acid and cold signal transduction in Arabidopsis. Plant Cell 2003, 15, 411–423. [Google Scholar] [CrossRef] [Green Version]
  79. Wang, K.; Bai, Z.Y.; Liang, Q.Y.; Liu, Q.L.; Zhang, L.; Pan, Y.Z.; Liu, G.L.; Jiang, B.B.; Zhang, F.; Jia, Y. Transcriptome analysis of chrysanthemum (Dendranthema grandiflorum) in response to low temperature stress. BMC Genom. 2018, 19, 319. [Google Scholar] [CrossRef] [Green Version]
  80. Xiang, Y.; Huang, Y.; Xiong, L. Characterization of stress-responsive CIPK genes in rice for stress tolerance improvement. Plant Physiol. 2007, 144, 1416–1428. [Google Scholar] [CrossRef] [Green Version]
  81. Luo, Q.; Wei, Q.; Wang, R.; Zhang, Y.; Zhang, F.; He, Y.; Yang, G.; He, G. Ectopic expression of BdCIPK31 confers enhanced low-temperature tolerance in transgenic tobacco plants. Acta Biochim. Biophys. Sin. 2018, 50, 199–208. [Google Scholar] [CrossRef] [Green Version]
  82. Deng, X.; Zhou, S.; Hu, W.; Feng, J.; Zhang, F.; Chen, L.; Huang, C.; Luo, Q.; He, Y.; Yang, G.; et al. Ectopic expression of wheat TaCIPK14, encoding a calcineurin B-like protein-interacting protein kinase, confers salinity and cold tolerance in tobacco. Physiol. Plant 2013, 149, 367–377. [Google Scholar] [PubMed]
  83. Chen, J.; Xue, B.; Xia, X.; Yin, W. A novel calcium-dependent protein kinase gene from Populus euphratica, confers both drought and cold stress tolerance. Biochem. Biophys. Res. Commun. 2013, 441, 630–636. [Google Scholar] [CrossRef] [PubMed]
  84. Saijo, Y.; Hata, S.; Kyozuka, J.; Shimamoto, K.; Izui, K. Over-expression of a single Ca2+-dependent protein kinase confers both cold and salt/drought tolerance on rice plants. Plant J. 2000, 23, 319–327. [Google Scholar] [CrossRef] [PubMed]
  85. Dubrovina, A.S.; Kiselev, K.V.; Khristenko, V.S.; Aleynova, O.A. VaCPK20, a calcium-dependent protein kinase gene of wild grapevine Vitis amurensis Rupr., mediates cold and drought stress tolerance. J. Plant Physiol. 2015, 185, 1–12. [Google Scholar] [CrossRef] [PubMed]
  86. Kasuga, M.; Liu, Q.; Miura, S.; Yamaguchi-Shinozaki, K.; Shinozaki, K. Improving plant drought, salt, and freezing tolerance by gene transfer of a single stress-inducible transcription factor. Nat. Biotechnol. 1999, 17, 287–291. [Google Scholar] [CrossRef]
  87. Mizoi, J.; Shinozaki, K.; Yamaguchi-Shinozaki, K. AP2/ERF family transcription factors in plant abiotic stress responses. Biochim. Biophys. Acta 2012, 1819, 86–96. [Google Scholar] [CrossRef]
  88. Gilmour, S.J.; Zarka, D.G.; Stockinger, E.J.; Salazar, M.P.; Houghton, J.M.; Thomashow, M.F. Low temperature regulation of the Arabidopsis CBF family of AP2 transcriptional activators as an early step in cold-induced COR gene expression. Plant J. 2001, 16, 433–442. [Google Scholar] [CrossRef]
  89. Jaglo-Ottosen, K.R.; Gilmour, S.J.; Zarka, D.G.; Schabenberger, O.; Thomashow, M.F. Arabidopsis CBF1 overexpression induces COR genes and enhances freezing tolerance. Science 1998, 280, 104–106. [Google Scholar] [CrossRef] [Green Version]
  90. Novillo, F.; Medina, J.; Salinas, J. Arabidopsis CBF1 and CBF3 have a different function than CBF2 in cold acclimation and define different gene classes in the CBF regulon. Proc. Natl. Acad. Sci. USA 2008, 104, 21002–21007. [Google Scholar] [CrossRef] [Green Version]
  91. Wang, M.Q.; Huang, Q.X.; Lin, P.; Zeng, Q.H.; Li, Y.; Liu, Q.L.; Zhang, L.; Pan, Y.Z.; Jiang, B.B.; Zhang, F. The Overexpression of a Transcription Factor Gene VbWRKY32 Enhances the Cold Tolerance in Verbena bonariensis. Front. Plant Sci. 2019, 10, 1746. [Google Scholar] [CrossRef]
  92. Zhang, L.; Zhao, T.; Sun, X.; Wang, Y.; Du, C.; Zhu, Z.; Gichuki, D.K.; Wang, Q.; Li, S.; Xin, H. Overexpression of VaWRKY12, a transcription factor from Vitis amurensis with increased nuclear localization under low temperature, enhances cold tolerance of plants. Plant Mol. Biol. 2019, 100, 95–110. [Google Scholar] [CrossRef] [PubMed]
  93. Chen, Y.; Chen, Z.; Kang, J.; Kang, D.; Gu, H.; Qin, G. AtMYB14 Regulates Cold Tolerance in Arabidopsis. Plant Mol. Biol. Rep. 2013, 31, 87–97. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  94. An, J.P.; Li, R.; Qu, F.J.; You, C.X.; Wang, X.F.; Hao, Y.J. R2R3-MYB transcription factor MdMYB23 is involved in the cold tolerance and proanthocyanidin accumulation in apple. Plant J. 2018, 96, 562–577. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  95. Shi, H.; Chan, Z. The cysteine2/histidine2-type transcription factor ZINC FINGER OF ARABIDOPSIS THALIANA 6-activated C-REPEAT-BINDING FACTOR pathway is essential for melatonin-mediated freezing stress resistance in Arabidopsis. J. Pineal Res. 2014, 57, 185–191. [Google Scholar] [CrossRef] [PubMed]
  96. Pertea, M.; Kim, D.; Pertea, G.M.; Leek, J.T.; Salzberg, S.L. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 2016, 11, 1650–1667. [Google Scholar] [CrossRef] [PubMed]
  97. Bastian, M.; Heymann, S.; Jacomy, M. Gephi: An Open Source Software for Exploring and Manipulating Networks. In Proceedings of the 3rd International Conference on Weblogs and Social Media, San Jose, CA, USA; 2009. [Google Scholar]
  98. Lu, P.; Magwanga, R.O.; Kirungu, J.N.; Hu, Y.; Dong, Q.; Cai, X.; Zhou, Z.; Wang, X.; Zhang, Z.; Hou, Y.; et al. Overexpression of Cotton a DTX/MATE Gene Enhances Drought, Salt, and Cold Stress Tolerance in Transgenic Arabidopsis. Front. Plant Sci. 2019, 10, 299. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Assessment of seedling chilling injury (CI) index and investigation of survival rate. (A) Evaluation of the cotyledon CI index of 74 upland cotton varieties. The seedlings growing to the seventh day were treated with cold stress (4 °C) for 36 h and then recovered for another week, after which the chilling phenotype was investigated. Error bar indicates the standard error (SE) of the mean (n = 3). (B) Chilling injury symptoms of H559 and YM21 seedlings recovered one week after 4 °C stress. (C) The seedling survival rate of H559 and YM21 after 48 h of chilling stress.
Figure 1. Assessment of seedling chilling injury (CI) index and investigation of survival rate. (A) Evaluation of the cotyledon CI index of 74 upland cotton varieties. The seedlings growing to the seventh day were treated with cold stress (4 °C) for 36 h and then recovered for another week, after which the chilling phenotype was investigated. Error bar indicates the standard error (SE) of the mean (n = 3). (B) Chilling injury symptoms of H559 and YM21 seedlings recovered one week after 4 °C stress. (C) The seedling survival rate of H559 and YM21 after 48 h of chilling stress.
Ijms 21 05095 g001
Figure 2. Physiological changes in young cotyledons of H559 and YM21 under chilling stress (4 °C) at different time points. These physiological indices include superoxide dismutase (SOD) activity, peroxidase (POD) activity, catalase (CAT) activity, malondialdehyde (MDA) content, soluble sugar (SS) content, soluble protein (SP) content, O2− content, H2O2 content, proline content, and superoxide anion scavenging (SAC) rate. Error bar indicates the standard error (SE) of the mean (n = 3). The asterisk indicates a significant difference level between H559 and YM21 (** p < 0.01).
Figure 2. Physiological changes in young cotyledons of H559 and YM21 under chilling stress (4 °C) at different time points. These physiological indices include superoxide dismutase (SOD) activity, peroxidase (POD) activity, catalase (CAT) activity, malondialdehyde (MDA) content, soluble sugar (SS) content, soluble protein (SP) content, O2− content, H2O2 content, proline content, and superoxide anion scavenging (SAC) rate. Error bar indicates the standard error (SE) of the mean (n = 3). The asterisk indicates a significant difference level between H559 and YM21 (** p < 0.01).
Ijms 21 05095 g002
Figure 3. Venn diagram analysis of expressed genes in all samples and analysis of constitutive transcriptional differences between H559 and YM21. (A) The number of expressed genes identified in different samples. The horizontal bar shows the total number of expressed genes identified in each sample. The colored dots connected by lines indicate the sharing of genes between samples, and the number of genes corresponds to the vertical bar above. (B) The upregulated or downregulated genes in S0 compared to T0. The colored dots represent genes, where the red dots represent upregulated genes, the green dots represent downregulated genes in T0, and the khaki dots only represent expressed genes. (C) Venn diagram of expressed and differentially expressed genes within and between varieties before cold stress.
Figure 3. Venn diagram analysis of expressed genes in all samples and analysis of constitutive transcriptional differences between H559 and YM21. (A) The number of expressed genes identified in different samples. The horizontal bar shows the total number of expressed genes identified in each sample. The colored dots connected by lines indicate the sharing of genes between samples, and the number of genes corresponds to the vertical bar above. (B) The upregulated or downregulated genes in S0 compared to T0. The colored dots represent genes, where the red dots represent upregulated genes, the green dots represent downregulated genes in T0, and the khaki dots only represent expressed genes. (C) Venn diagram of expressed and differentially expressed genes within and between varieties before cold stress.
Ijms 21 05095 g003
Figure 4. Nonredundant GO-terms related to biological processes obtained by enrichment analysis of the upregulated genes in H559 (A) and YM21 (B) before cold stress. The REVIGO program was used to de-redundant the significant GO-terms, and bubble charts were then drawn by using the R package to display the nonredundant GO-terms. Nonredundant GO-terms were clustered on a two-dimensional space according to the semantic similarity between them (adjacent bubbles are most closely connected). Bubble size is proportional to the uniqueness value of each GO-term, and the color represents the log10 (p-value).
Figure 4. Nonredundant GO-terms related to biological processes obtained by enrichment analysis of the upregulated genes in H559 (A) and YM21 (B) before cold stress. The REVIGO program was used to de-redundant the significant GO-terms, and bubble charts were then drawn by using the R package to display the nonredundant GO-terms. Nonredundant GO-terms were clustered on a two-dimensional space according to the semantic similarity between them (adjacent bubbles are most closely connected). Bubble size is proportional to the uniqueness value of each GO-term, and the color represents the log10 (p-value).
Ijms 21 05095 g004
Figure 5. Weighted gene coexpression network analysis (WGCNA) of 19,982 DEGs obtained from all pairwise comparisons. (A) Hierarchical clustering tree showing 11 coexpression modules identified by WGCNA. Different modules are marked with different colors. Each leaf of the cluster tree represents a gene. (B) The number of transcription factors and nontranscription factor genes for each module. The orange bar indicates transcription factors, and the green bar indicates nontranscription factor genes. (C) Module–sample correlations and corresponding p-values. Each row represents a specific module, and each column represents a sample. The heatmap on the right shows the Pearson correlation between module eigengenes and samples. The numbers in each cell represent the correlation coefficients and correlation significance levels (in parentheses). The color of the cell reflects the degree of correlation. (D) Module–trait correlations and corresponding p-values. These traits correspond to the eight physiological indexes mentioned above.
Figure 5. Weighted gene coexpression network analysis (WGCNA) of 19,982 DEGs obtained from all pairwise comparisons. (A) Hierarchical clustering tree showing 11 coexpression modules identified by WGCNA. Different modules are marked with different colors. Each leaf of the cluster tree represents a gene. (B) The number of transcription factors and nontranscription factor genes for each module. The orange bar indicates transcription factors, and the green bar indicates nontranscription factor genes. (C) Module–sample correlations and corresponding p-values. Each row represents a specific module, and each column represents a sample. The heatmap on the right shows the Pearson correlation between module eigengenes and samples. The numbers in each cell represent the correlation coefficients and correlation significance levels (in parentheses). The color of the cell reflects the degree of correlation. (D) Module–trait correlations and corresponding p-values. These traits correspond to the eight physiological indexes mentioned above.
Ijms 21 05095 g005
Figure 6. GO enrichment analysis of four significant modules. All significant GO-terms presented in the heatmap have been de-redundant with the REVIGO online platform. Each column represents a significant module, and each row represents a GO-term that is significantly (p-value < 0.01) enriched in at least one module. The color scale on the left represents specific intervals of significance level (corrected p-value).
Figure 6. GO enrichment analysis of four significant modules. All significant GO-terms presented in the heatmap have been de-redundant with the REVIGO online platform. Each column represents a significant module, and each row represents a GO-term that is significantly (p-value < 0.01) enriched in at least one module. The color scale on the left represents specific intervals of significance level (corrected p-value).
Ijms 21 05095 g006
Figure 7. Coexpression network construction of the blue module. (A) Gene coexpression heatmap of the blue module (upper panel) and expression level of the corresponding eigengene in each sample (lower panel). (B) Dot plot mining hub genes with higher module connectivity and expression levels at T24. The X-axis represents the connectivity values (Kme) of genes in the blue module, while the Y-axis represents the dominant expression value (SSCT24) of the genes at T24. Genes with Kme > 0.8 and SSCT24 > 0.5 were identified as hub genes (blue dots). (C) The coexpression network of hub genes in the blue module. The top 100 hub genes with weight value > 0.45 and higher expression levels at T24 were used to construct the network. The node size indicated gene degree, and the color was used to distinguish whether it was a transcription factor gene (orange node). Line thickness between nodes reflects the weight value.
Figure 7. Coexpression network construction of the blue module. (A) Gene coexpression heatmap of the blue module (upper panel) and expression level of the corresponding eigengene in each sample (lower panel). (B) Dot plot mining hub genes with higher module connectivity and expression levels at T24. The X-axis represents the connectivity values (Kme) of genes in the blue module, while the Y-axis represents the dominant expression value (SSCT24) of the genes at T24. Genes with Kme > 0.8 and SSCT24 > 0.5 were identified as hub genes (blue dots). (C) The coexpression network of hub genes in the blue module. The top 100 hub genes with weight value > 0.45 and higher expression levels at T24 were used to construct the network. The node size indicated gene degree, and the color was used to distinguish whether it was a transcription factor gene (orange node). Line thickness between nodes reflects the weight value.
Ijms 21 05095 g007
Figure 8. Coexpression network construction of the saddlebrown module. (A) Gene coexpression heatmap of the saddlebrown module (upper panel) and expression level of the corresponding eigengene in each sample (lower panel). (B) Dot plot mining hub genes with higher module connectivity and expression levels at T3. The X-axis represents the connectivity values (Kme) of genes in the saddlebrown module, while the Y-axis represents the dominant expression value (SSCT3) of the genes at T3. Genes with Kme > 0.8 and SSCT3 > 0.50 were identified as hub genes (saddlebrown dots). (C) The coexpression network of all coexpressed genes in the saddlebrown module. The node size indicated gene degree. Line thickness between nodes reflects the weight value.
Figure 8. Coexpression network construction of the saddlebrown module. (A) Gene coexpression heatmap of the saddlebrown module (upper panel) and expression level of the corresponding eigengene in each sample (lower panel). (B) Dot plot mining hub genes with higher module connectivity and expression levels at T3. The X-axis represents the connectivity values (Kme) of genes in the saddlebrown module, while the Y-axis represents the dominant expression value (SSCT3) of the genes at T3. Genes with Kme > 0.8 and SSCT3 > 0.50 were identified as hub genes (saddlebrown dots). (C) The coexpression network of all coexpressed genes in the saddlebrown module. The node size indicated gene degree. Line thickness between nodes reflects the weight value.
Ijms 21 05095 g008
Figure 9. The proportion of genes in the top 10 abundant transcription factor (TF) families in blue (left) and green (right) modules.
Figure 9. The proportion of genes in the top 10 abundant transcription factor (TF) families in blue (left) and green (right) modules.
Ijms 21 05095 g009
Figure 10. The DEGs of four significant modules in the plant hormone signaling pathway. (A) Distribution of genes from four modules on eight plant hormone pathways. The pink boxes denoted the genes enriched, and the colored dots below or to the right of the pink boxes indicated from which module the genes came. The gray boxes denoted the genes that had not been enriched. The meanings of different types of arrows have been strictly defined in KEGG databases (https://www.genome.jp/kegg/document/help_pathway.html). (B) The expression patterns of genes distributed on the eight hormone signaling pathways in the four modules. Gene expression data were normalized by log10 (FPKM + 0.01). The colored bar on the left of each heatmap represented the modules from which the corresponding genes on the right come from.
Figure 10. The DEGs of four significant modules in the plant hormone signaling pathway. (A) Distribution of genes from four modules on eight plant hormone pathways. The pink boxes denoted the genes enriched, and the colored dots below or to the right of the pink boxes indicated from which module the genes came. The gray boxes denoted the genes that had not been enriched. The meanings of different types of arrows have been strictly defined in KEGG databases (https://www.genome.jp/kegg/document/help_pathway.html). (B) The expression patterns of genes distributed on the eight hormone signaling pathways in the four modules. Gene expression data were normalized by log10 (FPKM + 0.01). The colored bar on the left of each heatmap represented the modules from which the corresponding genes on the right come from.
Ijms 21 05095 g010
Figure 11. Validation of the expression patterns of 20 genes in four significant WGCNA modules by using qRT-PCR. (A) The qRT-PCR results of 20 genes. Here we additionally showed the expression level of each gene at 6 h and 12 h, and only used the qRT-PCR results at 0 h, 3 h, and 6 h under cold stress to perform correlation analysis with RNA-seq data. The correlation coefficient (Cor) and correlation significance level (p-value) of the expression pattern was added above each histogram. (B) Validation of gene expression patterns in cold-sensitive variety YM21. Each point represents a value of fold change of expression level at S3 or S24 compared with that at S0 or S3. (C) Validation of gene expression patterns in cold-tolerant variety H559. Each point represents a value of fold change of expression level at T3 or T24 compared with that at T0 or T3. The fold change values were converted by log2 standardization.
Figure 11. Validation of the expression patterns of 20 genes in four significant WGCNA modules by using qRT-PCR. (A) The qRT-PCR results of 20 genes. Here we additionally showed the expression level of each gene at 6 h and 12 h, and only used the qRT-PCR results at 0 h, 3 h, and 6 h under cold stress to perform correlation analysis with RNA-seq data. The correlation coefficient (Cor) and correlation significance level (p-value) of the expression pattern was added above each histogram. (B) Validation of gene expression patterns in cold-sensitive variety YM21. Each point represents a value of fold change of expression level at S3 or S24 compared with that at S0 or S3. (C) Validation of gene expression patterns in cold-tolerant variety H559. Each point represents a value of fold change of expression level at T3 or T24 compared with that at T0 or T3. The fold change values were converted by log2 standardization.
Ijms 21 05095 g011
Table 1. Summary of sequencing data for different samples.
Table 1. Summary of sequencing data for different samples.
SampleRaw Data (bp)Clean Data (bp)Q30 (%)GC (%)Unmapped (%)Unique_Mapped (%)Total_Mapped (%)
T0-17.66 × 109 7.62 × 10991.20%46.09%3.42%91.90%96.58%
T0-27.20 × 1097.17 × 10992.04%46.10%3.24%92.05%96.76%
T0-37.95 × 1097.92 × 10991.27%45.98%3.61%91.62%96.39%
T3-17.59 × 1097.55 × 10991.08%46.40%3.83%90.17%96.17%
T3-27.56 × 1097.52 × 10992.15%46.15%3.19%91.04%96.81%
T3-37.04 × 1097.01 × 10991.68%46.47%3.08%91.00%96.92%
T24-18.19 × 1098.16 × 10992.32%46.28%3.41%91.45%96.59%
T24-27.49 × 1097.46 × 10991.79%46.07%3.08%91.97%96.92%
T24-38.69 × 1098.65 × 10991.75%46.26%3.39%91.31%96.61%
S0-16.53 × 1096.50 × 10992.02%46.17%3.18%92.14%96.82%
S0-27.95 × 1097.92 × 10991.17%46.00%3.51%91.84%96.49%
S0-38.07 × 1098.04 × 10991.45%45.83%3.38%92.01%96.62%
S3-17.05 × 1097.02 × 10992.09%46.78%2.96%91.63%97.04%
S3-28.39 × 1098.36 × 10991.23%46.03%3.90%90.84%96.10%
S3-38.98 × 1098.95 × 10992.33%46.25%2.53%91.83%97.47%
S24-16.89 × 1096.87 × 10992.16%45.93%2.56%92.27%97.44%
S24-28.77 × 1098.73 × 10991.36%46.16%3.39%91.13%96.61%
S24-38.60 × 1098.57 × 10991.55%46.08%7.48%87.60%92.52%
Table 2. The number of differentially expressed genes within and between H559 and YM21 during chilling stress.
Table 2. The number of differentially expressed genes within and between H559 and YM21 during chilling stress.
ComparisonTotalUpregulatedDownregulated
T0-vs-T31626954672
T0-vs-T2413,55845918967
T3-vs-T2411,90242247678
S0-vs-S320301056972
S0-vs-S2411,79545227273
S3-vs-S24898441524832
S0-vs-T01504964540
S3-vs-T327271886841
S24-vs-T2423161397919
Table 3. Significant Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of four significant WGCNA modules.
Table 3. Significant Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of four significant WGCNA modules.
ModulePathway IDPathwayNumber of Genesp-Value
Blueko04626Plant–pathogen interaction1141.60 × 10−14
ko04075Plant hormone signal transduction851.15 × 10−9
ko04016MAPK signaling pathway562.58 × 10−9
ko00062Fatty acid elongation121.66 × 10−3
ko00750Vitamin B6 metabolism72.65 × 10−3
ko00601Glycosphingolipid biosynthesis35.61 × 10−3
ko00330Arginine and proline metabolism159.27 × 10−3
Greenko04016MAPK signaling pathway251.35 × 10−8
ko04075Plant hormone signal transduction343.49 × 10−8
ko00730Thiamine metabolism61.12 × 10−4
ko04712Circadian rhythm72.02 × 10−3
Paleturquoiseko00071Fatty acid degradation87.08 × 10−6
ko00592alpha-Linolenic acid metabolism73.46 × 10−5
ko00280Valine, leucine, and isoleucine degradation51.03 × 10−3
ko00590Arachidonic acid metabolism32.58 × 10−3
ko00630Glyoxylate and dicarboxylate metabolism62.70 × 10−3
ko04075Plant hormone signal transduction142.74 × 10−3
ko00511Other glycan degradation33.29 × 10−3
ko00565Ether lipid metabolism43.86 × 10−3
ko01100Metabolic pathways509.82 × 10−3
ko00260Glycine, serine, and threonine metabolism59.92 × 10−3
Saddlebrownko00940Phenylpropanoid biosynthesis63.25 × 10−4
ko00480Glutathione metabolism42.02 × 10−3
ko00250Alanine, aspartate, and glutamate metabolism32.79 × 10−3
Table 4. The top five hub genes with high Z-score value in the dominant expression period of each module.
Table 4. The top five hub genes with high Z-score value in the dominant expression period of each module.
Module (Preponderant Stage)Gene IDSymbolDescriptionZ-ScoreKmeSSC
Blue (T24)GH_D12G2915CML27probable calcium-binding protein CML2725.890.990.54
GH_D09G1466NHL10protein YLS9-like17.860.990.55
GH_A10G2658WUN1wound-induced protein 1-like14.420.980.64
GH_D05G2128ZAT10zinc finger protein ZAT10-like13.540.980.52
GH_A08G1036TET7tetraspanin-8-like12.910.990.53
Green (S24)GH_A13G1006pcbABN-(5-amino-5-carboxypentanoyl)-L-cysteinyl-D-valine synthase2.750.990.58
GH_A08G0687AUX22Dauxin-induced protein 22D-like2.410.980.67
GH_D12G1406NIMIN-1protein NIM1-INTERACTING 1-like2.390.980.50
GH_A05G4197SIB2sigma factor binding protein 2, chloroplastic-like2.130.990.53
GH_D02G1771At1g60420probable nucleoredoxin 11.910.900.56
Paleturquoise (S3)GH_D06G1504DRM1auxin-repressed 12.5 kDa protein isoform X228.360.970.52
GH_A10G0258TSJT1stem-specific protein TSJT1-like11.380.920.55
GH_A07G0097OXS3Oxidative stress 310.420.940.58
GH_A07G0096PDD4Lprogrammed cell death protein 4-like5.360.940.55
GH_A11G3691TIP1-1aquaporin TIP1-1-like2.770.970.55
Saddlebrown (T3)GH_D04G1587BETV1Fpathogenesis-related protein STH-2-like1.040.910.67
GH_A12G2726CHIT3hevamine-A-like0.900.940.57
GH_A02G0099-uncharacterized protein LOC108475098 [Gossypium arboreum]0.510.900.67
GH_D02G0103-uncharacterized protein LOC105796351 [Gossypium raimondii]0.190.930.69
GH_D06G1735GPT2glucose-6-phosphate/phosphate translocator 2, chloroplastic-like0.120.900.67

Share and Cite

MDPI and ACS Style

Cheng, G.; Zhang, L.; Wang, H.; Lu, J.; Wei, H.; Yu, S. Transcriptomic Profiling of Young Cotyledons Response to Chilling Stress in Two Contrasting Cotton (Gossypium hirsutum L.) Genotypes at the Seedling Stage. Int. J. Mol. Sci. 2020, 21, 5095. https://doi.org/10.3390/ijms21145095

AMA Style

Cheng G, Zhang L, Wang H, Lu J, Wei H, Yu S. Transcriptomic Profiling of Young Cotyledons Response to Chilling Stress in Two Contrasting Cotton (Gossypium hirsutum L.) Genotypes at the Seedling Stage. International Journal of Molecular Sciences. 2020; 21(14):5095. https://doi.org/10.3390/ijms21145095

Chicago/Turabian Style

Cheng, Gongmin, Longyan Zhang, Hantao Wang, Jianhua Lu, Hengling Wei, and Shuxun Yu. 2020. "Transcriptomic Profiling of Young Cotyledons Response to Chilling Stress in Two Contrasting Cotton (Gossypium hirsutum L.) Genotypes at the Seedling Stage" International Journal of Molecular Sciences 21, no. 14: 5095. https://doi.org/10.3390/ijms21145095

APA Style

Cheng, G., Zhang, L., Wang, H., Lu, J., Wei, H., & Yu, S. (2020). Transcriptomic Profiling of Young Cotyledons Response to Chilling Stress in Two Contrasting Cotton (Gossypium hirsutum L.) Genotypes at the Seedling Stage. International Journal of Molecular Sciences, 21(14), 5095. https://doi.org/10.3390/ijms21145095

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop