- Research
- Open access
- Published:
Time-course transcriptome analysis reveals gene co-expression networks and transposable element responses to cold stress in cotton
BMC Genomics volume 26, Article number: 235 (2025)
Abstract
Background
Cold stress significantly challenges cotton growth and productivity, yet the genetic and molecular mechanisms underlying cold tolerance remain poorly understood.
Results
We employed RNA-seq and iterative weighted gene co-expression network analysis (WGCNA) to investigate gene and transposable element (TE) expression changes at six cold stress time points (0 h, 2 h, 4 h, 6 h, 12 h, 24 h). Thousands of differentially expressed genes (DEGs) were identified, exhibiting time-specific patterns that highlight a phase-dependent transcriptional response. While the A and D subgenomes contributed comparably to DEG numbers, numerous homeologous gene pairs showed differential expression, indicating regulatory divergence. Iterative WGCNA uncovered 125 gene co-expression modules, with some enriched in specific chromosomes or chromosomal regions, suggesting localized regulatory hotspots for cold stress response. Notably, transcription factors, including MYB73, ERF017, MYB30, and OBP1, emerged as central regulators within these modules. Analysis of 11 plant hormone-related genes revealed dynamic expression, with ethylene (ETH) and cytokinins (CK) playing significant roles in stress-responsive pathways. Furthermore, we documented over 15,000 expressed TEs, with differentially expressed TEs forming five distinct clusters. TE families, such as LTR/Copia, demonstrated significant enrichment in these expression clusters, suggesting their potential role as modulators of gene expression under cold stress.
Conclusions
These findings provide valuable insights into the complex regulatory networks underlying cold stress response in cotton, highlighting key molecular components involved in cold stress regulation. This study provides potential genetic targets for breeding strategies aimed at enhancing cold tolerance in cotton.
Background
As sessile organisms, plants must continuously adapt to environmental fluctuations, with cold stress is one of the most detrimental factors, especially as temperature extremes become more frequent due to climate change [1, 2]. Cold stress disrupts cellular homeostasis, causing membrane instability, oxidative damage, and impaired physiological processes [3,4,5]. Cotton (Gossypium spp.), a thermophilic crop native to tropical and subtropical regions, is particularly vulnerable to cold stress. Optimal growth occurs between 20-30°C, while exposure to temperatures below 15°C during critical stages, such as germination, flowering, and boll formation, can result in developmental delays, reduced fiber quality, and substantial yield losses [6,7,8]. Given the increasing unpredictability of weather patterns, understanding the physiological and molecular mechanisms underlying cotton cold stress response is essential for improving its resilience.
Plants have evolved intricate molecular pathways to mitigate the adverse effects of cold stress. Among these, the ICE-CBF-COR transcriptional cascade is one of the most extensively studied and is critical for cold response in plant species [9,10,11]. ICE1 (Inducer of CBF Expression 1) functions as an upstream regulator, activating CBF (C-repeat binding factors), a subclass of AP2/ERF family proteins. CBFs, in turn, bind to CRT/DRE elements in the promoters of cold-responsive (COR) genes, triggering the expression of genes that enhance cold tolerance. Overexpression of CBF genes has been shown to significantly enhance cold tolerance in important crops, including Brassica napus [12], potato [13], and rice [14, 15]. However, the cold stress response extends beyond the CBF pathway. Other transcription factors, such as ERF105, MYB4, and ZAT12, as well as hormone signaling pathways, play crucial roles in modulating plant response to cold stress [9, 16,17,18,19]. The integration of these diverse pathways highlights the complexity of cold tolerance mechanisms, enabling plants to balance growth and stress resilience.
While quantitative trait loci (QTL) mapping and genome-wide association studies (GWAS) have identified numerous genes associated with cold tolerance [3, 20,21,22], pinpointing specific resistance genes for complex traits like cold tolerance remains a challenge. RNA sequencing (RNA-seq) has emerged as a powerful tool for elucidating cold-responsive gene expression and unraveling the molecular mechanisms underlying cold tolerance across diverse plant species and developmental stages [23,24,25,26]. Recent advances in RNA-seq, when combined with weighted gene co-expression network analysis (WGCNA), have provided deeper insights into cold stress response in plant species [27,28,29,30]. WGCNA clusters genes with similar expression profiles into modules, enabling the detection of co-regulated gene patterns and the identification of key regulatory genes and their interactions in stress response networks. This method has been instrumental in identifying gene modules associated with cold stress and other environmental stressors. For example, in Arabidopsis, WGCNA has identified specific gene modules related to cold stress response, highlighting the role of certain transcription factors and signal transduction pathways [31, 32]. Despite these advancements, the comprehensive gene networks and specific regulatory pathways involved in cotton cold stress response remain insufficiently understood.
Transposable elements (TEs), which constitute a substantial portion of eukaryotic genomes, have been increasingly recognized for their involvement in stress responses. Under environmental stress conditions, TEs can become activated and mobilized within the genome, potentially altering the expression of nearby genes [33,34,35]. This mobilization often leads to the insertion of TEs near or within gene regulatory regions, where they can influence gene expression by creating new promoters or modifying chromatin structure, thereby activating or repressing stress-responsive genes [35,36,37]. Such changes in the genomic regulatory landscape can have profound effects on plant response to stress. Recent studies in Arabidopsis, maize, and other plants have highlighted the role of TEs in mediating responses to abiotic stresses [37,38,39]. These findings suggest that TEs not only contribute to genome plasticity but may also serve as key regulatory elements that enhance stress resilience. However, the specific roles of TEs in cold stress responses remain less understood, particularly in species like cotton. Given the significant number of TEs in the cotton genome, their activation under cold stress may play a pivotal role in regulating genes associated with cold tolerance, yet this area remains largely unexplored.
Here, we conducted a comprehensive analysis of the molecular mechanisms underlying cold stress response in cotton by integrating RNA-seq and iterative WGCNA. By identifying differentially expressed genes (DEGs) across multiple time points (2 h, 4 h, 6 h, 12 h, 24 h of cold stress), we characterized the temporal dynamics of gene expression during cold stress. Through the construction of co-expression networks, we identified gene modules and key hub genes that may play a role in cold stress regulation. Notably, we uncovered the dynamic expression patterns of TEs, suggesting their potential regulatory roles in cold stress response. Collectively, our findings provide valuable insights into the regulatory mechanisms of cold stress response and offer useful information for breeding strategies to enhance cold tolerance in cotton.
Results
Transcriptional dynamics of cotton under cold stress
To investigate the transcriptional response of cotton to cold stress, we performed RNA-seq experiments on the leaves of G. hirsutum at six different time points, including 0 h (control), 2 h, 4 h, 6 h, 12 h, and 24 h under cold stress. These time points were selected to represent key stages in the cold stress response and capture transcriptional changes associated with the stress [26, 40,41,42]. A total of 509.3 million read pairs were obtained, with an average of 28.3 million read pairs per sample (Table S1). The biological replicates showed high correlation (Pearson correlation coefficient > 0.95) and clustered together (Figure 1A, Figure S1), demonstrating the robustness and reliability of the data.
Transcriptional response of cotton to cold stress. A PCA analysis of RNA-seq samples. B Number of differentially expressed genes (DEGs) identified by comparing the cold stress samples to the control (0 h sample). C Venn diagram showing the overlap among DEGs across five cold stress time points. D GO enrichment analysis showing the enriched biological processes for DEGs at each time point under cold stress. The top 5 enriched results from each group are presented, and their enrichment patterns are compared across the groups. E Expression levels of GH_A12G2924 and GH_A13G1789 under cold stress. RNA-seq read counts were normalized using DESeq2, and results are shown as mean ± standard error (SE) (n = 3)
We identified differentially expressed genes (DEGs) by comparing the cold stress samples to the control (0 h sample). A total of 3,818-16,730 genes showed differential expression (log2(Fold Change) > 1 and adjusted p-value < 0.05) (Figure 1B). The number of DEGs increased with prolonged cold stress exposure, indicating a trend toward progressive transcriptional reprogramming in response to cold stress. Comparative analysis of these DEGs revealed that 9.59% (989/10,318) of the up-regulated genes maintained their expression across all five cold stress time points, while nearly half (47.32%, 4,882/10,318) were time-point-specific (Figure 1C). Similarly, 4.36% (483/11,072) of the down-regulated genes remained consistently down-regulated across all time points, and approximately 49.24% (5,452/11,072) were specific to individual time points (Figure 1C). This finding indicates the high temporal specificity of the transcriptional response to cold stress in cotton.
To validate the RNA-seq results, we performed quantitative real-time PCR (qRT-PCR) experiments on twelve randomly selected DEGs, confirming the RNA-seq findings and supporting the reliability of our data (Figure S2). We then performed GO and KEGG enrichment analyses to elucidate the functional implications of these DEGs (Figure 1D, Figure S3). The results revealed that these genes were involved in processes such as response to hydrogen peroxide, regulation of ethylene−activated signaling pathway, jasmonic acid mediated signaling pathway, abscisic acid metabolic process, and thiamine metabolism, which are crucial for plant cold stress response [43,44,45,46,47]. Notably, genes with sustained responses and those with time-point-specific responses were implicated in different biological processes (Figure S3B). For example, GH_A13G1789, a cotton homolog of the Arabidopsis gene NCED3, encodes 9-cis-epoxycarotenoid dioxygenase, consistently showed up-regulation at all five cold stress time points (Figure 1E). NCED3 is associated with carotenoid biosynthesis and contributes to cold stress tolerance in plants [48,49,50]. In contrast, GH_A12G2924, a cotton homolog of the Arabidopsis gene CAM5, encodes a calmodulin involved in environmental adaptation [51,52,53], showed specific up-regulation at 12 h (Figure 1E). These results highlight the dynamic transcriptional response of cotton to cold stress, where both sustained and time-specific gene expression changes are implicated in stress tolerance.
Differential subgenome responses to cold stress
As an allotetraploid species, G. hirsutum contains two distinct subgenomes, A and D. To investigate whether there are subgenomic biases in the transcriptional response to cold stress, we first examined the number of DEGs in each subgenome (Figure 2A). No significant difference was observed between the A and D subgenomes, suggesting that both subgenomes contribute similarly to the overall cold stress response. Next, we focused on homeologous gene pairs to explore potential subgenome-biased transcriptional regulation. Using reciprocal BLASTP searches, we identified 26,230 homeologous gene pairs between the A and D subgenomes. Among these, 2,386-9,218 pairs contained DEGs (Figure 2B). Notably, when comparing the fold change ratio between A-homeologues and D-homeologues, we found that a substantial proportion of homeologous gene pairs with DEGs (33.07%-44.41%) exhibited differential responses to cold stress, with a fold change ratio greater than 2 (Figure 2C). This indicates that both subgenomes are actively involved in the response to cold stress, but they display differential expression patterns.
Comparative analysis of cold-responsive genes between subgenomes. A The number of DEGs in the A subgenome and D subgenome during exposure to cold stress. B The number of homeologous gene pairs in which at least one gene is differentially expressed under cold stress. C Scatter plot showing the fold change in expression of A-homeologues and D-homeologues within homeologous gene pairs in which at least one gene is differentially expressed under cold stress. Based on the fold change ratio (>2) between A-homeologues and D-homeologues, the gene pairs were categorized into five groups: (1) A-homeologue dominant, where A-homeologues showed a greater magnitude of upregulation or downregulation compared to D-homeologues; (2) D-homeologue dominant, where D-homeologues exhibited a stronger transcriptional response; (3) static A-homeologues with responsive D-homeologues, where A-homeologues were non-differentially expressed while D-homeologues were differentially expressed; (4) static D-homeologues with responsive A-homeologues, where D-homeologues were non-differentially expressed while A-homeologues were differentially expressed; and (5) discordant responses, where A and D homeologues exhibited opposite responses to cold stress. The numbers in parentheses indicate the count of gene pairs in each category. D Representative GO terms enriched in differentially responding homeologues versus non-differentially responding homeologues. E The number of cold-responsive genes unique to the A subgenome and D subgenome at each cold stress time point
To further characterize these differential responses, we classified the differentially responding homeologous gene pairs into five distinct categories (Figure 2C): (1) A-homeologues dominant, where A-homeologues showed a greater magnitude of upregulation or downregulation compared to D-homeologues; (2) D-homeologues dominant, where D-homeologues exhibited a stronger transcriptional response; (3) static A-homeologues with responsive D-homeologues; (4) static D-homeologues with responsive A-homeologues; and (5) discordant responses, where A and D homeologues responded oppositely to cold stress. The most prevalent categories were static A-homeologues with responsive D-homeologues (369-1,483 pairs) and static D-homeologues with responsive A-homeologues (337-1,396 pairs), while the discordant responses were the least frequent (4-72 pairs) (Figure 2C). The predominance of static homeologues with responsive counterparts suggests that the A and D subgenomes may have evolved complementary roles in cold stress response. The rarity of discordant responses underscores the coordination between subgenomes. Furthermore, GO and KEGG enrichment analyses revealed that differentially responding homeologues and non-differentially responding homeologues are involved in distinct biological processes (Figure 2D, Figure S4). For example, the former were linked to processes such as cellular response to reactive oxygen species, cellular response to reactive nitrogen species, cellular response to nitric oxide, and thiamine metabolism. In contrast, the latter were enriched in processes such as negative regulation of signaling, negative regulation of signal transduction, and negative regulation of cell communication.
We also identified 3,997 subgenome-unique genes (1,951 in the A subgenome and 2,046 in the D subgenome). Among these, 155-617 genes were differentially expressed under cold stress (Figure 2E). The absence of a significant numerical bias between A and D subgenome-specific DEGs further supports the idea that both subgenomes contribute equally to the cold stress response, albeit through different gene sets and mechanisms. These findings highlight the coordinated subgenome-specific transcriptional regulation in cotton during cold stress response.
Distinct gene co-expression modules in cotton response to cold stress
To gain a deeper molecular understanding of the transcriptome changes and uncover patterns of gene expression during cold stress, we employed iterative Weighted Gene Co-expression Network Analysis (WGCNA) to construct gene co-expression networks from the RNA-seq datasets. This iterative approach provided more refined clustering compared to standard WGCNA [54], enabling us to partition 83.79% (17,492/20,877) of the DEGs into 125 distinct modules. These modules varied in size, ranging from 27 to 2,429 genes, and exhibited distinct expression patterns in response to cold stress (Figure 3A, Figure S5, Table S2). For example, module M5 comprises genes that are highly expressed under normal conditions but undergo rapid downregulation following cold exposure, maintaining low expression levels throughout the stress period. This suggests that these genes may be involved in processes that are suppressed under cold stress conditions. In contrast, module M30 contains genes that are weakly expressed in the absence of stress but are activated upon cold stress, suggesting their involvement in cold-induced pathways. Similarly, module M1 is enriched with genes that are specifically activated at the 24 h cold stress time point, suggesting their role in later stages of the cold stress response. Modules M18 and M86 are characterized by genes that are uniquely upregulated at the 2 h and 4 h time points, respectively, suggesting that these modules may represent early response genes. The distinct expression patterns observed across the modules highlight the dynamic and time-specific nature of the transcriptional reprogramming in cotton under cold stress.
Identification of gene co-expression modules during cold stress. A Example modules illustrating the progression of gene expression across different time points under cold stress. Gene expression levels are presented as TPM values. B GO enrichment analysis of representative modules. The top 10 most significantly enriched GO terms are displayed
Next, to obtain a biologically oriented description of each module, we performed GO and KEGG enrichment analyses. Among the 125 identified modules, 68 showed significant enrichment for specific GO terms, and 45 for KEGG pathways (adjusted p-value < 0.05) (Figure 3B, Table S3). For example, module M5 showed the most significant enrichment for GO terms related to the jasmonic acid-mediated signaling pathway, suggesting that these genes may be involved in modulating the plant response to stress through hormone signaling. In contrast, module M30 was most significantly enriched for GO terms associated with the response to zinc ion, suggesting a potential role in metal ion homeostasis or signaling during cold stress. Module M1 was most enriched for GO terms related to the response to chitin, suggesting its potential involvement in defense-related processes that may become critical during prolonged cold stress. Similarly, module M18 was most enriched for GO terms related to protein folding, indicating its role in maintaining protein stability and function under initial stress conditions. Module M86, enriched for GO terms associated with the response to red or far-red light, suggests a potential interaction between cold stress responses and light signaling pathways, which may be crucial for optimizing energy use and growth under stress conditions. Together, these findings highlight the intricate regulatory networks potentially employed by cotton in response to cold stress.
Chromosomal distribution of gene co-expression modules
To systematically assess whether any modules were over-represented on specific chromosomes, we performed a binomial test for each module on each chromosome. This analysis identified 18 modules with significant enrichment on particular chromosomes (p-value < 0.01) (Figure 4A). For example, module M1 showed notable over-representation of genes located on chromosomes A05 and D05, suggesting that these chromosomes disproportionately contribute to the gene content of this module. Similarly, module M4, which includes genes involved in the early cold response (Figure S5), displayed significant enrichment on chromosome D07. Notably, chromosomes A13, D02, D05, D07, and D12 were recurrently enriched across multiple modules, suggesting that these chromosomes might be involved in coordinated transcriptional programs that could be important for the stress responses. The observed diversity in the degree of enrichment across different modules and chromosomes underscores the complexity and heterogeneity of the transcriptional response, with each chromosomes contributing differently to the gene co-expression landscape.
Chromosomal distribution characteristics of identified gene co-expression modules. A Identification of chromosome-biased modules. Gene counts for each module on each chromosome were compared to the genome-wide background. A binomial test was performed for each module on each chromosome, with p-values less than 0.01 indicating significant chromosomal bias. Fold enrichment was calculated as the ratio of observed to expected values. B Scatter plot showing the significance of gene clustering for modules M1 and M81 across chromosomes. Chromosomes were divided into 500 kb windows with a 100 kb step, with each point representing a window. Gene counts for each module within each window were compared to the chromosome-wide background using a binomial test. The red dashed line indicates the p-value threshold of 0.01. C Representative example of the 21.7-22.2 Mb window on chromosome A05. Arrows indicate genes belonging to module M1, with a zoomed-in view highlighting four representative genes. RNA-seq tracks illustrate changes in signal intensity during cold stress exposure
We sought to determine whether genes within these chromosome-biased modules were clustered in specific chromosomal regions or dispersed across the chromosome. To this end, each chromosome was divided into 500 kb windows with a 100 kb step, and a binomial test was performed for each module within each window. This analysis found that 14 out of the 18 chromosome-biased modules exhibited a non-random distribution of genes across the chromosomes (p-value < 0.01) (Figure 4B, Figure S6). For example, in module M1, significant gene clustering was observed in nine windows on chromosome A05 and eight windows on chromosome D05. Further examination of these clustered regions revealed their functional significance. For example, the 21.7-22.2 Mb window on chromosome A05 contains eight genes from module M1 (Figure 4C), including GH_A05G2275, which encodes the calcium transporter protein ANN1, known for its role in cold stress tolerance in Arabidopsis by regulating Ca2+ signals after cold shock [55, 56]. Additionally, GH_A05G2252 and GH_A05G2253 encode the EXORDIUM protein, which is involved in the cold stress response in maize [57]. GH_A05G2247 encodes the FBW2 protein, a key player in ABA signaling [58]. ABA is known to regulate the plant cold response by modulating various stress pathways [59, 60]. These findings suggest that these regions might represent hotspots of co-regulated genes, highlighting the potential existence of localized regulatory networks operating within specific chromosomal regions.
Transcription factor regulation in gene co-expression modules
Given the pivotal role of transcription factors (TFs) in regulating gene expression [61, 62], we sought to identify TF genes within each gene co-expression module. This analysis identified that 123 out of the 125 modules contained TF genes (except for modules M94 and M125), with the proportion of TFs ranging from 1.30% to 24.14% (Figure S7), suggesting the widespread involvement of TFs in these gene networks. Notably, individual modules often contained multiple TF families (Figure 5A), suggesting potential coordinated regulatory interactions among these families. Comparative analysis between modules showed that the ERF family was the most prevalent, found in 65 modules, followed closely by the bHLH and MYB families, present in 57 and 56 modules, respectively (Figure 5A, Figure S8). This widespread distribution highlights the potential diverse regulatory roles of these three TF families. Further enrichment analysis, employing the hypergeometric test to compare TF family distributions within modules to the genomic background, revealed distinct patterns of TF family enrichment across modules (p-value < 0.01) (Figure 5A). For example, module M5 showed significant enrichment for the C3H, CO-like, LBD, and MYB TF families. In contrast, module M14 was predominantly enriched for ERF and NAC TF families, while module M18 exhibited exclusive enrichment for the Trihelix TF family. These findings suggest that specific TF families might be preferentially involved in regulating distinct gene modules, reflecting the complexity and specificity of transcriptional regulation under cold stress conditions in cotton.
Identification of transcription factors (TFs) in gene co-expression modules. A Bubble plot illustrating the proportion of different TF families within each module. Enrichment analysis, using the hypergeometric test to compare TF family distributions in modules to the genomic background, is indicated by red points (enriched, p-value < 0.01) and gray points (not enriched). B Bubble plot illustrating the proportion of genes regulated by identified hub TFs in each module. Larger and darker points represent higher proportions of regulated genes. C Gene regulatory network involving MYB73, ERF017, MYB30, and OBP1. FIMO was used to predict the regulatory relationships between transcription factors (TFs) and their target genes by analyzing the promoter regions (1 kb upstream of the transcription start site) of genes within each module. A gene was predicted to be regulated by a TF if the TF's binding motifs were present in its promoter region. Each node represents a gene, with different colors indicating different modules. The gray lines represent potential regulatory relationships
To further investigate the master regulators within each module, we performed FIMO analysis to predict the regulatory relationships between TFs and their target genes [63]. This analysis focused on the promoter regions (1 kb upstream of the transcription start site) of genes within each module. TFs were predicted to regulate a gene if their binding motifs were present in the gene’s promoter region. The TF with the highest number of target genes in each module was designated as the hub TF, leading to the identification of 77 hub TFs across the modules (Figure 5B). For example, in module M1, 91.89% (2232/2429) of gene promoters contained the DOF1 binding motif, identifying DOF1 as a critical regulator of this module. Similarly, in module M2, 87.25% (1006/1153) of gene promoters contained the OBP3 binding motif, and in module M3, 88.89% (824/927) contained the TCX2 binding motif. These findings suggest that these TFs may play significant roles in regulating their respective gene networks. Notably, some TFs were involved in multiple modules (Figure 5B, Figure S9). Among these, MYB73 was the most prevalent, with homologous genes found in 11 different modules, serving as the regulator in each, followed by ERF017, MYB30, and OBP1, each regulating genes in six modules. We constructed a regulatory network for these four TFs and assessed their connectivity degree within the network, which ranged from 379 to 2836 (Figure 5C). This widespread involvement highlights the importance of these TFs in the transcriptional response to cold stress, potentially orchestrating gene expression across multiple co-expression modules.
Differential expression of hormone-related genes in cotton under cold stress
Previous studies in Arabidopsis have highlighted the key roles of plant hormones such as salicylic acid (SA), ethylene (ETH), jasmonic acid (JA), and strigolactones (SLs) in mediating plant responses to cold stress [64,65,66]. To investigate the involvement of hormone-related genes in cotton during cold stress, we examined the differential expression of genes associated with 11 plant hormones, including SA, ETH, JA, SLs, auxin (AUX), cytokinin (CK), abscisic acid (ABA), gibberellin (GA), brassinosteroids (BRs), melatonin (MT), and polypeptide (PEP) hormones. These genes were identified using lists from the Plant Hormone Gene Database (PHGD) [67]. We found that a substantial proportion of hormone-related genes exhibit differential expression during cold stress, with the percentages ranging from 4.82% to 39.39% (Figure 6A, Figure S10). The progressive increase in DEGs over time corresponds with the overall genomic trend of increasing gene expression changes under cold stress (Figure 1B). Notably, the proportion of DEGs varied significantly among different hormones, with the highest proportions observed for ETH and CK (Figure 6A). Importantly, the number of upregulated genes in ETH and CK significantly exceeded the number of downregulated genes, suggesting that these two hormones might play a prominent role in orchestrating the activation of key stress-responsive pathways in cotton. In contrast, hormones such as PEP and SLs exhibited a moderate response, with consistently lower proportions of DEGs throughout the cold stress exposure, suggesting a potentially limited role in the cold stress response. Interestingly, no DEGs were observed for MT at 2 h of cold exposure, but DEGs emerged at 4 h and beyond (Figure S10), suggesting a delayed response of MT-related genes to cold stress.
Expression dynamics of hormone-related genes in cotton under cold stress. A Proportions of DEGs in hormone-related genes across different time points under cold stress. B Heatmap showing the expression dynamics of hormone-related DEGs over the course of cold stress. Gene expression values are presented as TPM and normalized by row-wise z-scores
To further understand the hormonal regulation underlying cotton response to cold stress, we examined the expression patterns of hormone-related genes categorized by their roles in metabolism, signaling, synthesis, transport, and receptor activity. Distinct expression patterns were observed among hormone-related genes (Figure 6B, Figure S11). For example, in CK-related genes, receptor genes were downregulated under cold stress, while signaling genes were upregulated, particularly at the 24 h of exposure. In contrast, ETH receptor and signaling genes both exhibited upregulation under cold stress. JA synthesis genes were specifically upregulated at 24 h of cold exposure, whereas AUX synthesis genes were suppressed at the same time point (Figure S11). ABA transport genes showed high expression under normal conditions but were inhibited by cold stress, while ETH transport genes exhibited low expression in the absence of stress. Together, these findings highlight the complexity and dynamism of the hormonal network in the cold stress response.
Dynamic TE expression in cotton during cold stress
As transposable elements (TEs) are highly responsive to environmental changes [68, 69], we investigated their expression during cold stress in cotton. To evaluate TE expression, we employed the Telescope tool [70], which enabled us to analyze individual TE loci by quantifying the reads mapping to these loci. Given the pervasive nature of TEs throughout the genome and their potential integration into gene transcripts, distinguishing whether a TE-mapping read is derived from a gene or a TE may present challenges. To address this, we categorized TEs based on their positional relationship to genes, focusing on reads that mapped to TEs without overlapping genes for further analysis. This approach revealed that 15,352 TEs were expressed (with TE-mapping read >10) in at least one sample, suggesting the significant involvement of TEs in the transcriptional response to cold stress in cotton.
We next conducted differential expression analysis to identify TEs that are transcriptionally regulated during cold stress in cotton. By comparing cold stress samples to control samples, we observed that 1,421-7,075 TEs exhibited differential expression (log2(Fold Change) > 1 and adjusted p-value < 0.05) (Figure 7A). qRT-PCR experiments on eight randomly selected differentially expressed TEs confirmed the RNA-seq findings (Figure S12). Notably, the number of differentially expressed TEs increased with prolonged cold stress exposure, with a greater prevalence of up-regulated TEs compared to down-regulated ones. Further analysis revealed that 10.39% (512/4,928) of the up-regulated TEs maintained elevated expression across all five cold stress time points, while nearly one-third (31.96%, 1,575/4,928) exhibited time-point-specific up-regulation (Figure 7B). Similarly, 6.96% (250/3,592) of the down-regulated TEs remained consistently down-regulated across all time points, with approximately 42.34% (1,521/3,592) showing time-point-specific down-regulation. To further elucidate their response patterns, we performed a clustering analysis of TE expression profiles, which identified five distinct temporal expression patterns (Figure 7C). Specifically, Cluster 1 (1,506 TEs, 17.77%) exhibited significant up-regulation starting at 12 h of cold stress. Cluster 2 (679 TEs, 8.01%) was notably activated at 4 h but subsequently suppressed. Cluster 3 (2,853 TEs, 33.67%) remained inactive until being activated at 24 h of cold stress. Cluster 4 (1,988 TEs, 23.46%) was expressed under normal conditions but was immediately suppressed upon cold stress exposure. In contrast, Cluster 5 (1,447 TEs, 17.08%) showed expression under normal conditions that progressively suppressed with increasing cold stress duration. These findings suggest that TE expression is tightly regulated in a time-dependent manner during the cotton cold stress response.
TEs response to cold stress in cotton. A Number of differentially expressed transposable elements (TEs) identified by comparing the cold stress samples to the control (0 h sample). B Venn diagram showing the overlap among differentially expressed TEs across five cold stress time points. C Clustering analysis of differentially expressed TEs. The heatmap displays normalized read counts for each TE locus across all samples, following z-score transformation. Clusters were identified using k-means clustering. D TE family-specific enrichment analysis per expression clusters. Only significantly enriched TE families are shown (p-value < 0.01, Fisher’s exact test). The odds ratio represents the degree of enrichment, with larger values indicating stronger enrichment and smaller values suggesting under-representation. Both p-values and odds ratios were calculated using the fisher.test function in R. Class I elements (retrotransposons): LINE, LTR/Copia, LTR/Gypsy, LTR/unknown. Class II elements (DNA transposons): hAT (DTA), CACTA (DTC), PIF–Harbinger (DTH), Mutator (DTM), Tc1–Mariner (DTT), Helitron
TEs are typically classified into two major classes: retrotransposons (Class I) and DNA transposons (Class II), each comprising various families. To uncover the roles of different TE families during cold stress, we performed an enrichment analysis to detect over- or under-represented TE families within each expression cluster (Fisher's exact test). Distinct enrichment patterns were observed across the clusters (Figure 7D). For example, the retrotransposon LTR/Copia family was significantly enriched in all five clusters (p-value < 0.01), suggesting its extensive involvement in the cold stress response. LINE elements were enriched in four clusters (excluding Cluster 2), with varying degrees of enrichment as reflected by odds ratio values. In contrast, the LTR/Gypsy family was conspicuously absent in Clusters 1, 2, and 3, suggesting its limited role in the later stages of the cold stress response. Among DNA transposons, the DNA/DTA family was consistently enriched across all five clusters, whereas the DNA/DTT family showed no significant enrichment in any of the clusters. Notably, the dynamic expression patterns of differentially expressed TEs and their closest neighboring genes exhibited a moderate positive correlation (mean Pearson correlation coefficients ranging from 0.153 to 0.430) (Figure S13), suggesting a potential regulatory relationship between TEs and adjacent genes under cold stress. Together, these results suggest that different TE families may play distinct roles in cotton response to cold stress.
Discussion
Cold stress presents a significant challenge to cotton, adversely affecting its growth and productivity. Despite extensiveresearch into cold tolerance mechanisms, a comprehensive understanding of the genetic and molecular networks that govern cold tolerance in cotton remains limited. Previous studies have identified several transcription factors and pathways involved in cold responses in plant species, but the specific regulators and gene networks in cotton are not fully characterized. Furthermore, while transposable elements (TEs) are increasingly recognized for their role in environmental stress adaptation, their contribution to cold tolerance in cotton remains largely unexplored. In this study, using RNA-seq and iterative WGCNA, we captured the dynamic changes in gene and TE expression across multiple time points during cold stress. This approach allowed us to unveil intricate regulatory networks crucial for cotton response to cold stress.
We identified thousands of differentially expressed genes (DEGs) across multiple cold stress time points, with a significant proportion showing time-specific expression patterns (Figure 1C). This temporal specificity suggests that cotton employs a finely tuned, phase-dependent transcriptional response to cold stress, similar to observations in Arabidopsis and rice. In these species, distinct waves of gene expression are triggered at different stages of cold exposure, coordinating immediate stress perception with later tolerance [71,72,73,74]. Importantly, our findings highlight the differential contribution of the A and D subgenomes to the cold stress response in G. hirsutum (Figure 2). While the overall number of DEGs was comparable between the two subgenomes, 33.07%-44.41% of homeologous gene pairs with DEGs exhibited subgenome-biased expression patterns, suggesting that the A and D subgenomes respond differentially to cold stress. This observation is consistent with reports from other polyploid crops, such as wheat and Brassica napus, where subgenomes exhibit distinct transcriptional responses to environmental stress [75,76,77]. Notably, the observation that one homeologue remains relatively static while the other exhibits dynamic responsiveness to stress is widespread (Figure 2C), suggesting a divergence in homologous gene regulation. This divergence likely reflects a strategic division of labor between subgenomes. Indeed, such regulatory divergence has been well documented in polyploid species [78,79,80], indicating that this separation of roles may contribute to enhanced adaptability and resilience under stress conditions. Furthermore, we identified 3,997 genes unique to the subgenome, with no significant numerical bias between A and D subgenome-specific differentially expressed genes (Figure 2E). These subgenome-specific genes likely participate in specialized regulatory pathways that complement the stress responses mediated by homeologous genes. This dual strategy, integrating both homeologous and subgenome-specific genes, likely contributes to the evolutionary advantage of polyploidy, enhancing the plant capacity for stress tolerance.
Compared to traditional WGCNA, iterative WGCNA enabled us to identify more functionally cohesive modules [54], offering a more nuanced view of the transcriptional response to cold stress. This method revealed distinct temporal dynamics in gene expression, reflecting plant adaptive strategies at various stages of stress exposure. In total, we identified 125 modules, each representing specific aspects of the stress response (Figure 3A, Figure S5). The clustering of early-response genes into specific modules, such as M18 and M86, indicates a rapid transcriptional shift shortly after cold exposure. These modules are enriched with genes involved in protein folding, sesquiterpene metabolism, and light responses, which have been implicated in cold stress responses [81,82,83,84]. These genes likely serve as key regulators, orchestrating the initial activation of cold-responsive pathways and preparing the plant for subsequent stress responses. In contrast, modules active at later time points, such as M1, which is enriched in genes active at 24 h of cold stress, are involved in processes such as the response to chitin (Figure 3). The activation of chitin-responsive pathways suggests cross-talk between abiotic and biotic stress responses, as chitin-related signaling is typically associated with pathogen defense but may also play a role in cold stress tolerance [85,86,87]. Interestingly, our analysis identified significant enrichment of some modules on specific chromosomes or chromosomal regions. For example, module M1 was notably over-represented on chromosomes A05 and D05 (Figure 4A), suggesting that these chromosomes play a crucial role in stress response . Additionally, the non-random distribution of modules, particularly their clustering in specific chromosomal regions (Figure 4B), points to the presence of localized regulatory networks. Such clustering of stress-responsive genes likely reflects their co-regulation, enhancing their ability to respond efficiently to environmental stresses [88].
A recent study on cotton cold stress response reports the physiological responses and transcriptomic changes in cotton leaves under varying cold stress temperatures (15°C, 10°C, and 4°C) over 24 hours, highlighting the role of transcription factors such as GhPLATZ, which acts as a positive regulator of cold tolerance in cotton [89]. Our analysis similarly highlights the significant role of transcription factors (TFs) in regulating gene co-expression modules associated with cold stress in cotton. TF families such as ERF, bHLH, and MYB are notably represented within these modules (Figure 5A). This prominence is consistent with their established roles in stress responses, including cold stress, where they regulate various aspects of stress tolerance by modulating the expression of downstream genes [2, 90, 91]. Moreover, the identification of hub TFs within each gene co-expression module further clarifies the key regulatory nodes driving these networks. For example, DOF1 in module M1, OBP3 in module M2, and TCX2 in module M3 were pinpointed as central regulators due to the presence of their binding motifs in the majority of gene promoters within these modules (Figure 5B). This suggests that these TFs are critical in orchestrating gene expression during cold stress, serving as pivotal nodes within their respective networks. Notably, the extensive involvement of TFs such as MYB73, ERF017, MYB30, and OBP1 across multiple modules indicates their multifaceted functions (Figure 5B), coordinating diverse stress response pathways and contributing to a robust regulatory framework for managing cold stress. Indeed, previous studies have demonstrated the involvement of several of these TFs, including DOF1, MYB73, and MYB30, in cold stress responses in plant species. For example, overexpression of Dof1 in upland cotton not only altered seed oil composition but also significantly enhanced tolerance to both salt and cold stress [92]. MYB73, a CBF-independent first-wave transcription factor, is activated during cold acclimation [73, 93], while MYB30 has been identified as a key regulator of cuticle synthesis in response to cold stress [94, 95].
The dynamic expression of TEs during cold stress in cotton reveals an intriguing layer of gene regulation that may contribute to stress tolerance. Our results show that over 15,000 TEs are expressed in response to cold, with a notable increase in the number of differentially expressed TEs as stress duration progresses (Figure 7A). This trend, particularly the prevalence of upregulated TEs, suggests that TEs could serve as important players in shaping the transcriptional landscape during stress response, a phenomenon previously reported in other species [96,97,98]. This raises the possibility that TEs act as regulatory elements or even stress sensors, modulating nearby gene expression in response to environmental cues. Such a mechanism could be particularly important during stress conditions, where rapid shifts in gene expression are essential for maintaining cellular homeostasis. We identified five distinct temporal expression patterns of TEs throughout the cold stress response (Figure 7C). This temporal stratification suggests that TEs may be activated at specific stages of stress, contributing to a phase-specific regulatory network. Notably, the enrichment of specific TE families, such as LTR/Copia, across distinct clusters underscores their potentially multifaceted roles in cotton response to cold stress (Figure 7D). LTR/Copia elements are known for their ability to influence gene expression by providing regulatory sequences [99,100,101]. Their recurrent presence in multiple clusters suggests a broader role in initiating or amplifying stress-responsive transcriptional programs. In contrast, the absence of significant enrichment of LTR/Gypsy elements in these clusters observed during cold stress suggests that these TEs may not play a prominent role in the immediate transcriptional adjustments necessary for rapid stress responses (Figure 7D). LTR/Gypsy elements are often concentrated in heterochromatic regions, particularly around centromeres, where they primarily contribute to maintaining genome stability [33, 102]. This spatial organization likely limits their involvement in the dynamic gene regulation needed during environmental stress.
Tissue-specific responses to environmental stress are well-documented, with different tissues exhibiting distinct transcriptional profiles under stress conditions [90, 103,104,105]. Our study focused on leaves, which, due to their direct exposure to the external environment, are highly sensitive to temperature fluctuations, making them ideal for capturing rapid transcriptional changes in response to cold stress. Other tissues, such as roots, may respond differently under cold stress, given their unique physiological roles and relatively insulated environment within the soil [105, 106]. For example, in Norway spruce, a bHLH101-like transcription factor was found to be co-expressed with a root-specific subset of genes within the gene-regulatory network under cold stress, highlighting the importance of tissue-specific responses in cold hardiness [106]. Furthermore, recent advancements in single-cell sequencing have revealed substantial heterogeneity in stress responses even among different cell types within the same tissue [107,108,109], emphasizing the complexity and specificity of stress tolerance mechanisms. Future studies integrating transcriptomic analyses across multiple tissues and employing single-cell resolution approaches will be essential to unravel the comprehensive regulatory networks underlying cotton response to cold stress.
In summary, our study proposes a potential mechanism by which cotton responds and acclimates to cold stress. First, cold stress is likely sensed by receptors, leading to the activation of signal transduction pathways, as evidenced by the upregulation of receptor kinase-related genes (Figure S5, Table S3). These pathways may initiate cold-responsive signaling cascades. Second, transcriptional regulation appears to play a pivotal role, with significant enrichment of transcription factors such as ERF, bHLH, and MYB in the co-expression gene modules (Figure 5A), which may modulate downstream stress-responsive genes to enhance cold tolerance. Third, metabolic reprogramming appears to occur, with altered expression of genes involved in carbon and nitrogen metabolism, including those related to sucrose and amino acid metabolism (Table S3), which may provide energy and osmotic protective compounds. Additionally, hormonal signaling pathways, such as abscisic acid, jasmonic acid, and ethylene, which are known to be related to cold acclimation [9, 110, 111], appear to be differentially regulated under cold stress (Figure 6, Figure S10-S11). These findings collectively suggest a comprehensive framework of the signaling, transcriptional, metabolic, and hormonal adjustments that may underpin cotton response and acclimation to cold stress.
Conclusions
This study presents a comprehensive exploration of the molecular mechanisms potentially underlying cold stress response in cotton. By integrating RNA-seq with iterative WGCNA, we identified co-expression networks and characterized key regulators, including MYB73, ERF017, MYB30, and OBP1, which may coordinate multiple stress response networks. Furthermore, we observed the dynamic expression of TEs during the stress response, particularly highlighting the potential contribution of LTR/Copia elements in modulating gene expression. These findings offer valuable insights into the regulatory mechanisms of cold stress response and may inform breeding strategies aimed at developing cold-tolerant cotton varieties.
Materials and methods
Plant materials and growth conditions
Cotton (G. hirsutum, TM-1) seeds were germinated and cultivated in a controlled greenhouse environment (13/11 h light/dark, 28/26℃ light/dark, 60% humidity). Yan Dai and Jinlei Han undertook the formal identification of the plant material used in the study, and no voucher specimens were collected. At the three-leaf and one-heart stage, seedlings were subjected to cold stress by transferring them to an incubator set at 4°C for varying durations (0 h, 2 h, 4 h, 6 h, 12 h, and 24 h). Leaf tissues were harvested at each time point, immediately frozen in liquid nitrogen, and stored at -80°C for further analysis.
RNA-seq and data analysis
Total RNA was extracted from leaf samples (three biological replicates per time point) using the Omega Plant RNA kit (Omega Bio-tek, Cat. no. R6827-01) according to the manufacturer's protocol. RNA-seq libraries were prepared using the Illumina TruSeq RNA Kit (NEB, Cat. no. E7530) and sequenced on an Illumina NovaSeq 6000 system, generating 150 bp paired-end reads. Clean reads were mapped to the G. hirsutum reference genome (ZJU-v2.1) using Tophat2 v.2.1.1 [112] with default parameters. Genome sequence and annotation files were downloaded from CottonGen (https://www.cottongen.org/). Gene expression levels were calculated using StringTie v.2.1.5 [113] as transcripts per million (TPM). For differential expression analysis, read counts per gene were determined using FeatureCounts v.2.0.1 [114] and analyzed with DEseq2 [115], with significant differentially expressed genes (DEGs) identified based on an adjusted p-value < 0.05 and |log2(fold change)| > 1. Principal component analysis (PCA) was performed using the DESeq2 plotPCA function. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was conducted using clusterProfiler [116], with significance determined by adjusted p-value < 0.05.
Identification of homeologous gene pairs and subgenome-unique genes
Homeologous gene pairs were identified as previously described [117, 118]. Protein sequence from the A and D subgenomes were used as queries in reciprocal BLAST searches. Gene pairs exhibiting the best reciprocal BLAST hits, with criteria of E-value < 1e-5, sequence coverage ≥ 50%, and sequence identity ≥ 50%, were classified as homeologous genes. Conversely, genes with little or no sequence similarities between the A and D subgenomes (coverage < 50% or identity < 50%) were categorized as subgenome-unique genes. A total of 26,230 homeologous gene pairs and 3,997 subgenome-unique genes (1,951 A subgenome-unique genes and 2,046 D subgenome-unique genes) were identified.
Co-expression network identification and chromosomal enrichment analysis
Co-expression networks were constructed using iterative weighted gene co-expression network analysis (WGCNA) (https://github.com/cstoeckert/iterativeWGCNA) with default parameters. The resulting co-expression networks, termed modules, were analyzed for chromosomal distribution. To identify whether module genes were significantly enriched on specific chromosomes, the number of genes within a given module located on a specific chromosome was compared to the total number of genes on that chromosome across the entire genome, using a binomial test (R function: binom.test). A p-value < 0.01 was considered statistically significant for enrichment. Additionally, to evaluate the enrichment of module-associated genes within chromosomal regions, each chromosome was divided into 500 kb windows with a 100 kb step, and a binomial test was applied to each window using the same statistical method.
Transcription factor identification and enrichment analysis
Transcription factors (TFs) were identified by performing BLAST searches of cotton protein sequences against the PlantTFDB database (http://planttfdb.gao-lab.org/). To assess the overrepresentation of TF families within specific co-expression modules, we applied a hypergeometric test (R function: phyper). This test compared the observed distribution of TF families within each module to their expected distribution across the entire genome. Enrichment of TF families was considered significant if the p-value was less than 0.01, highlighting TF families disproportionately associated with specific co-expression modules.
Defining predicted binding sites for TFs
TF DNA binding motifs were obtained from PlantTFDB database. Potential TF binding sites were determined by scanning for these motifs within the regions of interest using FIMO v.5.4.1 [119]. Candidate target genes for each TF were determined based on the presence of the TF binding motif in their promoter regions, defined as 1 kb upstream of the transcription start site (TSS).
Quantitative RT‐PCR analysis
RNA was first extracted from the samples and reverse transcribed using the StarScript II RT Mix with gDNA Remover (Genstar, Cat#A224-10). qRT-PCR was performed with ChamQ SYBR qPCR Master Mix (Vazyme, Q3331-02) on a CFX Connect Real-Time PCR Detection System (Bio-Rad). The 2−ΔΔCT method was used for normalization and calculation of relative expression levels with three replicates. Cotton ubiquitin 14 (UBQ14) was used as an internal control. Primer sequences are detailed in Table S4.
TE annotation and differential expression analysis
Transposable elements (TEs) in the G. hirsutum genomes were annotated using the EDTA v.1.9.6 pipeline [120]. The resulting annotation table was parsed and hierarchically classified into two major classes: retrotransposons (including Gypsy, Copia, etc.) and DNA transposons (including hAT, CACTA, etc.). For differential expression analysis, clean reads were mapped to the reference genome using STAR v.2.5.2b [121] with the following parameters: ‐‐chimSegmentMin 10 ‐‐winAnchorMultimapNmax 200 ‐‐outFilterMultimapNmax 100. The alignment files were then sorted and indexed using Samtools v1.9 [122]. Read counts for each TE were obtained using the Telescope tool [70]. Differential expression analysis of TE loci was conducted with DESeq2, applying a threshold of adjusted p-value < 0.05 and |log2(fold change)| > 1. TEs with fewer than 10 reads in any sample were excluded from the analysis.
TE clustering and enrichment analysis
The normalized count matrix for differentially expressed TE loci was standardized using z-score transformation. The resulting matrix was then clustered using k-means clustering with the R function kmeans, applying the following parameters: iter.max = 500, nstart = 50 and algorithm = “Lloyd”. The number of clusters was determined to be five based on visual inspection. The clustered matrix was visualized using a heatmap generated by the R function pheatmap. To assess the enrichment of TE families within each cluster, Fisher's exact test was conducted using the R function fisher.test.
Data visualization
For visualization, the filtered, sorted, and indexed BAM files were converted to bigwig format using the bamCoverage function in deepTools v.3.1.3 [123] with a bin size of 10 bp and RPKM normalization. Genome browser images were generated using the Integrative Genomics Viewer (IGV) v.2.3.92 with bigwig files processed as described above.
Data availability
The RNA–seq data described in this work have been deposited to the Genome Sequence Archive (GSA) database (http://gsa.big.ac.cn/) under the accession number CRA019230.
References
Kim JS, Kidokoro S, Yamaguchi-Shinozaki K, Shinozaki K. Regulatory networks in plant responses to drought and cold stress. Plant Physiol. 2024;195(1):170–89.
Kidokoro S, Shinozaki K, Yamaguchi-Shinozaki K. Transcriptional regulatory network of plant cold-stress responses. Trends in plant science. 2022;27(9):922–35.
Zhou X, Muhammad I, Lan H, Xia C. Recent Advances in the Analysis of Cold Tolerance in Maize. Front Plant Sci. 2022;13:866034.
Xu Z, Zhang J, Wang X, Essemine J, Jin J, Qu M, Xiang Y, Chen W. Cold-induced inhibition of photosynthesis-related genes integrated by a TOP6 complex in rice mesophyll cells. Nucleic Acids Res. 2023;51(4):1823–42.
Tiwari M, Kumar R, Subramanian S, Doherty CJ, Jagadish SVK. Auxin-cytokinin interplay shapes root functionality under low-temperature stress. Trends in plant science. 2023;28(4):447–59.
Li ZB, Zeng XY, Xu JW, Zhao RH, Wei YN. Transcriptomic profiling of cotton Gossypium hirsutum challenged with low-temperature gradients stress. Scientific data. 2019;6(1):197.
Dhaliwal LK, Shim J, Auld D, Angeles-Shim RB. Fatty acid unsaturation improves germination of upland cotton (Gossypium hirsutum) under cold stress. Front Plant Sci. 2024;15:1286908.
Kaur Dhaliwal L, Gannaban RB, Shrestha A, Shim J, Kaur Mangat P, Singleton JJ, Angeles-Shim RB. Integrated morpho-biochemical and transcriptome analyses reveal multidimensional response of upland cotton (Gossypium hirsutum L.) to low temperature stress during seedling establishment. Plant-environment interactions (Hoboken, NJ). 2021;2(6):290–302.
Shi Y, Ding Y, Yang S. Molecular regulation of CBF signaling in cold acclimation. Trends in plant science. 2018;23(7):623–37.
Hwarari D, Guan Y, Ahmad B, Movahedi A, Min T, Hao Z, Lu Y, Chen J, Yang L. ICE-CBF-COR Signaling Cascade and Its Regulation in Plants Responding to Cold Stress. Int J Mol Sci. 2022;23(3):1549.
Gusain S, Joshi S, Joshi R. Sensing, signalling, and regulatory mechanism of cold-stress tolerance in plants. Plant physiology and biochemistry : PPB. 2023;197:107646.
Jaglo KR, Kleff S, Amundsen KL, Zhang X, Haake V, Zhang JZ, Deits T, Thomashow MF. Components of the Arabidopsis C-repeat/dehydration-responsive element binding factor cold-response pathway are conserved in Brassica napus and other plant species. Plant Physiol. 2001;127(3):910–7.
Pino MT, Skinner JS, Park EJ, Jeknić Z, Hayes PM, Thomashow MF, Chen TH. Use of a stress inducible promoter to drive ectopic AtCBF expression improves potato freezing tolerance while minimizing negative effects on tuber yield. Plant Biotechnol J. 2007;5(5):591–604.
Ito Y, Katsura K, Maruyama K, Taji T, Kobayashi M, Seki M, Shinozaki K, Yamaguchi-Shinozaki K. Functional analysis of rice DREB1/CBF-type transcription factors involved in cold-responsive gene expression in transgenic rice. Plant & cell physiology. 2006;47(1):141–53.
Byun MY, Cui LH, Lee J, Park H, Lee A, Kim WT, Lee H. Identification of Rice Genes Associated With Enhanced Cold Tolerance by Comparative Transcriptome Analysis With Two Transgenic Rice Plants Overexpressing DaCBF4 or DaCBF7, Isolated From Antarctic Flowering Plant Deschampsia antarctica. Front Plant Sci. 2018;9:601.
Illgen S, Zintl S, Zuther E, Hincha DK, Schmülling T. Characterisation of the ERF102 to ERF105 genes of Arabidopsis thaliana and their role in the response to cold stress. Plant Mol Biol. 2020;103(3):303–20.
Vogel JT, Zarka DG, Van Buskirk HA, Fowler SG, Thomashow MF. Roles of the CBF2 and ZAT12 transcription factors in configuring the low temperature transcriptome of Arabidopsis. Plant J. 2005;41(2):195–211.
Vannini C, Locatelli F, Bracale M, Magnani E, Marsoni M, Osnato M, Mattana M, Baldoni E, Coraggio I. Overexpression of the rice Osmyb4 gene increases chilling and freezing tolerance of Arabidopsis thaliana plants. Plant J. 2004;37(1):115–27.
Waadt R, Seller CA, Hsu PK, Takahashi Y, Munemasa S, Schroeder JI. Plant hormone regulation of abiotic stress responses. Nature reviews Molecular cell biology. 2022;23(10):680–94.
Li C, Liu J, Bian J, Jin T, Zou B, Liu S, Zhang X, Wang P, Tan J, Wu G, et al. Identification of cold tolerance QTLs at the bud burst stage in 211 rice landraces by GWAS. BMC Plant Biol. 2021;21(1):542.
Schläppi MR, Jessel AR, Jackson AK, Phan H, Jia MH, Edwards JD, Eizenga GC. Navigating rice seedling cold resilience: QTL mapping in two inbred line populations and the search for genes. Front Plant Sci. 2023;14:1303651.
Chen Y, Liu Z, Han D, Yang Q, Li C, Shi X, Zhang M, Yang C, Qiu L, Jia H, et al. Cold tolerance SNPs and candidate gene mining in the soybean germination stage based on genome-wide association analysis. TAG Theoretical and applied genetics Theoretische und angewandte Genetik. 2024;137(8):178.
Zhong Y, Luo Y, Sun J, Qin X, Gan P, Zhou Z, Qian Y, Zhao R, Zhao Z, Cai W, et al. Pan-transcriptomic analysis reveals alternative splicing control of cold tolerance in rice. The Plant cell. 2024;36(6):2117–39.
Yu G, Ye F, Zhang X, Cai J, Zhu W, Zhang H, Chen S, Han J, Wang K. Characterization of open chromatin in response to cold reveals transcription factor association with preferred binding distances in cassava. Industrial Crops and Products. 2023;202:117055.
Cai P, Lan Y, Gong F, Li C, Xia F, Li Y, Fang C. Comparative physiology and transcriptome response patterns in cold-tolerant and cold-sensitive varieties of Solanum melongena. BMC Plant Biol. 2024;24(1):256.
Wang S, Shen Y, Deng D, Guo L, Zhang Y, Nie Y, Du Y, Zhao X, Ye X, Huang J, et al. Orthogroup and phylotranscriptomic analyses identify transcription factors involved in the plant cold response: A case study of Arabidopsis BBX29. Plant communications. 2023;4(6):100684.
Cheng M, Pan Z, Cui K, Zheng J, Luo X, Chen Y, Yang T, Wang H, Li X, Zhou Y, et al. RNA sequencing and weighted gene co-expression network analysis uncover the hub genes controlling cold tolerance in Helictotrichon virescens seedlings. Front Plant Sci. 2022;13:938859.
Luo Z, Zhou Z, Li Y, Tao S, Hu ZR, Yang JS, Cheng X, Hu R, Zhang W. Transcriptome-based gene regulatory network analyses of differential cold tolerance of two tobacco cultivars. BMC Plant Biol. 2022;22(1):369.
Wang X, Liu Y, Han Z, Chen Y, Huai D, Kang Y, Wang Z, Yan L, Jiang H, Lei Y, et al. Integrated Transcriptomics and Metabolomics Analysis Reveal Key Metabolism Pathways Contributing to Cold Tolerance in Peanut. Front Plant Sci. 2021;12:752474.
Zheng L, Li B, Zhang G, Zhou Y, Gao F. Jasmonate enhances cold acclimation in jojoba by promoting flavonol synthesis. Hortic Res. 2024;11(7):uhae125.
Zhang Z, Zhang T, Ma L. Analysis of basic pentacysteine6 transcription factor involved in abiotic stress response in Arabidopsis thaliana. Frontiers in genetics. 2023;14:1097381.
Li F, Hu Q, Chen F, Jiang JF. Transcriptome analysis reveals Vernalization is independent of cold acclimation in Arabidopsis. BMC genomics. 2021;22(1):462.
Baduel P, Leduque B, Ignace A, Gy I, Gil J Jr, Loudet O, Colot V, Quadrana L. Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana. Genome Biol. 2021;22(1):138.
Lisch D. Epigenetic regulation of transposable elements in plants. Annu Rev Plant Biol. 2009;60:43–66.
Negi P, Rai AN, Suprasanna P. Moving through the Stressed Genome: Emerging Regulatory Roles for Transposons in Plant Stress Response. Front Plant Sci. 2016;7:1448.
Quesneville H. Twenty years of transposable element analysis in the Arabidopsis thaliana genome. Mobile DNA. 2020;11:28.
Roquis D, Robertson M, Yu L, Thieme M, Julkowska M, Bucher E. Genomic impact of stress-induced transposable element mobility in Arabidopsis. Nucleic Acids Res. 2021;49(18):10431–47.
Makarevitch I, Waters AJ, West PT, Stitzer M, Hirsch CN, Ross-Ibarra J, Springer NM. Transposable elements contribute to activation of maize genes in response to abiotic stress. PLoS Genet. 2015;11(1):e1004915.
Yasuda K, Ito M, Sugita T, Tsukiyama T, Saito H, Naito K, Teraishi M, Tanisaka T, Okumoto Y. Utilization of transposable element mPing as a novel genetic tool for modification of the stress response in rice. Molecular breeding : new strategies in plant improvement. 2013;32(3):505–16.
Prerostova S, Černý M, Dobrev PI, Motyka V, Hluskova L, Zupkova B, Gaudinova A, Knirsch V, Janda T, Brzobohatý B, et al. Light Regulates the Cytokinin-Dependent Cold Stress Responses in Arabidopsis. Front Plant Sci. 2020;11:608711.
Pagter M, Alpers J, Erban A, Kopka J, Zuther E, Hincha DK. Rapid transcriptional and metabolic regulation of the deacclimation process in cold acclimated Arabidopsis thaliana. BMC genomics. 2017;18(1):731.
Wu H, Zheng L, Qanmber G, Guo M, Wang Z, Yang Z. Response of phytohormone mediated plant homeodomain (PHD) family to abiotic stress in upland cotton (Gossypium hirsutum spp.). BMC Plant Biol. 2021;21(1):13.
Wang L, Zhao R, Zheng Y, Chen L, Li R, Ma J, Hong X, Ma P, Sheng J, Shen L. SlMAPK1/2/3 and Antioxidant Enzymes Are Associated with H(2)O(2)-Induced Chilling Tolerance in Tomato Plants. Journal of agricultural and food chemistry. 2017;65(32):6812–20.
Ren C, Fan P, Li S, Liang Z. Advances in understanding cold tolerance in grapevine. Plant Physiol. 2023;192(3):1733–46.
An JP, Liu ZY, Zhang XW, Wang DR, Zeng F, You CX, Han Y. Brassinosteroid signaling regulator BIM1 integrates brassinolide and jasmonic acid signaling during cold tolerance in apple. Plant Physiol. 2023;193(2):1652–74.
An S, Liu Y, Sang K, Wang T, Yu J, Zhou Y, Xia X. Brassinosteroid signaling positively regulates abscisic acid biosynthesis in response to chilling stress in tomato. Journal of integrative plant biology. 2023;65(1):10–24.
Konecny T, Nikoghosyan M, Binder H. Machine learning extracts marks of thiamine’s role in cold acclimation in the transcriptome of Vitis vinifera. Front Plant Sci. 2023;14:1303542.
Esmaeili N, Shen G, Zhang H. Genetic manipulation for abiotic stress resistance traits in crops. Front Plant Sci. 2022;13:1011985.
Luo X, Bai X, Zhu D, Li Y, Ji W, Cai H, Wu J, Liu B, Zhu Y. GsZFP1, a new Cys2/His2-type zinc-finger protein, is a positive regulator of plant tolerance to cold and drought stress. Planta. 2012;235(6):1141–55.
Chao L, Kim Y, Gilmour SJ, Thomashow MF. Temperature modulation of CAMTA3 gene induction activity is mediated through the DNA binding domain. Plant J. 2022;112(1):235–48.
Kang X, Zhao L, Liu X. Calcium Signaling and the Response to Heat Shock in Crop Plants. Int J Mol Sci. 2023;25(1):324.
Al-Quraan NA, Locy RD, Singh NK. Expression of calmodulin genes in wild type and calmodulin mutants of Arabidopsis thaliana under heat stress. Plant physiology and biochemistry : PPB. 2010;48(8):697–702.
Dal Santo S, Palliotti A, Zenoni S, Tornielli GB, Fasoli M, Paci P, Tombesi S, Frioni T, Silvestroni O, Bellincontro A, et al. Distinct transcriptome responses to water limitation in isohydric and anisohydric grapevine cultivars. BMC genomics. 2016;17(1):815.
Osipovich AB, Dudek KD, Greenfest-Allen E, Cartailler JP, Manduchi E, Potter Case L, Choi E, Chapman AG, Clayton HW, Gu G, et al. A developmental lineage-based gene co-expression network for mouse pancreatic β-cells reveals a role for Zfp800 in pancreas development. Development (Cambridge, England). 2021;148(6):dev196964.
Köster P, DeFalco TA, Zipfel C. Ca(2+) signals in plant immunity. The EMBO journal. 2022;41(12):e110741.
Xu T, Niu J, Jiang Z. Sensing Mechanisms: Calcium Signaling Mediated Abiotic Stress in Plants. Front Plant Sci. 2022;13:925863.
Frey FP, Pitz M, Schön CC, Hochholdinger F. Transcriptomic diversity in seedling roots of European flint maize in response to cold. BMC genomics. 2020;21(1):300.
Kumar N, Mishra BK, Liu J, Mohan B, Thingujam D, Pajerowska-Mukhtar KM, Mukhtar MS. Network Biology Analyses and Dynamic Modeling of Gene Regulatory Networks under Drought Stress Reveal Major Transcriptional Regulators in Arabidopsis. Int J Mol Sci. 2023;24(8):7349.
Tian J, Ma Y, Chen Y, Chen X, Wei A. Plant Hormone Response to Low-Temperature Stress in Cold-Tolerant and Cold-Sensitive Varieties of Zanthoxylum bungeanum Maxim. Front Plant Sci. 2022;13:847202.
Zhang H, Zhu J, Gong Z, Zhu JK. Abiotic stress responses in plants. Nat Rev Genet. 2022;23(2):104–19.
Lambert SA, Jolma A, Campitelli LF, Das PK, Yin Y, Albu M, Chen X, Taipale J, Hughes TR, Weirauch MT. The Human Transcription Factors. Cell. 2018;172(4):650–65.
Blanc-Mathieu R, Dumas R, Turchi L, Lucas J, Parcy F. Plant-TFClass: a structural classification for plant transcription factors. Trends in plant science. 2024;29(1):40–51.
Bailey TL, Johnson J, Grant CE, Noble WS. The MEME Suite. Nucleic Acids Res. 2015;43(W1):W39–49.
Wang X, Li Z, Shi Y, Liu Z, Zhang X, Gong Z, Yang S. Strigolactones promote plant freezing tolerance by releasing the WRKY41-mediated inhibition of CBF/DREB1 expression. The EMBO journal. 2023;42(19):e112999.
Chen QF, Xu L, Tan WJ, Chen L, Qi H, Xie LJ, Chen MX, Liu BY, Yu LJ, Yao N, et al. Disruption of the Arabidopsis Defense Regulator Genes SAG101, EDS1, and PAD4 Confers Enhanced Freezing Tolerance. Mol Plant. 2015;8(10):1536–49.
Li S, He L, Yang Y, Zhang Y, Han X, Hu Y, Jiang Y. INDUCER OF CBF EXPRESSION 1 promotes cold-enhanced immunity by directly activating salicylic acid signaling. The Plant cell. 2024;36(7):2587–606.
Feng S, Liu Z, Chen H, Li N, Yu T, Zhou R, Nie F, Guo D, Ma X, Song X. PHGD: An integrative and user-friendly database for plant hormone-related genes. iMeta. 2024;3(1):e164.
Benoit M, Drost HG, Catoni M, Gouil Q, Lopez-Gomollon S, Baulcombe D, Paszkowski J. Environmental and epigenetic regulation of Rider retrotransposons in tomato. PLoS Genet. 2019;15(9):e1008370.
Sun L, Jing Y, Liu X, Li Q, Xue Z, Cheng Z, Wang D, He H, Qian W. Heat stress-induced transposon activation correlates with 3D chromatin organization rearrangement in Arabidopsis. Nat Commun. 2020;11(1):1886.
Bendall ML, de Mulder M, Iñiguez LP, Lecanda-Sánchez A, Pérez-Losada M, Ostrowski MA, Jones RB, Mulder LCF, Reyes-Terán G, Crandall KA, et al. Telescope: Characterization of the retrotranscriptome by accurate estimation of transposable element expression. PLoS computational biology. 2019;15(9):e1006453.
Liu Y, Cai Y, Li Y, Zhang X, Shi N, Zhao J, Yang H. Dynamic changes in the transcriptome landscape of Arabidopsis thaliana in response to cold stress. Front Plant Sci. 2022;13:983460.
Dasgupta P, Das A, Datta S, Banerjee I, Tripathy S, Chaudhuri S. Understanding the early cold response mechanism in IR64 indica rice variety through comparative transcriptome analysis. BMC genomics. 2020;21(1):425.
Park S, Lee CM, Doherty CJ, Gilmour SJ, Kim Y, Thomashow MF. Regulation of the Arabidopsis CBF regulon by a complex low-temperature regulatory network. Plant J. 2015;82(2):193–207.
Zeng Z, Zhang S, Li W, Chen B, Li W. Gene-coexpression network analysis identifies specific modules and hub genes related to cold stress in rice. BMC genomics. 2022;23(1):251.
Wang M, Li Z, Zhang Y, Zhang Y, Xie Y, Ye L, Zhuang Y, Lin K, Zhao F, Guo J, et al. An atlas of wheat epigenetic regulatory elements reveals subgenome divergence in the regulation of development and stress responses. The Plant cell. 2021;33(4):865–81.
Graham CA, Paajanen P, Edwards KJ, Dodd AN. Genome-wide circadian gating of a cold temperature response in bread wheat. PLoS Genet. 2023;19(9):e1010947.
Tossi VE, Martínez Tosar LJ, Laino LE, Iannicelli J, Regalado JJ, Escandón AS, Baroli I, Causin HF, Pitta-Álvarez SI. Impact of polyploidy on plant tolerance to abiotic and biotic stresses. Front Plant Sci. 2022;13:869423.
Cheng F, Wu J, Cai X, Liang J, Freeling M, Wang X. Gene retention, fractionation and subgenome differences in polyploid plants. Nat Plants. 2018;4(5):258–68.
Li M, Sun W, Wang F, Wu X, Wang J. Asymmetric epigenetic modification and homoeolog expression bias in the establishment and evolution of allopolyploid Brassica napus. New Phytol. 2021;232(2):898–913.
Xie Y, Ying S, Li Z, Zhang Y, Zhu J, Zhang J, Wang M, Diao H, Wang H, Zhang Y, et al. Transposable element-initiated enhancer-like elements generate the subgenome-biased spike specificity of polyploid wheat. Nat Commun. 2023;14(1):7465.
Ahmad P, Abdel Latef AA, Rasool S, Akram NA, Ashraf M, Gucel S. Role of Proteomics in Crop Stress Tolerance. Front Plant Sci. 2016;7:1336.
Zhao M, Zhang N, Gao T, Jin J, Jing T, Wang J, Wu Y, Wan X, Schwab W, Song C. Sesquiterpene glucosylation mediated by glucosyltransferase UGT91Q2 is involved in the modulation of cold stress tolerance in tea plants. New Phytol. 2020;226(2):362–72.
Kopecká R, Kameniarová M, Černý M, Brzobohatý B, Novák J. Abiotic Stress in Crop Production. Int J Mol Sci. 2023;24(7):6603.
Wang F, Chen X, Dong S, Jiang X, Wang L, Yu J, Zhou Y. Crosstalk of PIF4 and DELLA modulates CBF transcript and hormone homeostasis in cold response in tomato. Plant Biotechnol J. 2020;18(4):1041–55.
Labunskyy VM, Gerashchenko MV, Delaney JR, Kaya A, Kennedy BK, Kaeberlein M, Gladyshev VN. Lifespan extension conferred by endoplasmic reticulum secretory pathway deficiency requires induction of the unfolded protein response. PLoS Genet. 2014;10(1):e1004019.
Yeh S, Moffatt BA, Griffith M, Xiong F, Yang DS, Wiseman SB, Sarhan F, Danyluk J, Xue YQ, Hew CL, et al. Chitinase genes responsive to cold encode antifreeze proteins in winter cereals. Plant Physiol. 2000;124(3):1251–64.
Ali M, Muhammad I, Ul Haq S, Alam M, Khattak AM, Akhtar K, Ullah H, Khan A, Lu G, Gong ZH. The CaChiVI2 Gene of Capsicum annuum L. Confers Resistance Against Heat Stress and Infection of Phytophthora capsici. Front Plant Sci. 2020;11:219.
Zhang T, Zhou L, Pu Y, Tang Y, Liu J, Yang L, Zhou T, Feng L, Wang X. A chromosome-level genome reveals genome evolution and molecular basis of anthraquinone biosynthesis in Rheum palmatum. BMC Plant Biol. 2024;24(1):261.
Li Y, Zhu J, Xu J, Zhang X, Xie Z, Li Z. Effect of cold stress on photosynthetic physiological characteristics and molecular mechanism analysis in cold-resistant cotton (ZM36) seedlings. Front Plant Sci. 2024;15:1396666.
Yu G, Zhang B, Chen Q, Huang Z, Zhang B, Wang K, Han J. Dynamic DNA methylation modifications in the cold stress response of cassava. Genomics. 2024;116(4):110871.
He J, Yao L, Pecoraro L, Liu C, Wang J, Huang L, Gao W. Cold stress regulates accumulation of flavonoids and terpenoids in plants by phytohormone, transcription process, functional enzyme, and epigenetics. Critical reviews in biotechnology. 2023;43(5):680–97.
Su Y, Liang W, Liu Z, Wang Y, Zhao Y, Ijaz B, Hua J. Overexpression of GhDof1 improved salt and cold tolerance and seed oil content in Gossypium hirsutum. Journal of plant physiology. 2017;218:222–34.
Wang S, Guo T, Wang Z, Kang J, Yang Q, Shen Y, Long R. Expression of Three Related to ABI3/VP1 Genes in Medicago truncatula Caused Increased Stress Resistance and Branch Increase in Arabidopsis thaliana. Front Plant Sci. 2020;11:611.
Zhao Y, Liu X, Wang M, Bi Q, Cui Y, Wang L. Transcriptome and physiological analyses provide insights into the leaf epicuticular wax accumulation mechanism in yellowhorn. Hortic Res. 2021;8(1):134.
He J, Tang S, Yang D, Chen Y, Ling L, Zou Y, Zhou M, Xu X. Chemical and Transcriptomic Analysis of Cuticle Lipids under Cold Stress in Thellungiella salsuginea. Int J Mol Sci. 2019;20(18):4519.
Garrigues JM, Tsu BV, Daugherty MD, Pasquinelli AE. Diversification of the Caenorhabditis heat shock response by Helitron transposable elements. Elife. 2019;8:e51139.
Marasca F, Gasparotto E, Polimeni B, Vadalà R, Ranzani V, Bodega B. The Sophisticated Transcriptional Response Governed by Transposable Elements in Human Health and Disease. Int J Mol Sci. 2020;21(9):3201.
Papolu PK, Ramakrishnan M, Mullasseri S, Kalendar R, Wei Q, Zou LH, Ahmad Z, Vinod KK, Yang P, Zhou M. Retrotransposons: How the continuous evolutionary front shapes plant genomes for response to heat stress. Front Plant Sci. 2022;13:1064847.
Butelli E, Licciardello C, Zhang Y, Liu J, Mackay S, Bailey P, Reforgiato-Recupero G, Martin C. Retrotransposons control fruit-specific, cold-dependent accumulation of anthocyanins in blood oranges. The Plant cell. 2012;24(3):1242–55.
Bolger A, Scossa F, Bolger ME, Lanz C, Maumus F, Tohge T, Quesneville H, Alseekh S, Sørensen I, Lichtenstein G, et al. The genome of the stress-tolerant wild tomato species Solanum pennellii. Nat Genet. 2014;46(9):1034–8.
Fiol A, García S, Dujak C, Pacheco I, Infante R, Aranzana MJ. An LTR retrotransposon in the promoter of a PsMYB10.2 gene associated with the regulation of fruit flesh color in Japanese plum. Hortic Res. 2022;9:uhac206.
Orozco-Arias S, Isaza G, Guyot R. Retrotransposons in Plant Genomes: Structure, Identification, and Classification through Bioinformatics and Machine Learning. Int J Mol Sci. 2019;20(15):3837.
Han J, Dai Y, Zhou J, Tian J, Chen Q, Kou X, Raza G, Zhang B, Wang K. Tissue-specific chromatin accessibility and transcriptional regulation in maize cold stress response. Genomics. 2024;117(1):110981.
Han J, Wang P, Wang Q, Lin Q, Chen Z, Yu G, Miao C, Dao Y, Wu R, Schnable JC, et al. Genome-Wide Characterization of DNase I-Hypersensitive Sites and Cold Response Regulatory Landscapes in Grasses. The Plant cell. 2020;32(8):2457–73.
Ambroise V, Legay S, Guerriero G, Hausman JF, Cuypers A, Sergeant K. The Roots of Plant Frost Hardiness and Tolerance. Plant & cell physiology. 2020;61(1):3–20.
Vergara A, Haas JC, Aro T, Stachula P, Street NR, Hurry V. Norway spruce deploys tissue-specific responses during acclimation to cold. Plant, cell & environment. 2022;45(2):427–45.
Sun X, Feng D, Liu M, Qin R, Li Y, Lu Y, Zhang X, Wang Y, Shen S, Ma W, et al. Single-cell transcriptome reveals dominant subgenome expression and transcriptional response to heat stress in Chinese cabbage. Genome Biol. 2022;23(1):262.
Liu Q, Ma W, Chen R, Li T, Wang Q, Wei C, Hong Y, Sun HX, Cheng Q, Zhao J, et al. Multiome in the Same Cell Reveals the Impact of Osmotic Stress on Arabidopsis Root Tip Development at Single-Cell Level. Advanced science (Weinheim, Baden-Wurttemberg, Germany). 2024;11(24):e2308384.
Li P, Liu Q, Wei Y, Xing C, Xu Z, Ding F, Liu Y, Lu Q, Hu N, Wang T, et al. Transcriptional Landscape of Cotton Roots in Response to Salt Stress at Single-cell Resolution. Plant communications. 2023;5(2):100740.
Jarošová J, Prerostova S, Černý M, Dobrev P, Gaudinova A, Knirsch V, Kobzová E, Müller K, Fiala R, Benczúr K, et al. Hormonal responses of rice to organ-targeted cold stress. Environmental and Experimental Botany. 2024;222:105739.
Qiu P, Liu T, Xu Y, Ye C, Zhang R, Wang Y, Jin Q. Multi-omic dissection of the cold resistance traits of white water lily. Hortic Res. 2024;11(6):uhae093.
Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics : a journal of integrative biology. 2012;16(5):284–7.
Zhang T, Hu Y, Jiang W, Fang L, Guan X, Chen J, Zhang J, Saski CA, Scheffler BE, Stelly DM, et al. Sequencing of allotetraploid cotton (Gossypium hirsutum L. acc. TM-1) provides a resource for fiber improvement. Nat Biotechnol. 2015;33(5):531–7.
Han J, Yu G, Zhang X, Dai Y, Zhang H, Zhang B, Wang K. Histone Maps in Gossypium darwinii Reveal Epigenetic Regulation Drives Subgenome Divergence and Cotton Domestication. Int J Mol Sci. 2023;24:10607.
Grant CE, Bailey TL, Noble WS. FIMO: scanning for occurrences of a given motif. Bioinformatics. 2011;27(7):1017–8.
Ou S, Su W, Liao Y, Chougule K, Agda JRA, Hellinga AJ, Lugo CSB, Elliott TA, Ware D, Peterson T, et al. Benchmarking transposable element annotation methods for creation of a streamlined, comprehensive pipeline. Genome Biol. 2019;20(1):275.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Subgroup GPDP. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, Heyne S, Dündar F, Manke T. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 2016;44(W1):W160–5.
Acknowledgements
We are grateful for the funding support.
Funding
This work was supported by the National Natural Science Foundation of China (32070544 and 32200543), the Natural Science Foundation of Jiangsu Province (BK20220604), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX23_3387), and the Zhejiang Provincial Natural Science Foundation of China (LQ20C130005).
Author information
Authors and Affiliations
Contributions
JH and KW planned and designed the research. YD, JZ, and DZ performed the experiments. JH and YD analyzed the data. JH, KW, YD, and BZ wrote and revised the manuscript. All authors read and approved the final version for publication.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
The seeds of TM-1 plants were provided by Prof. Kai Wang at Nantong University. The appropriate permissions were obtained for all materials used in this study. We complied with all relevant institutional, national, and international guidelines and legislation.
Consent for publication
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Dai, Y., Zhou, J., Zhang, B. et al. Time-course transcriptome analysis reveals gene co-expression networks and transposable element responses to cold stress in cotton. BMC Genomics 26, 235 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11433-z
Received:
Accepted:
Published:
DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11433-z