Comparative Transcriptome Analysis Provides Insights into the Seed Germination in Cotton in Response to Chilling Stress

Gossypium hirsutum L., is a widely cultivated cotton species around the world, but its production is seriously threatened by its susceptibility to chilling stress. Low temperature affects its germination, and the underlying molecular mechanisms are rarely known, particularly from a transcriptional perspective. In this study, transcriptomic profiles were analyzed and compared between two cotton varieties, the cold-tolerant variety KN27-3 and susceptible variety XLZ38. A total of 7535 differentially expressed genes (DEGs) were identified. Among them, the transcripts involved in energy metabolism were significantly enriched during germination based on analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, such as glycolysis/gluconeogenesis, tricarboxylic acid cycle (TCA cycle), and glyoxylate cycle (GAC). Results from further GO enrichment analysis show the earlier appearance of DNA integration, meristem growth, cotyledon morphogenesis, and other biological processes in KN27-3 compared with XLZ38 under chilling conditions. The synthesis of asparagine, GDP-mannose, and trehalose and the catabolic process of raffinose were activated. DEGs encoding antioxidants (spermidine) and antioxidase (CAT1, GPX4, DHAR2, and APX1) were much more up-regulated in embryos of KN27-3. The content of auxin (IAA), cis-zeatin riboside (cZR), and trans-zeatin riboside (tZR) in KN27-3 are higher than that in XLZ38 at five stages (from 12 h to 54 h). GA3 was expressed at a higher level in KN27-3 from 18 h to 54 h post imbibition compared to that in XLZ38. And abscisic acid (ABA) content of KN27-3 is lower than that in XLZ38 at five stages. Results from hormone content measurements and the related gene expression analysis indicated that IAA, CTK, and GA3 may promote germination of the cold-tolerant variety, while ABA inhibits it. These results expand the understanding of cottonseed germination and physiological regulations under chilling conditions by multiple pathways.


Introduction
Cotton (Gossypium hirsutum L.; AADD, 2n = 4× = 52) is the world's most important crop for natural fiber, but it is sensitive to chilling stress [1]. With the increasing needs of cotton, the cotton planting regions have expanded into the areas of high latitude and altitude with unsuitable climates [2]. Planting normally starts at spring when the night temperatures are well below 12 • C for most cotton-growing regions, including southeastern Turkey [3] and northwestern China [4]. Low temperature leads to delayed growth and development, especially at the germination and seedling stages [5]. In the meantime, more frequent chilling stress further affects cotton production [6].
Plants respond to low temperature with a series of physiological and biochemical reactions controlled by a finely-tuned regulatory network [7]. Through the network, genes encoding protective proteins such as osmotin, mRNA binding proteins, and key enzymes for biosynthesis of metabolites are activated for osmotic adjustment [8]. The accumulation of metabolites like glycine betaine, proline, soluble sugars, and polyols can adjust osmotic homeostasis against chilling stress [9,10]. Furthermore, genes encoding stress responsive proteins, including various transcription factors, proteins kinases, and other signaling molecules, are also activated [11]. C-REPEAT BINDING FACTOR (CBFs)-controlled gene regulation mode is the most well-known mechanism associated with cold stress, but it only controls a share of the cold-responsive transcriptome. There are other CBF-independent regulatory mechanisms in plants. For example, ZAT10 is a negative regulator of COR genes in which activity is repressed by the enolase LOW EXPRESSION OF OSMOTICALLY RESPONSIVE GENES 2/C-MYC BINDING PROTEIN (LOS2/AtMBP-1) under cold conditions [12].
Hormone changes also play an important role in response to chilling stress. Auxin (IAA) is required for plant root development. Exogenous IAA can promote root development [13]. A high level of IAA can maintain normal starch hydrolysis and sugar consumption to promote seed germination under low temperature [14]. Bio-active gibberellins (GAs) increase after seed imbibition, and this increase occurs just before radicle protrusion [15]. It has been reported that the reduction of GA accumulation along with the increase of abscisic acid (ABA) in seeds inhibits seed germination under chilling injury [16]. In addition, low temperature treatment decreases the content of cytokinin (CTK) and CTK oxidase in wheat [17]. Thus, better cold tolerance may be obtained by altering hormone levels.
Seed germination is considered as the most critical developmental phase in the life cycle of seed plants. It requires coordination of complex cellular and biochemical processes, including cell cycle activation, hormonal balance, energy metabolism, and repair responses [18]. According to the water absorption rate, seed germination is divided into three phases. Phase I is characterized with rapid water absorption, and it is the initial stage of germination. At phase I, the cell membrane turns from gel phase to liquid crystal status, and the respiration of seeds is rapidly enhanced [19]. Phase II is named as the absorbing platform. It is the critical phase of germination because gene expression and metabolic activity are largely enhanced [20,21]. Phase III is indicated with the appearance of the hypocotyl and continued water absorption and seedling growth. It represents the post-germination period [22]. Normal seed germination and seedling development assure a good beginning of subsequent cultivation and harvest [23]. Therefore, it is important to understand the regulatory mechanisms of seed germination.
Seed germination has been studied by using transcriptomic and proteomic analysis, either independently or in combination, in legumes [24,25], Arabidopsis thaliana [26], cereals [27], and other species [28]. In our study, two varieties of cotton exhibiting a significant difference in chilling response during germination were selected. The hypocotyl of KN27-3 (cold-tolerant variety) grew faster compared to XLZ38 (cold-susceptible variety) at 16 • C /4 • C. Their embryos from five hours post imbibition (hpi) of germination were harvested for transcriptome analysis and endogenous hormones test. The results characterized the major features of the chilling response during imbibition and identified low temperature related genes that can be used for molecular breeding.

Morphological Characteristics of KN27-3 and XLZ38 during Germination under Chilling Stress
Twenty-one different varieties of cotton were investigated for their hypocotyl length during germination under chilling conditions (Supplementary File S1). The biggest difference in seed germination speed was found between KN27-3 and XLZ38 (Figure 1a). Under normal conditions (28 • C/25 • C), there was no difference in either 100-grain weight (Figure 1b) or hypocotyl length (Figure 1c). Under chilling conditions (16 • C /4 • C), the hypocotyl length of KN27-3 was 0.94 cm, which was about 30% longer than that of XLZ38 on the 6 th day after imbibition. The length difference got larger on the 8 th day and 10 th day. On the 12 th day, the hypocotyl length of KN27-3 was about 45% longer than that of XLZ38 (Figure 1d), indicating that the KN27-3 variety was more cold-tolerant.

Morphological Characteristics of KN27-3 and XLZ38 during Germination under Chilling Stress
Twenty-one different varieties of cotton were investigated for their hypocotyl length during germination under chilling conditions (Supplementary File S1). The biggest difference in seed germination speed was found between KN27-3 and XLZ38 (Figure 1a). Under normal conditions (28 ℃/25 ℃), there was no difference in either 100-grain weight (Figure 1b) or hypocotyl length ( Figure   1c). Under chilling conditions (16 °C /4 °C), the hypocotyl length of KN27-3 was 0.94 cm, which was about 30% longer than that of XLZ38 on the 6 th day after imbibition. The length difference got larger on the 8 th day and 10 th day. On the 12 th day, the hypocotyl length of KN27-3 was about 45% longer than that of XLZ38 (Figure 1d), indicating that the KN27-3 variety was more cold-tolerant. Bars marked with stars "*" are the samples showing difference (p < 0.05) according to Duncan's multiple range test using SPSS software. "NS" represents no difference.

RNA-Seq Analysis of the Two Varieties at the Early Germination Stage under Chilling Conditions and Identification of Differentially Expressed Genes (DEGs)
Under the 16 ℃/4 ℃ conditions, cottonseeds began to absorb moisture at a faster rate for the first 42 h than 54 h, belonging to Phase I. From 42 h to 54 h, seeds barely absorbed water, belonging to Phase II (Supplementary Figure S1). In order to understand the reasons for the different chilling resistance in the two cotton varieties, samples at five imbibition time points were harvested for RNA extraction and RNA-Seq analysis. At each time point, three biological replicates were designed in both two varieties, and embryos at five imbibition time points were collected. A total of thirty samples were obtained.
Approximately 1.41 × 10 9 raw reads were generated for the thirty samples via RNA-Seq whereas, 98.28% clean reads were obtained (Supplementary File S2). By anchoring the clean reads to the cotton reference genome, we obtained the expression data of 76,331 genes (Supplementary File S3). The Arabidopsis thaliana database was then used to determine the functional annotations of cotton genes using the BLAST algorithm. To further identify the genes closely correlated with the phenotypes of the two varieties, differentially expressed genes (DEGs) were further obtained by Bars marked with asterisks indicate differences change between the same hour (* p < 0.05) according to Duncan's multiple range test using SPSS software. "NS" represents no difference.

RNA-Seq Analysis of the Two Varieties at the Early Germination Stage under Chilling Conditions and Identification of Differentially Expressed Genes (DEGs)
Under the 16 • C/4 • C conditions, cottonseeds began to absorb moisture at a faster rate for the first 42 h than 54 h, belonging to Phase I. From 42 h to 54 h, seeds barely absorbed water, belonging to Phase II (Supplementary Figure S1). In order to understand the reasons for the different chilling resistance in the two cotton varieties, samples at five imbibition time points were harvested for RNA extraction and RNA-Seq analysis. At each time point, three biological replicates were designed in both two varieties, and embryos at five imbibition time points were collected. A total of thirty samples were obtained. Approximately 1.41 × 10 9 raw reads were generated for the thirty samples via RNA-Seq whereas, 98.28% clean reads were obtained (Supplementary File S2). By anchoring the clean reads to the cotton reference genome, we obtained the expression data of 76,331 genes (Supplementary File S3). The Arabidopsis thaliana database was then used to determine the functional annotations of cotton genes using the BLAST algorithm. To further identify the genes closely correlated with the phenotypes of the two varieties, differentially expressed genes (DEGs) were further obtained by comparing their expression levels in KN27-3 and XLZ38 at five imbibition stages. The results show 1631, 1597, 1510, 1161, and 933 genes significantly up-regulated XLZ38 relative to KN27-3 at stages one to five, while 1318, 1149, 1191, 685, and 994 genes were down-regulated. There were a total of 7535 genes expressed differentially ( Figure 2). There were more up-regulated genes than down-regulated ones from stage one to stage four.  Figure 2). There were more up-regulated genes than down-regulated ones from stage one to stage four.

KEGG Pathway Analysis of DEGs
To further elucidate the functions of the DEGs in response to chilling stress, Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to classify DEGs into the corresponding pathways. Twenty KEGG pathways were significantly enriched in five comparisons ( Figure 3a, Supplementary File S4). Carbon metabolism was an essential process in plants to produce both energy sources and structural components of cells. In our study, carbon metabolism was the most enriched pathway containing 194 DEGs involved in "glyoxylate and dicarboxylate metabolism" (47 unigenes), "glycolysis/gluconeogenesis" (71 unigenes), "citrate cycle" (38 unigenes), and "galactose metabolism" (38 unigenes). As shown in figure 3b, the expression changes of most DEGs involved in glycolysis/gluconeogenesis exhibited the same dynamic trend in both varieties, but the increase was significantly stronger in KN27-3 at 42 and 54 h than in XLZ38 (Figure 3b, Supplementary File S5). These genes encode fructose-bisphosphate aldolase (FBA), phosphoglycerate kinase (PGK), and phosphoglycerate mutase (PGM). In addition, the genes encoding alcohol dehydrogenase (ADH) with a function in ethanol production under anaerobic conditions was significantly more up-regulated at 12, 18, and 30 hpi in XLZ38.

KEGG Pathway Analysis of DEGs
To further elucidate the functions of the DEGs in response to chilling stress, Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to classify DEGs into the corresponding pathways. Twenty KEGG pathways were significantly enriched in five comparisons (Figure 3a, Supplementary File S4). Carbon metabolism was an essential process in plants to produce both energy sources and structural components of cells. In our study, carbon metabolism was the most enriched pathway containing 194 DEGs involved in "glyoxylate and dicarboxylate metabolism" (47 unigenes), "glycolysis/gluconeogenesis" (71 unigenes), "citrate cycle" (38 unigenes), and "galactose metabolism" (38 unigenes). As shown in Figure 3b, the expression changes of most DEGs involved in glycolysis/gluconeogenesis exhibited the same dynamic trend in both varieties, but the increase was significantly stronger in KN27-3 at 42 and 54 h than in XLZ38 ( (from 42 to 54 hpi), the energy came from glycolysis and the TCA cycle for hypocotyl elongation in KN27-3.

GO Function Enrichment Analyses of DEGs
DEGs obtained from the five comparisons were used to search against the Gene Ontology (GO) database and the analysis of GO categories was performed. The results showed that the biological processes (BP) related to seed germination were significantly enriched earlier in KN27-3, indicating that KN27-3 developed faster than XLZ38. Figure 4 shows GO terms with up-regulated genes enrichment in KN27-3, which were not present in the down-regulated genes enrichment results (Supplementary File S7). After the seeds absorbed water for 12 hours, DNA integration and RNA assembly were enriched. Anaerobic respiration and glyoxylate metabolism were enriched in the embryo. At 18 hpi, DNA was repaired, vacuole formed, the seed switched from anaerobic Relative levels of DEG expression are shown by a color gradient from low (blue) to high (red); (e) fragments per kilobase of transcripts per million mapped fragments (FPKM) of different major genes involved in the above metabolic pathways at different hours post imbibition. K, KN27-3. X, XLZ38. The numbers along the X-axis represent the sampling hours post imbibitions. The numbers in the scale bars indicate the log2 (fold changes) in gene expression. Bars marked with asterisks indicate differences change between the same hour (* p < 0.05, ** p < 0.01) according to Duncan's multiple range test using SPSS software. "NS" represents no difference.

GO Function Enrichment Analyses of DEGs
DEGs obtained from the five comparisons were used to search against the Gene Ontology (GO) database and the analysis of GO categories was performed. The results showed that the biological processes (BP) related to seed germination were significantly enriched earlier in KN27-3, indicating that KN27-3 developed faster than XLZ38. Figure 4 shows GO terms with up-regulated genes enrichment in KN27-3, which were not present in the down-regulated genes enrichment results (Supplementary File S7). After the seeds absorbed water for 12 hours, DNA integration and RNA assembly were enriched. Anaerobic respiration and glyoxylate metabolism were enriched in the embryo. At 18 hpi, DNA was repaired, vacuole formed, the seed switched from anaerobic respiration to aerobic, meristematic cells started differentiation to form new organs, DNA began to replicate, and mitosis was activated at 30 hpi. In addition, the synthesis of chlorophyll and photomorphogenesis were also significantly enriched. At 42 h, hypocotyls were developing, anther and gynoecium development emerged in KN27-3. At the last 54 hpi, meristem continued to grow and the cotyledon began to form in KN27-3 ( Figure 4). XLZ38 had a similar growth pattern throughout but was delayed compared to  respiration to aerobic, meristematic cells started differentiation to form new organs, DNA began to replicate, and mitosis was activated at 30 hpi. In addition, the synthesis of chlorophyll and photomorphogenesis were also significantly enriched. At 42 h, hypocotyls were developing, anther and gynoecium development emerged in KN27-3. At the last 54 hpi, meristem continued to grow and the cotyledon began to form in KN27-3 ( Figure 4). XLZ38 had a similar growth pattern throughout but was delayed compared to KN27-3.

DEGs Related to Osmoregulation
After the onset of chilling, plants often increase the accumulation of osmotic-related substances to reduce water loss. Therefore, the osmoregulation associated genes were analyzed. According to the GO enrichment analysis, asparagine biosynthesis occurred in KN27-3 at all five stages analyzed (Supplementary File S8), but it occurred in XLZ38 only at 30 hpi. Glutamine-dependent asparagine synthetase 1 (ASN1) is an important enzyme in the nitrogen metabolism cycle of plants and a key enzyme for asparagines synthesis. ASN1 (Gh_A09G2475) was much more expressed from 12 to 30 h in KN27-3 than in XLZ38 but showed a slight decrease afterward while maintaining a low level in XLZ38 (Figure 5a). ASP1 catalyzes transamination between aspartate and oxaloacetate, the expression of aspartate aminotransferase 1 (ASP1) (Gh_A09G0035) showed a notable increase in both KN27-3 and XLZ38 at 42 hpi and remained at similar levels at 54 h, but its fragments per kilobase of transcripts per million mapped fragments (FPKM) in KN27-3 was nearly twice that in XLZ38 at 42 and 54 h (Figure 5b).
Seed germination was accompanied by soluble sugar synthesis and catabolism. GDP-mannose biosynthetic process (GO: 0009298) and trehalose biosynthetic process (GO: 0005992) were enriched in KN27-3 at 12 hpi accompanied by the high up-regulation of the genes involved (Supplementary File S9). CYT1 encoded a GDP-mannose pyrophosphorylase/mannose-1-pyrophosphatase and TPS11 encoded an enzyme putatively involved in trehalose biosynthesis. CYT1 (Gh_A10G2006) and TPS11 (Gh_D08G0936) had higher expression levels at 12 hpi in KN27-3, and then lowered with the prolonged germination time of two varieties (Figure 5c,d). Raffinose catabolic process (GO: 0034484) was found to be active in KN27-3 at all analyzed stages except stage 4. The enriched genes were also up-regulated in XLZ38 (Supplementary File S9). SIP2 (Gh_D03G0964) encodes a raffinose-specific alpha-galactosidase that catalyzes the breakdown of raffinose into alpha-galactose and sucrose. Its expression level in XLZ38 was higher than in KN27-3 during the entire germination process (Figure 5e). The expression of both anabolic and catabolism-related genes showed higher expression levels from 12 to 30 hpi. Thus, the early up-regulation of the genes involved in the biosynthesis of mannose, trehalose, and raffinose in KN27-3 may provide useful protection against chilling stress.

DEGs Related to Osmoregulation
After the onset of chilling, plants often increase the accumulation of osmotic-related substances to reduce water loss. Therefore, the osmoregulation associated genes were analyzed. According to the GO enrichment analysis, asparagine biosynthesis occurred in KN27-3 at all five stages analyzed (Supplementary File S8), but it occurred in XLZ38 only at 30 hpi. Glutamine-dependent asparagine synthetase 1 (ASN1) is an important enzyme in the nitrogen metabolism cycle of plants and a key enzyme for asparagines synthesis. ASN1 (Gh_A09G2475) was much more expressed from 12 to 30 h in KN27-3 than in XLZ38 but showed a slight decrease afterward while maintaining a low level in XLZ38 (Figure 5a). ASP1 catalyzes transamination between aspartate and oxaloacetate, the expression of aspartate aminotransferase 1 (ASP1) (Gh_A09G0035) showed a notable increase in both KN27-3 and XLZ38 at 42 hpi and remained at similar levels at 54 h, but its fragments per kilobase of transcripts per million mapped fragments (FPKM) in KN27-3 was nearly twice that in XLZ38 at 42 and 54 h (Figure 5b). Values are the FPKM means of three independent replicates ± SE. The internodes in the same hour act as a comparison. Differences are indicated by "*", "NS" represents no difference.

DEGs Involved in Synthesis of Antioxidants and Antioxidative Enzymes
Antioxidants play an important role in plant resistance to low temperatures. Therefore, the expression of the genes involved in antioxidant biosynthesis was analyzed. The major genes are shown in Figure 6 and other related genes in Supplementary File S10. The spermine biosynthetic process (GO: 0006597) and polyamine biosynthetic process (GO: 0006595) were found to be significantly enriched in KN27-3 at 54 hpi. For example, S-adenosylmethionine decarboxylase (SAMDC) is a key enzyme in the polyamine metabolic pathway. Its expression was remarkably more up-regulated at 54 h in KN27-3 compared to XLZ38 (Figure 6a). In addition, Gh_A08G0862 encoding spermine synthase (SPMS) was expressed at a higher level in KN27-3 at 42 and 54 hpi Values are the FPKM means of three independent replicates ± SE. Bars marked with asterisks indicate differences change between the same hour (* p < 0.05, ** p < 0.01) according to Duncan's multiple range test using SPSS software. "NS" represents no difference.
Seed germination was accompanied by soluble sugar synthesis and catabolism. GDP-mannose biosynthetic process (GO: 0009298) and trehalose biosynthetic process (GO: 0005992) were enriched in KN27-3 at 12 hpi accompanied by the high up-regulation of the genes involved (Supplementary File S9). CYT1 encoded a GDP-mannose pyrophosphorylase/mannose-1-pyrophosphatase and TPS11 encoded an enzyme putatively involved in trehalose biosynthesis. CYT1 (Gh_A10G2006) and TPS11 (Gh_D08G0936) had higher expression levels at 12 hpi in KN27-3, and then lowered with the prolonged germination time of two varieties (Figure 5c,d). Raffinose catabolic process (GO: 0034484) was found to be active in KN27-3 at all analyzed stages except stage 4. The enriched genes were also up-regulated in XLZ38 (Supplementary File S9). SIP2 (Gh_D03G0964) encodes a raffinose-specific alpha-galactosidase that catalyzes the breakdown of raffinose into alpha-galactose and sucrose. Its expression level in XLZ38 was higher than in KN27-3 during the entire germination process (Figure 5e). The expression of both anabolic and catabolism-related genes showed higher expression levels from 12 to 30 hpi. Thus, the early up-regulation of the genes involved in the biosynthesis of mannose, trehalose, and raffinose in KN27-3 may provide useful protection against chilling stress.

DEGs Involved in Synthesis of Antioxidants and Antioxidative Enzymes
Antioxidants play an important role in plant resistance to low temperatures. Therefore, the expression of the genes involved in antioxidant biosynthesis was analyzed. The major genes are shown in Figure 6 and other related genes in Supplementary File S10. The spermine biosynthetic process (GO: 0006597) and polyamine biosynthetic process (GO: 0006595) were found to be significantly enriched in KN27-3 at 54 hpi. For example, S-adenosylmethionine decarboxylase (SAMDC) is a key enzyme in the polyamine metabolic pathway. Its expression was remarkably more up-regulated at 54 h in KN27-3 compared to XLZ38 (Figure 6a). In addition, Gh_A08G0862 encoding spermine synthase (SPMS) was expressed at a higher level in KN27-3 at 42 and 54 hpi compared to that in XLZ38, which remained unchanged (Figure 6b).

Endogenous Hormone Content and Related Gene Expression
Plant hormones are important for seed development and response to chilling stress; therefore, the contents of endogenous hormones at five germination stages were measured via high performance liquid chromatography-mass spectrometry (HPLC-MS). As shown in figure 7a, higher IAA content was detected in KN27-3 at all time points with a maximum level of 157.40 ng/g at 30 hpi, which was about four-fold higher than that of XLZ38. There was also more accumulation of cZR and tZR in KN27-3 compared to XLZ38 (Figure 7b,c). The GA3 content (Figure 7d) in the embryos of the two varieties showed a large difference at 42 hpi, which was about three times higher in KN27-3 than that of XLZ38. While ABA content in XLZ38 was higher than in KN27-3 at all five germination stages with a significant difference at 12, 42, and 54 hpi (Figure 7e).
DEGs involved in CKT, GA, and ABA metabolism were identified, and their changes during seed germination were consistent with the accumulation of hormones. TP/ADP isopentenyl transferase 2 (IPT2) catalyzes the rate-limiting step of CTK biosynthesis. Its encoding gene The Y axis shows the FPKM means of three independent replicates ± SE. Bars marked with asterisks indicate differences change between the same hour (* p < 0.05, ** p < 0.01) according to Duncan's multiple range test using SPSS software. "NS" represents no difference.

Endogenous Hormone Content and Related Gene Expression
Plant hormones are important for seed development and response to chilling stress; therefore, the contents of endogenous hormones at five germination stages were measured via high performance liquid chromatography-mass spectrometry (HPLC-MS). As shown in Figure 7a, higher IAA content was detected in KN27-3 at all time points with a maximum level of 157.40 ng/g at 30 hpi, which was about four-fold higher than that of XLZ38. There was also more accumulation of cZR and tZR in KN27-3 compared to XLZ38 (Figure 7b,c). The GA3 content (Figure 7d) in the embryos of the two varieties showed a large difference at 42 hpi, which was about three times higher in KN27-3 than that of XLZ38. While ABA content in XLZ38 was higher than in KN27-3 at all five germination stages with a significant difference at 12, 42, and 54 hpi (Figure 7e).

Discussion
Plants growing in the natural environment are constantly subjected to a variety of abiotic stresses, such as drought, salinity, heat, cold, and heavy metals. Plant stress response and adaptation are extremely complex, involving molecular, cellular and physiological changes [29]. Based on the present study, seed germination rates of varieties, KN27-3 and XLZ38 showed a significant difference at low temperatures, indicating different responses to chilling (Supplementary File S1). Our results demonstrated that early and rapid response to low temperature in KN27-3 played an important role in promoting its hypocotyl elongation. KN27-3 ensured the stability of the intracellular environment to reduce chilling damages. It is known that enhancing the biosynthesis of different types of compatible organic solutes is one of the most common responses of plants against stress [30]. Asparagines as the main form of transport of amino acids in cotton was not only an important carrier for transporting new synthetic nitrogen in plants but also a potential osmotic adjustment [31,32]. Stewart found that asparagine might have a function similar to proline-regulating cell penetration [31,33]. ASP1 and ASN1 had high expression levels during the germination of KN27-3. ASP1 and ASN1 catalyze the production of aspartate and asparagine, DEGs involved in CKT, GA, and ABA metabolism were identified, and their changes during seed germination were consistent with the accumulation of hormones. TP/ADP isopentenyl transferase 2 (IPT2) catalyzes the rate-limiting step of CTK biosynthesis. Its encoding gene Gh_D09G0911 had a higher FPKM value in KN27-3 than in XLZ38. Its expression reached its peak at 30 hpi and declined later (Figure 7f). Moreover, Gh_D05G1813, encoding oxidase/dehydrogenase, which is involved in CTK degradation, was up-regulated in XLZ38 and was significantly higher than in KN27-3 at 12 and 42 hpi (Figure 7g). Gh_D09G0042 encoding GA20 oxidase homolog GA20ox1B was expressed during germination. It is involved in GA biosynthesis. GA20ox1B had significantly higher expression in KN27-3 than that in XLZ38 at 12 and 18 hpi. In KN27-3, its expression was stable but significantly down-regulated in XLZ38 at 54 h compared to 42 h (Figure 7h). Abscisic acid 8'-hydroxylase 1 (CYP707A1) is an important enzyme for ABA inactivation. Gh_A08G1344 was down-regulated in XLZ38 from 30 h to 42 h, while it was not changed in KN27-3. Its expression first increased and then decreased by 18 h in both varieties, although the decrease was faster in XLZ38 than in KN27-3 ( Figure 7i). In summary, the expression dynamics for these genes was consistent with the hormone accumulative dynamics in the two varieties. IAA, GA3, and CTK played a positive regulatory role, and ABA played a negative regulatory role in KN27-3 during the low temperature germination.

Discussion
Plants growing in the natural environment are constantly subjected to a variety of abiotic stresses, such as drought, salinity, heat, cold, and heavy metals. Plant stress response and adaptation are extremely complex, involving molecular, cellular and physiological changes [29]. Based on the present study, seed germination rates of varieties, KN27-3 and XLZ38 showed a significant difference at low temperatures, indicating different responses to chilling (Supplementary File S1). Our results demonstrated that early and rapid response to low temperature in KN27-3 played an important role in promoting its hypocotyl elongation. KN27-3 ensured the stability of the intracellular environment to reduce chilling damages. It is known that enhancing the biosynthesis of different types of compatible organic solutes is one of the most common responses of plants against stress [30]. Asparagines as the main form of transport of amino acids in cotton was not only an important carrier for transporting new synthetic nitrogen in plants but also a potential osmotic adjustment [31,32]. Stewart found that asparagine might have a function similar to proline-regulating cell penetration [31,33]. ASP1 and ASN1 had high expression levels during the germination of KN27-3. ASP1 and ASN1 catalyze the production of aspartate and asparagine, respectively. In our results, more aspartic acid and asparagine were produced in the seed germination of KN27-3 to regulate the osmotic pressure of cells, which may protect the germination process from the chilling conditions in addition to producing more ammonium.
The biosynthesis process of trehalose and GDP-mannose were enriched during the germination of the two varieties. CYT1 encoding GDP-mannose pyrophosphorylase/mannose-1-pyrophosphatase had a higher expression level in KN27-3 at 18 h. Over-expression of tomato GMPase in tobacco reduced H 2 O 2 and O 2− content in plants under high temperature or low temperature stress [34]. GMPase was significantly up-regulated in KN27-3, which may reduce the stress-induced peroxidative damage and produce the structural component, GDP-mannose, for the cell wall. The synthesis of carbohydrates was also the substrate for the synthesis of trehalose [35]. Chickpea with accumulative trehalose performed better under both optimal temperature conditions and chilling stress because trehalose protects the plants from oxidative damage and helps to maintain carbon assimilation and seedling growth [36]. In the present study, four TPS11/9 homologous genes had higher expression levels from 12 to 30 h in KN27-3 similar to ASN1 and CYT1. The TaTPS11 transcription level was positively correlated with wheat cold tolerance [37]. The expression of AtTPS6, 7, 9, and 10 persists towards the root base and stops at the root/hypocotyl junction, promoting hypocotyl elongation in Arabidopsis [38]. Cold stress-induced raffinose synthesis can regulate plant permeability for cold resistance [39]. Raffinose family oligosaccharides (RFO) facilitate chickpea seed germination [40]. The catabolism of raffinose was enriched to almost the entire imbibition process. The degree of raffinose decomposition in XLZ38 was higher than that of KN27-3, which may lead to more raffinose in KN27-3 and improve stress tolerance (Figure 5e). Therefore, during the imbibition stage, more GDP-mannose and trehalose may be synthesized, while more raffinose may be broken down in the early germination stage (from 12 to 30 hpi) of KN27-3, which helps with osmotic adaptation under low temperatures during seed germination.
The embryo development is often associated with the accumulation of reactive oxygen species (ROS) [41]. ROS is also generated in tissues under chilling stress. And more ROS accumulation can cause different degrees of damage to proteins, lipids, nucleic acids, and other biological macromolecules, which may lead to the change of biofilm fluidity and inactivation of biological enzymes, hinder protein synthesis and degradation, and cause DNA damage, etc. [42]. Therefore, an effective antioxidant mechanism is regarded as a prerequisite for obtaining vigorous seed. The damages caused by chilling stress are correlated with unbalanced intracellular oxidative systems in plants [43]. All polyamines have the ability to bind DNA for protection. They also contribute to plant stress tolerance by minimizing oxidative damage and maintaining membrane structure [44]. Spermidine (Spd) was shown to enhance clover seed germination under low temperature [45]. In our study, DEGs involved in SAMDC and SPMS were expressed in both cultivars, but higher expression levels were observed in KN27-3 than in XLZ38 during seed germination. Multiple studies have demonstrated that Spd is generally correlated with the hypocotyl elongation rate, such as in chickpea [46], corn [47], and in mung bean [48]. Thus, higher expression of SAMDC and SPMS may increase the accumulation of Spd in germinating KN27-3 seeds under chilling conditions. Stress, for example, water shortage [49] and K deficiency [50], can increase the activities of antioxidant enzymes for scavenging ROS, while low temperatures may inhibit enzyme activity [51]. Thus, it is likely that the fine-tuned regulation of antioxidant enzyme activities was involved in promoting hypocotyl elongation under chilling conditions in KN27-3. The expression of CAT1, GPX4, DHAR2, and APX1 in KN27-3 was higher than that in XLZ38 ( Figure 6). Rapid and efficient expression of CAT1 and GPX4 was more pronounced at 12 hpi, while APX1 was more strongly expressed at later time points (42 and 54 h), and DHAR2 was expressed throughout the germination stage (12, 42, and 54 h). GPX4 is present in various cell compartments and it closely cooperates with enzymes, such as SOD and CAT1, to coordinate and eliminate excessive free radicals in plants and improve plant resistance [52]. Under drought and cold stress conditions, glutathione (GSH) and APX are responsible for clearing drought and cold-induced ROS detoxification, and their increased expression promotes higher germination rates, increased root length and higher accumulation of fresh weight [53]. Studies have shown that higher CAT activity can reduce H 2 O 2 content in plants [54] and that elevated CAT activity can improve rice tolerance and seed germination rate at low temperatures [55]. Our data suggest that KN27-3 has a higher ability to alleviate the harmful effects of ROS under cold stress than XLZ38, confirming that rapid accumulation of antioxidants is an important factor in determining cottonseed cold stress tolerance.
Various studies have shown that plant hormones regulate seed germination. ABA and GAs play the most important roles, while IAA is present in the seedling radicle tip during and after germination, and CTK is activated during germination [56,57]. Our data suggest that IAA and CTK accumulated during germination of KN27-3 under chilling stress, which serves to ensure the growth of embryos (Figure 7a-c). When plants are stressed, they maintain growth by balancing both IAA and CTK hormonal pathways along with ROS signals [58]. IPT2 catalyzes the decomposition of prenylated tRNA, which is a key step in cZT synthesis. The occurrence of zeatin cis-trans isomerase activity suggested that the tRNA-mediated pathway might also contribute to the synthesis of tZ-type CKs via cZ-type CKs [59]. CKX enzymes catalyze their reversible degradation of CTK by oxidative side chain cleavage [58]. Our results showed that the transcripts Gh_D09G0911, encoding IPT2, and Gh_D05G1813, encoding CKX, were up-regulated and down-regulated, respectively, to produce more CTK in KN27-3. During the germination of dicotyledons, GA mainly stimulates germination by promoting radicle elongation and penetration of the seed coat [60]. Our data suggest that the accumulation of GA3 in KN27-3 began to increase at 30 h, and was significantly higher than that of XLZ38 at 42 h. Expression of the gene encoding GA20ox1B also significantly increased in KN27-3 at 42 h, and there was a significant difference at 54 hpi (Figure 7d,h). This suggests that the synthesis of GA3 positively regulated the germination process of KN27-3 under chilling stress. ABA is a key phytohormone that modulates plant growth and development as well as abiotic and biotic stress responses. ABA accumulates in the developing embryo and regulates seed development, seed maturation, and seed dormancy [61]. ABA endogenous content was higher in XLZ38 (Figure 7e). CYP707A1 is one of the two important ABA catabolic genes [62]. Gh_A08G1344 transcription levels were down-regulated during germination in XLZ38 (Figure 7i). We speculated that ABA inhibits seed germination while IAA, CTK, and GA3 positively regulates seed germination, suggesting that hormones and transcription of metabolic genes together regulate the germination process of seeds.
Antioxidase and catalytic enzymes involved in the metabolism of the two varieties have different expression levels of transcription under the same temperature. In the next research, we will continue to study the reasons for the different antioxidant enzyme activity of different varieties at the same temperature, and explore the relationship between the level of antioxidant enzyme activity and the content of H 2 O 2 and O 2− content in cottonseed under chilling conditions.

Plant Material and Seed Germination Conditions
Two varieties were used in the experiments: a cold-susceptible variety XinLuZao38 (XLZ38) and a cold-tolerant material KenN27-3 (KN27-3). Both varieties are suitable for planting in the northwest inland cotton region (Xinjiang, China). Their seeds were harvested in Xinjiang in 2017.
The experiments were performed in a germination box (185 × 145 × 62 mm) with 12 holes. Plants were cultured in sand with a water content of 13%. Five seeds were placed in each hole, a total of 20 seeds per variety were planted. Three replicates were performed.
Cotton seeds were grown (28 • C/25 • C, 12 h/12 h, day/night) for 4 days to analyze the seed germination phenotype. The hypocotyls lengths were measured at 2, 3, and 4 days after sowing. Germination tests under chilling stress were performed in a growth chamber with~150 µmol m −2 s −1 of fluorescent light, a 12 h day/12 h night cycle, and 16 • C/4 • C, and hypocotyl length was determined 6, 8, 10, and 12 days after the initiation. Fresh seed weight was measured the weight every 6 h, and water absorption per unit mass of cotton seed (∆M) was calculated as ∆M = (M i − M 0 )/M 0 , in which M 0 is the dry weight of 20 seeds, M i is the weight of 20 seeds at different germination times. Seed embryos of KN27-3 and XLZ38 (named K1-5, X1-5) were collected at five germination stages: 12, 18, 30, 42, and 54 h (Supplemental Figure S2), frozen in liquid nitrogen and stored at −80 • C for RNA-Seq and hormone determination.

RNA Isolation, Library Construction, and RNA-Seq
After total RNA was extracted, eukaryotic mRNA was enriched by Oligo (dT) beads, while prokaryotic mRNA was enriched by removing rRNA by Ribo-Zero TM Magnetic Kit (Epicentre). Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP, and buffer. Then, the cDNA fragments were purified with a QIAquick PCR extraction kit, end repaired, poly(A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq TM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Sequencing Analysis and Differential Expression Analysis
Clean reads were obtained by removing adaptor sequences, more than 10% N bases and low-quality (Q ≤ 20) reads with more than 50% bases from each data set to gain more reliable results. The reads were mapped onto the transcriptome assembly by TopHat 2 (V2.1.1) to map the reads to the Gossypium hirsutum L. genome [63]. Read counts per gene were expressed as the expected number of fragments per kilobase of transcripts per million mapped fragments (FPKM), and unigene abundance differences between the samples were calculated based on the ratio of the FPKM values and false discovery rate (FDR). Genes with FDR ≤ 0.05 and FPKM ≥ 10 were considered DEGs.
GO classification was performed via WEGO, (http://wego.genomics.org.cn/cgi-bin/wego/index.pl), and the GO distributions of DEGs were then obtained from three levels: biological process and cellular component, molecular function. For each Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, the numbers of DEGs were compared to the entire reference gene set by hypergeometric tests to find out the pathways enriched with regulated genes. KEGG enrichment analyses were adjusted using the Bonferroni correction, and a corrected p. adjust ≤ 0.05 was chosen as the threshold value for determining significantly enriched GO terms and KEGG enrichment pathways.

Statistical Analysis
The statistical analysis was conducted with SPSS software 22 (SPSS. Inc) for Windows. One-Way ANOVA analysis and Tukey's honestly significant difference test were performed to test for differences between the means of the embryo data for KN27-3 and XLZ38 at different stages using a significance level of p < 0.05.

Conclusions
In total, 7535 DEGs were identified and analyzed for their potential role between KN27-3 and XLZ38 at five imbibition stages using GO enrichment and KEGG pathway analysis. Elevated levels of IAA, CTK, and GA and reduced ABA in KN27-3 (tolerant genotype) might interact with one another to alter energy metabolism to maintain nitrogen and carbon balance, which might contribute to seed germination. Ultimately, asparagine, GDP-mannose, trehalose, and raffinose metabolic process, the rapid accumulation of spermine, and higher antioxidase activity of genotype KN27-3, may explain its improved performance under chilling stress.

Conflicts of Interest:
The authors declare no conflicts of interest.