ZM36 is provided by the Cotton Molecular Breeding Laboratory of Shihezi University. Cotton seedlings are cultured in an artificial climate box. These are cotton seeds that have undergone germination in a nutrient bowl. The substrate is peat and vermiculite, with a ratio of 3:1. The plants were cultivated under 24/22°C (day/night) and a photoperiod scheme of 16/8 h of light/darkness. Seedlings at the two-leaf stage were processed at low temperature in an incubator with adjustable temperature settings.

The processing method is as follows: cotton seedlings with consistent growth are grown for 24 h at 25°C (T25, control), 15°C (T15), 10°C (T10), and 4°C (T4), then sampled, frozen, and stored at -80°C. The experiment setup had three replicates, with five identical cotton seedlings as one replicate, and each treatment setup had three biological replicates.

Portable photosynthetic instrument LI-6400XT (LI-COR, USA) is used to measure the photosynthetic performance parameters, with the instrument’s built-in red and blue light source selected and the light intensity set to 1,000 μ Mol m^-2 s^-1. HandyPEA-100 (UK) was used to measure the fluorescence parameters, and these were measured at the same leaf position of seedlings with similar growth. The conductivity meter method was used to measure the relative conductivity, and the acidic ninhydrin colorimetric method was employed to measure the proline content.

RNA Purification Kit (Tiangen, Beijing) was used to isolate RNA from the leaves of 12 samples (ZM36) according to the manufacturer’s instructions. The RNA isolated from each sample was then used to construct RNA-seq libraries using NEBNext Ultra RNA Library Prep Kit. RNA-seq was conducted on an Illumina Hiseq 4000 platform with 150-bp paired-end reads (Novogene, Tianjin, China).

qRT-PCR was carried out by using SYBR Green (Roche, Rotkreuz, Switzerland) on Light Cycler 480II (Roche) with default parameters. All primers used for the validation experiments were designed with Primer5 software and are shown in Supplementary Table S1 . The GhUBQ7 (DQ116441.1) gene served as an internal control to normalize differences between samples. The qRT-PCR conditions were as follows: initial denaturation at 96°C for 5 min, denaturation at 96°C for 15 s for a total of 41 cycles, annealing at 62°C for 16 s, and extension at 70°C for 18 s. The relative expression levels of genes from three biologically independent experiments were calculated using the 2^-ΔΔCTmethod (Livak and Schmittgen, 2001).

Clean reads were mapped to the reference genome of G. hirsutum (Hu et al., 2019) using TopHat (v2.0.12). The level of gene expression was measured by fragments per kilobase of exon model per million mapped fragments (FPKM) (Trapnell et al., 2010). DESeq2 (v1.18.0) (Wang et al., 2010) was used to identify differentially expressed genes (DEGs) with the criteria of an adjusted log₂ (fold change) ≥3. The STEM software was used to classify the gene expression patterns at various points into differential clusters. Using the cluster Profiler package (version 3.18.1), DEG enrichment analysis was conducted using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG).

We used Blast (blast x) for the sequences of DEGs with the genomes of related species to obtain the predicted PPI for these DEGs (the protein interactions exist in the STRING database: http://string-db.org/. Then, the PPI of these DEGs was visualized using Cytoscape (Shannon et al., 2003).

Tobacco rattle virus (TRV) vectors, pTRV1 and pTRV2, were used in the VIGS experiments, and TRV:: GhCHLI was used as a positive control as previously reported (Li et al., 2023b). A 300-bp fragment specific to GhPLATZ (GH_A09G2161) was amplified by PCR from low-resistance ZM36 with gene-specific primers ( Supplementary Table S2 ). VIGS was performed with the same procedures as previously described (Gao et al., 2013). Approximately 12 days after infiltration, the leaves of five TRV::00 and TRV:: GhPLATZ plants were collected to analyze the expression level of GhPLATZ by qRT-PCR. A total of 60 TRV::00 and TRV:: GhPLATZ plants at the two-leaf-stage were subjected to low temperature treatment to compare their phenotypic response.

ZM36 was tested for Pn, Gs, Ci, Tr, Fv/Fm, Piabs, Fv/Fo, and Fo/Fm in four different temperature points (T25, T15, T10, and T4) ( Figure 1A ). Pn, Gs, Tr, Fv/Fm, Fv/Fo, and Piabs significantly decreased, but Ci and Fo/Fm significantly increased under cold stress. The highest value of Pn is 12.75 at T25, the minimum is 1.18 at T4, reduced by 10.8 times; Gs decreased by 26.28 times from T25 (0.3278) to T4 (0.0124); Tr decreased by 20.75 times from T25 (4.9022) to T4 (0.2362); Fv/Fm was 0.7951 and 0.311 at T25 and at T4, respectively, reduced by 2.56 times; and Piabs decreased by 37.81 times from T25 (4.2996) to T4 (0.1137). However, Ci increased by 1.64 times from T25 to T4, and Fo/Fm also increased by 3.25 times from T25 (0.2049) to T4 (0.6659) ( Figure 1B ; Supplementary Table S2 ). This fully demonstrates that low temperature can have a significant impact on cotton growth, photosynthesis, and yield of cotton, thereby affecting its yield and quality.

The number of raw reads of the samples was between 40,476,800 and 49,397,772, and the number of clean reads after filtering was 20,238,400–24,698,886 ( Supplementary Table S3 ). The GC percentage of each sample was between 43.46% and 44.34%, while Q30 was between 90.16% and 95.04% ( Supplementary Table S3 ). Pearson correlation coefficient was used to perform a correlation test on the samples, and the correlation between the three replicates was greater than 0.9, indicating the reliability of the data ( Supplementary Figure S1A ). Principal component analysis (PCA) was used to test the similarity between samples. The first principal component (PC1) accounted for 43%, while the second principal component (PC2) accounted for 26.3% ( Supplementary Figure S1B ). The volatility of log10^FPKM is similar across all samples ( Supplementary Figure S1C ), indicating that the sequencing quality was high and suitable for subsequent analysis. Using fold change ≥3 and FDR <0.05 as the threshold, 6,096 DEGs (2,430 upregulated and 3,666 downregulated), 8,545 DEGs (4,473 upregulated and 4,072 downregulated), and 3,322 DEGs (1,659 upregulated and 1,663 downregulated) were screened in T15 vs. T25, T10 vs. T25, and T4 vs. T25, respectively, indicating that gene expression levels in cotton varied after different low temperature treatments ( Figures 2A, C ). Using Wayne analysis, a total of 13,183 DEGs were found, of which 2,365 DEGs were specifically responsive to 15°C, 4,600 DEGs were specifically responsive to 10°C, and 1,972 DEGs were specifically responsive to 4°C ( Figure 2B ). The DEGs were divided into eight significant clusters by K-means clustering analysis ( Figure 2D ). The genes of clusters 1, 2, and 3 showed a trend of rising first and then remaining unchanged or decreasing slightly in response to low temperature stress. However, the genes of clusters 5, 6 and 7 had an expression trend opposite to that of clusters 1, 2, and 3. The genes of cluster 4 showed a decreasing trend, and cluster 8 had an expression trend opposite to that of cluster 4 ( Figure 2D ). The DEGs of these clusters would be expected to be linked to the low temperature resistance of ZM36.

The expression patterns of 12 genes were validated using qRT-PCR, including six significant downregulation expressed genes (Gh_D09G2404, Gh_A02G0898, Gh_A06G0948, Gh_D10G2061, Gh_A07G1351, and Gh_D10G1486), six significant upregulation expressed genes (Gh_D05G2845, Gh_D06G1877, Gh_A12G1885, Gh_A13G1204, Gh_A01G1052, and Gh_A05G1366). All these genes showed the same expression trend between qRT-PCR ( Figure 3 ) and RNA-Seq ( Figure 3 ) at T25, T15, T10, and T4, confirming the reliability of the RNA-seq result.

GO enrichment analysis of DEGs was performed to determine the functions of the distinct transcripts differentially expressed in ZM36 after low temperature stress. In T15 vs. T25, the GO terms such as response to ATP binding (GO:0005524), response to oxygen-containing compound (GO:1901700), cell wall organization (GO:0071555), microtubule-based process (GO:0007017) were commonly enriched ( Figure 4A ). The GO terms such as response to ubiquitin–protein transferase activity (GO:0004842), protein ubiquitination (GO:0016567), response to external stimulus (GO:0009605), ubiquitin protein ligase activity (GO:0061630), and response to water deprivation (GO:0009414) were commonly enriched in T10 vs. T25 ( Figure 4B ). In T4 vs. T25, the DEGs are mainly enriched in ATP binding (GO:0005524), plasma membrane (GO:0005886), transcription factor activity, sequence-specific DNA binding (GO:0003700), protein kinase activity (GO:0004672), and regulation of transcription, DNA-templated (GO:0006355) ( Figure 4C ).

The KEGG enrichment analysis of all the DEGs resulted in 2,454 pathways. There were many significant changes in pathways related to circadian rhythm—plants, beta-alanine metabolism, zeatin biosynthesis, fatty acid metabolism, biosynthesis of amino acids, MAPK signaling pathway–plant, glutathione metabolism, and plant–pathogen interaction ( Figures 5A–C ), which are mainly associated with plant growth and development and response to stresses. The KEGG pathways enriched in T4 vs. T25 included MAPK signaling pathway–plant, plant–pathogen interaction, photosynthesis–antenna proteins, photosynthesis, flavonoid biosynthesis, and plant hormone signal transduction ( Figure 5C ). The enriched pathways also suggested that genes related to MAPK cascade signaling might be related to the cotton response to low-temperature stress. The pathways related to MAPK signaling pathway–plant were enriched in two temperature points (T10 and T4) ( Figures 5B, C ). The genes involved included those encoding MAPKKK18 (GH_A03G0386), serine/threonine protein kinase OXI1 (GH_A07G2243), abscisic acid receptor PYR1 (GH_A12G2288), ultraviolet-B receptor UVR8 (GH_A05G3668), and PLATZ transcription factor family protein (GH_A09G2161). The above-mentioned results indicated that the response of cotton to low temperature stress is governed by a complex gene network that regulates multiple metabolic pathways.

In T15 vs. T25, T10 vs. T25, and T4 vs. T25 groups, 471 transcription factors (49 transcription factor families), 767 transcription factors (59 transcription factor families), and 341 transcription factors (41 transcription factor families) were found to be involved in the response to different degrees of cold stress, respectively. To predict the transcription factors of all DEGs, 1,119 transcription factors were found in 59 transcription factor families ( Supplementary Figure S2A ). Among them, different members of major gene families such as AP2/ERF-ERF (133 genes), MYB (108 genes), bHLH (78 genes), and NAC (73 genes) were upregulated or downregulated under different cold stresses ( Supplementary Figures S2B–E ).

The phosphatidylinositol signaling system regulates many physiological processes, such as growth, cytoskeleton rearrangement, and membrane transport. The DEGs of the phosphatidylinositol signaling system (ko04070) were identified to be 32 genes (21 upregulated, 11 downregulated), 40 genes (29 upregulated, 11 downregulated), and 22 genes (18 upregulated, four downregulated) in T15, T10, and T4, respectively ( Supplementary Figure S3 ). Among them, all others were upregulated except for GH_D10G0541 in two 1-phosphatidylinositol-3-phosphate 5-kinase (FAB1) after cold stress ( Supplementary Figure S3A ). GH_A10G0515 was upregulated 4.60, 3.88, and 5.66 times after T15, T10, and T4, respectively ( Supplementary Figure S3A ). Six phosphatidylinositol-specific phospholipases (PLCD) were downregulated except for GH_A06G1884 ( Supplementary Figure S3A ). Three diacylglycerol kinases (DGK) were upregulated under different cold stresses, except for GH_A12G2233, and GH_D12G2002 was upregulated 4.40, 15.03, and 2.40 times after T15, T10, and T4, respectively ( Supplementary Figure S3B ). Most of the 25 calmodulin (CALM) genes were upregulated, especially GH_A04G1762 and GH_D12G1967. They were upregulated by 104.81, 17.55, 4.12 times and 74.50, 76.91, and 2.80 times after T15, T10, and T4, respectively ( Supplementary Figure S3B ).

Photosynthesis is sensitive to cold stress. KEGG indicate that many DEGs are involved in photosynthesis metabolic pathways. In this study, 34 differentially expressed genes were enriched under different cold stresses in photosynthesis (ko00195) ( Figure 6A ). One psaA, one psaL, and one psaO gene were found in the photosystem I complex. Among them, psaA (GH_D03G0877) was upregulated 1.93, 1.20, and 2.99 times after T15, T10, and T4, respectively. Two psbA, one psbD, three psbB, and one psbH gene in the photosystem II complex were upregulated after T15, T10, and T4, while four psb27 genes were downregulated. One petB gene and one petA gene were found in the cytochrome b6-f complex, which were upregulated under T4 for 3.13, 3.52, and 4.61 times, respectively. Seven DEGs related to photosynthetic electron transfer were downregulated after T15, T10, and T4, indicating that cold stress has a significant inhibitory effect on electron transfer in photosynthesis. The 27 DEGs (five LHC I and 22 LHC II) were found in the photosynthesis antenna proteins (ko00196) ( Figure 6B ). Except for three Lhcb4 genes that were downregulated after T15 and T10 and upregulated after T4, most genes were upregulated after T15 and T10 and downregulated after T4. The 24 genes were upregulated after T15 and T10 and downregulated after T4, and three Lhcb4 genes were exactly the opposite. A total of 34 DEGs were found in carbon fixation in photosynthetic organisms (ko00710) pathway ( Figure 6C ). The C4 pathway is enriched to 15, and the Calvin cycle is enriched to 19. Among them, one ribose 5-phosphate isomerase A (GH_A01G1655) in the Calvin cycle was upregulated by 2.00, 2.75, and 6.72 times after T15, T10, and T4 stress, respectively, and one phosphoglycerate kinase (GH_A03G0406) was upregulated 2.03, 20.07, and 4.49 times after T15, T10, and T4 stress, respectively. At the same time, we also discovered 35 DEGs ( Figure 6D ). Six genes were in carotenoid biosynthesis (ko00906) β. Most of the beta carotene isomerase (DWARF27) was upregulated, especially GH_A07G2351 and GH_D07G2294, which were upregulated 4.99, 8.27, and 4.70 times and 4.06, 6.08, and 7.14 times after T15, T10, and T4 stress, respectively. Beta carotene 3-hydroxylase (crtZ) was upregulated in expression. Three zeaxanthin epoxidase (ZEP) genes were upregulated, with GH_D01G1934 upregulated 7.35, 40.20, and 1.51 times after T15, T10, and T4 stress, respectively. Four 9-cis-epoxycarotene dioxygenase (NCED) genes were upregulated, with GH_D13G1744 upregulated 12.84, 58.35, and 2.99 times after T15, T10, and T4 stress, respectively.

The 526 DEGs were detected under three different cold stress conditions, indicating that the DEGs may have been involved in adapting to stress at different low temperature points, and participating in the same pathway. Therefore, we defined 526 co-expressed DEGs as key cold-resistant genes. Further utilizing the STRING database to predict the interrelationships between 526 proteins, Cytoscape software was used for visualization processing, and five hub genes were identified based on the criteria of degree value ≥16 ( Figure 7 ). Four genes were significantly upregulated under cold stress (GH_D01G0766, GH_A05G3668, GH_A01G0766, and GH_A09G2161), while one gene was significantly downregulated under cold stress (GH_A02G1336).

To get insight on the potential function of GhPLATZ, we knocked down the expression level of GhPLATZ in ZM36 using VIGS with cotton seedlings treated by TRV:: GhCHLI as positive control of the VIGS experiment ( Figure 8A ). Compared to TRV2::00 plants, TRV:: GhPLATZ plants had a significantly low expression level of GhPLATZ ( Figure 8B ; Supplementary Table S4 ), suggesting the successful inhibition of GhPLATZ by VIGS. After about 10 days of cold stress, the TRV:: GhCHLI plants showed a yellowing phenotype ( Figure 8A ), indicating that the VIGS system was functioning properly. To confirm the effectiveness of the silencing system, cotton seedlings were subjected to 15°C (low temperature treatment) to verify the resistance of cotton to low temperature stress after inhibiting GhPLATZ expression. After 48 h of treatment at 15°C, it was found that TRV:: GhPLATZ plants were more sensitive to low temperature compared to TRV::00 plants, with more significant leaf wilting and dehydration ( Figure 8C ; Supplementary Table S4 ), indicating that inhibiting the expression of GhPLATZ would reduce cotton’s resistance to low temperature. The electrical conductivity of TRV : GhPLATZ plants was 49.82% ± 0.03, significantly higher than that of TRV:00 plants ( Figure 8D ; Supplementary Table S4 ). Consequently, the proline of TRV : GhPLATZ plants was higher than that of TRV:00 plants at 87 and 56 ug^-1 g, respectively ( Figure 8E ; Supplementary Table S4 ).

Cold stress can induce photoinhibition, reducing the absorption and capture capacity of the light system, electron transfer efficiency, and fixed quantum efficiency of carbon dioxide. In this study, the photosynthetic system was significantly damaged, inhibiting the normal progress of photosynthesis and leading to the generation of photoinhibition effects. Especially under 4°C, the photoinhibition effect was more severe. The downregulation of LHCA and LHCB genes under cold stress reduces the ability to absorb and capture light energy (Wilson et al., 2006). The downregulation of key enzyme activity and gene transcriptional expression levels in the Calvin cycle may be the reason for the decrease in plant carbon assimilation efficiency and RuBP regeneration rate caused by low temperature (Hussain et al., 2021). In this study, 27 DEGs were enriched in the photosynthesis antenna protein pathway under cold stress, while LHCA and LHCB genes were upregulated, which may be an adaptive protective response. The 34 DEGs are involved in the carbon fixation pathway in photosynthetic organisms, of which 19 are related to the Calvin cycle. The DEGs related to photosynthetic electron transfer are downregulated, indicating severe damage to the photosynthetic system, and reduce the efficiency of carbon assimilation under cold stresses. However, some DEGs were upregulated, although the Pn and Fv/Fm values decreased in T4. It is speculated that ZM36 can resist low temperature stress and alleviate oxidative damage by regulating the expression of photosynthetic-related genes. In addition, proline is a common osmoregulatory substance in plants, which can maintain the osmotic balance between the protoplast and the environment, thereby alleviating the damage caused by low temperature. Plants accumulate higher levels of proline when subjected to stress (Charest and Ton Phan, 1990; Shen et al., 2020). This study found that the proline content in silenced plants decreased by 30.1% compared to the control plants after low temperature stress, indicating a significant decrease in the cold tolerance of silenced plants. This suggests that enhancing the expression of photosynthesis-related genes is not enough to maintain a certain level of photosynthetic ability under cold stress. The specific mechanism values need further in-depth research.

When plants are subjected to cold stress, transcription factors (TFs) bind to specific cis regulatory elements in the promoter to regulate target genes related to cold resistance, thereby enhancing cold tolerance (Knight and Knight, 2012; Zhao et al., 2015). In recent years, TF families such as AP2/ERF, MYB, bHLH, and NAC have received widespread attention as key regulatory factors in plant stress response (Mehrotra et al., 2020; He et al., 2023). AP2/ERF plays a crucial regulatory role in response to low temperature stress. In addition, AtMYB15 plays a negative regulatory role in regulating cold resistance in Arabidopsis (Wang et al., 2019). The NAC family is the largest specific transcription factor family in plants, playing a crucial regulatory role in plant growth, development, and response to abiotic stress (Tran et al., 2010). In Arabidopsis, overexpression of HuNAC20 and HuNAC25 enhances tolerance to cold stress by altering the expression of cold response genes in transgenic plants (Hu et al., 2022). In this study, the number of members belonging to the AP2/ERF transcription factor family was the highest. At the same time, MYB transcription factors also play a crucial role in low temperature response. Most members of the MYB family show a downward trend under low temperature stress, and most genes in the NAC family show varying degrees of upregulation under different low temperature stress conditions, such as GH_D02G1383 and GH_A03G1198, which were upregulated hundreds of times under low temperature stress, highlighting the important role of the NAC family in cotton’s resistance to low temperature stress.

Cold resistance is a complex physiological and biochemical process involving gene regulation. Plants improve cold resistance by coordinating the expression of multiple genes (Umer et al., 2020) UVR8 is a photoreceptor that specifically absorbs UV-B light (Kliebenstein et al., 2002). The bZIP transcription factor HY5 is activated by UVR8 binding to COP1, thereby inducing the expression of multiple metabolic pathway genes (Favory et al., 2009). The tomato UVR8 gene participates in UV-B-induced cold tolerance by upregulating CuZnSOD, FeSOD, and CAT1 genes (Jiang et al., 2022). UVR8 regulates plant response to UV-B light by interacting with various transcription factors such as COP1 and HY5, including light morphogenesis, secondary metabolism, and adaptability to environmental stress. These transcription factors may affect the expression of genes related to ascorbic acid synthesis and metabolism, thereby affecting the level of ascorbic acid in plants (Lin et al., 2020; Podolec et al., 2021). The zinc finger transcription factor PLATZ is widely present in plants and plays an important role in regulating plant growth and development and responding to abiotic stress. The transgenic Arabidopsis thaliana with GhPLATZ1 may promote seed germination and seedling formation under salt stress by increasing the GA and ethylene content and reducing the ABA content (Han et al., 2022). It was found that PhePLATZ23 and PLATZ27 are highly responsive to cold stress and play an extremely important role in regulating bamboo’s response to external environmental stimuli (Zhang et al., 2022). UVR8 mainly involves the perception and transduction of light signals, while the PLATZ transcription factor is more related to plant transcriptional regulation and stress response. Although both play important roles in plant life activities, their mechanisms of action and biological functions are different, and there is currently no clear evidence to suggest a direct interaction or functional connection between them. Future research may reveal whether there are some unknown connections between them or whether they have cross-functional pathways in plant life activities. In this study, the five key candidate genes in response to cold stress were identified, including GhUVR8 (GH_A05G3668) and GhPLATZ (GH_A09G2161). Meanwhile, it was also found that GhPLATZ was subjected to low temperature treatment after transient silencing, and the silenced plants were more sensitive to low temperature and suffered more severe damage than the control plants. It is speculated that GhPLATZ is positively regulating cotton’s tolerance to low temperature.