Chromatin Topology Reconfiguration Orchestrates Thermotolerant Male Fertility via GhAL5 in Cotton
Yanlong Li, Jing Yang, Weiran Wang, Chunyang Zuo, Liuling Pei, Yizan Ma, Rui Zhang, Yaru Fan, Huanhuan Ma, Yawei Li, Ruizhen Liu, Shuangxia Jin, Longfu Zhu, Jie Kong, Xianlong Zhang, Ling Min
TL;DR
The study reveals how chromatin structure changes help cotton plants maintain male fertility under heat stress, with a key gene, GhAL5, playing a central role in this adaptation.
Contribution
The study identifies GhAL5 as a chromatin-regulated transcription factor that enhances thermotolerance in cotton and other crops through 3D genome reorganization.
Findings
Heat-tolerant cotton lines exhibit controlled chromatin dynamics, while heat-sensitive lines show structural hyperactivation.
GhAL5 overexpression improves thermotolerance in cotton and rice, indicating cross-species functionality.
Dynamic chromatin topology reorganization is central to male fertility under heat stress.
Abstract
Cotton, a globally vital crop, faces severe yield losses due to heat‐induced male sterility. To decipher the thermotolerance mechanisms, we conducted multi‐omics analyses (3D chromatin architecture, transcriptome, and epigenome profiling) on heat‐tolerant (84021) and heat‐sensitive (H05) lines across critical anther developmental stages. We identified subgenome homoeologous gene expression bias linked to thermotolerance, driven by high temperature (HT)‐induced dynamic chromatin topology reorganization. The sensitive line exhibited aberrant 3D structural hyperactivation during anther dehiscence, causing deleterious gene overexpression. Central to this regulation is GhAL5, an Alfin‐like transcription factor modulated through chromatin loop dynamics and TAD‐like boundary reorganization under heat stress. Functional studies confirmed the pivotal role of GhAL5: overexpression enhanced…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5
FIGURE 6
FIGURE 7- —Development Fund for Xinjiang Talents
- —Project of Fund for Biological Breeding of Stress tolerant and High Yield Cotton Varieties
- —Stable Support to Agricultural Sci‐Tech Renovation
- —Hubei Provincial Natural Science Foundation of China for Distinguished Young Scholars
- —Finance Science and Technology Project of Xinjiang Uyghur Autonomous Region
- —Multi‐Resistant and Machine‐Harvested Cotton Bio‐Breeding Innovation Team
- —China Agriculture Research System
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsPlant Reproductive Biology · Plant Molecular Biology Research · Photosynthetic Processes and Mechanisms
Introduction
1
Cotton stands as an important strategic material and the principal raw material for the textile industry [1, 2]. In the face of escalating global warming, the incidence of high temperature (HT, defined as daytime temperatures above 35°C and nighttime temperatures above 28°C) events has risen, especially the extreme HT episodes in July and August, which are highly synchronized with the critical boll‐setting stage of cotton. These temperature extremes seriously affect the development of male reproductive organs in cotton, leading to diminished pollen viability, anther indehiscence, and other sterility phenotypes. These issues hinder subsequent fertilization processes, culminating in a marked reduction in cotton yield [3, 4]. Therefore, studying the molecular mechanisms underlying cotton anther response to HT stress, mining key genes conferring heat tolerance, and harnessing them to breed HT‐tolerant cotton varieties represent a viable strategy to bolster cotton's resilience to HT and to secure the stability of the cotton industry.
Chromatin is a complex, linear structure within the eukaryotic cell nucleus, composed of DNA, histones, non‐histone proteins, and small amounts of RNA [5, 6]. It exhibits a hierarchical organization, ranging from the fundamental repeating units, nucleosomes, to higher‐order structures [7, 8, 9]. A nucleosome is formed by wrapping of about 147 base pairs of double‐stranded DNA around histones [7]. Nucleosomes further interact to coil and compact DNA, forming chromatin loops. Tthe aggregation of multiple chromatin loops forms stable topologically associated domains (TADs), thereby achieving higher‐order chromatin structure [9]. Studies in animals have revealed that distant cis‐regulatory elements often require chromatin structural proteins, such as CCCTC binding factor (CTCF) and cohesin, to mediate loop formation, bringing them into close spatial proximity with gene promoters to regulate transcription [10, 11]. Mounting evidence suggests that the intranuclear chromatin folding mechanism is crucial for determining the relationship between genome structure and important biological functions, such as DNA damage and repair, gene expression fine‐tuning, nuclear stability maintenance, and DNA replication [12, 13, 14, 15, 16, 17]. Recently, 3D genomics has emerged as a powerful tool for exploring the spatial structure of the genome and its impact on gene expression and regulation [18, 19, 20].
The advent of 3D genomics technologies has furnished a wealth of high‐quality 3D genomic data, significantly advancing our understanding of the chromatin folding state within the cell nucleus [21, 22, 23]. It has been discovered that the plant genome presents a hierarchical structure at varying base resolutions, ranging from low to high: chromatin territories (CTs), chromatin compartments (A/B compartments), TADs, and chromatin loops [24]. In plants, 3D genomics has not only been extensively applied to studies of growth and development [25, 26] but has also been employed to investigate responses to abiotic stresses [27, 28, 29, 30]. For instance, drought‐tolerant cotton lines have demonstrated greater chromatin elasticity compared to drought‐sensitive lines under drought stress [27]. Cold stress in rice has been associated with weakened long‐range cis‐chromatin interactions [28]. HT stress studies in Arabidopsis have revealed a connection between the 3D genome reorganization and the activation of heterochromatic transposons induced by heat stress [29]. In rice, two varieties with distinct thermotolerance levels exhibited changes in the 3D structure following HT treatment [30]. These findings showed that HT stress affects higher‐order chromatin structures at multiple levels, including A/B compartment switching, enlarged TAD size, and the loss of proximal cis‐interactions, which are associated with alterations in chromatin accessibility and gene expression [30]. Despite the widespread application of 3D genomics in abiotic stresses response research across various plant species, there remains a paucity of studies focusing on the 3D genomics of reproductive organs in response to abiotic stress.
Our previous research has established that under HT stress, cotton anthers undergo a complex interplay of epigenetic, transcriptional, and physiological metabolic regulations to counteract HT stress [3, 31, 32, 33, 34]. However, in cotton anthers, it is unclear how the chromatin state coordinates epigenetic modification changes to regulate gene transcription under HT and affects the mechanism of male fertility. Therefore, this study investigated the dynamic alterations in the anther transcriptome, epigenome, and chromatin spatial structure under HT stress in extremely heat‐tolerant/sensitive cotton lines. We depicted the dynamic landscape of the 3D genome in anthers under HT, dissected the molecular mechanism of cotton heat tolerance from the perspective of high‐order chromatin structure cooperatively regulating epigenetic modifications to control transcriptional shifts under HT, identified a pivotal heat response gene ALFIN‐LIKE 5 (GhAL5), modulated by dynamic chromatin changes, and perform functional analysis on this important gene. Notably, we demonstrated that GhAL5 overexpression not only enhanced cotton's heat tolerance but also conferred significant thermotolerance to rice, indicating its potential as a universal heat protector across different plant species. This groundbreaking finding highlights the critical role of chromatin‐regulated genes in plant stress responses and opens up new avenues for developing heat‐tolerant crops through genetic manipulation.
Results
2
Unbiased Expression of Subgenomic Homologous Genes Improves HT Tolerance in Cotton
2.1
To investigate the mechanisms underlying high temperature (HT) tolerance in cotton, we analyzed the transcriptomes of anthers at three HT sensitive stages: tetrad stage (TS), tapetal degradation stage (TDS), anther dehiscence stage (ADS). We compared two cotton lines 84021 (heat‐tolerant) and H05 (heat‐sensitive), which have showed significant differences in HT tolerance in previous studies (Figure 1A) [3]. By counting the differentially expressed genes (DEGs) under HT stress, we observed that the HT‐sensitive line H05 exhibited a higher proportion of upregulated DEGs at both TS and ADS, especially at ADS (Figure 1A; Table S1), suggesting that H05 undergoes more pronounced transcriptional activation under HT stress compared to 84021. To further mine key HT‐responsive factors, we constructed a co‐expression network using all DEGs under HT from both cotton lines, and found that this co‐expression network was divided into nine major modules (Figure 1B; Table S2), with many core genes identified as key stress‐responsive genes, such as CALMODULIN 8 (CAM8) [35] and CLASS A HEAT SHOCK FACTOR 1B (HSFA1b) [36] (Figure 1B).
Expression bias of homologous gene pairs from two subgenomes in cotton anthers under HT. (A) Number of differentially expressed genes (DEGs) in 84021 and H05 under HT. The numbers above the bars represent the proportion of up‐ and down‐regulated genes among all DEGs. (B) Co‐expression network constructed using DEGs under HT, different colors represent different modules, and the size of the dots represents the weight values of the genes. CAM8, CALMODULIN 8; HSFA1b, CLASS A HEAT SHOCK FACTOR 1B; SPL/NZZ, SPOROCYTELESS/NOZZLE. (C) Number of genes expressed in the At and Dt subgenomes in different cotton anther samples. The expressed genes include homologous genes (blue) and non‐homologous genes (purple). (D) Proportion of homologous gene pairs with different expression patterns under HT. (E) Statistics on the biased expression of homologous gene pairs in different samples. (F) Dynamic number of homologous genes with expression bias at TS, TDS, and ADS under HT. At bias refers to the expression level of homologous genes in the At subgenome is higher than that of homologous genes in the Dt subgenome (Fold Change, FC ≥ 2; False Discovery Rate, FDR ≤ 0.05). meaning of Dt bias versus At bias. The dynamic category indicates the dynamic change of expression bias toward At or Dt subgenome. The red line indicates genes with At‐biased expression under HT, the blue line indicates genes with Dt‐biased expression under HT, and the green line indicates genes without biased expression under HT. (G) Statistics on the biased expression of homologous gene pairs in different anther cells. TS: tetrads stage; TDS: tapetal degradation stage; ADS, anther dehiscence stage.; NT, normal temperature; HT, high temperature; EP, epidermis; EN, endothecium; ML, middle layer; T, tapetum; MAM, microspore after meiotic; V, vascular region; C, connective; MC, meiotic cell.
Subgenomic expression bias has been implicated in driving physiological and biochemical changes in polyploid plants, leading to new phenotypes [37, 38]. Cotton is an allopolyploid containing the At and Dt subgenomes, provides an excellent model to study these effects. We further explored the roles of the two subgenomes in responding to HT stress in cotton anthers at the expression level. First, in terms of the total number of expressed genes, the Dt subgenome had a slightly higher number of expressed genes and homologous gene pairs than the At subgenome in cotton anthers (Figure 1C). We then classified the expression patterns of homologous gene pairs under HT into three types: single change (only one gene in the homologous pair showed differential expression in the At or Dt subgenome during the HT response in cotton anthers), common change (both homologs were up‐ or down‐regulated), and opposite change (the expression of the two homologs changed in opposite directions) (Figure 1D). The results showed that the single change pattern was the most common in the all three stages of 84021 and H05 under HT stress, while 84021 had more homologous gene pairs exhibiting the same change pattern than H05 at TS and ADS (Figure 1D). Further analysis revealed that at ADS, 84021 had more homologous gene pairs transitioning from biased expression to unbiased expression compared to H05 under HT (Figure 1E,F), while H05 exhibited a generally higher level of bias (Figure S1), suggesting that unbiased expression of homologous gene pairs may contribute to HT tolerance. To test whether this pattern of subgenomic expression balance also holds true at the single cell level, we revisited our previously published single‐cell transcriptome dataset of tetrad‐stage anthers under HT stress [31]. We analyzed the expression bias of homoeologous gene pairs within each cell cluster. Notably, heat‐insensitive cell types such as epidermal cells (C1) and endothecial cells (C2, C4, C8, C9) maintained highly balanced expression between subgenomes under HT conditions, exhibiting minimal subgenomic bias.In contrast, heat‐sensitive cell types, including microspores (C14) and tapetal cells (C11, C5, C7, C18), showed pronounced subgenomic expression bias under the same conditions (Figure 1G; Figure S2). These results highlight a cell‐type‐specific pattern of subgenomic regulation, suggesting that imbalanced expression between subgenomes may be a molecular signature of heat‐sensitive cell states. Together, these results reveal a negative correlation between subgenomic expression bias and thermotolerance in cotton anthers, both at the tissue and single‐cell levels. Importantly, they suggest that maintaining balanced expression between the A and D subgenomes is may particularly beneficial in cell types that are otherwise vulnerable to heat stress.
Dynamic 3D Genomic Structure of Cotton Anthers in Response to HT
2.2
To explore the chromatin topological structure's role in HT tolerance in cotton anthers, we employed in situ Hi‐C with the same samples used for RNA‐seq to uncover the dynamic alterations in chromatin conformation under HT stress. The high correlation between biological replicates for each developmental stage (Figure S3A) confirmed the reliability of our Hi‐C dataset (Figure S3A). Our analysis yielded at least 600 million unique mapped contact reads for each stage (Table S3), achieving a resolution of 5 Kb (Figure S3B) following established methods [39].
Hi‐C analysis at various resolutions revealed distinct higher‐order chromatin architectures (Figure 2A; Figure S4). At different resolutions, we observed different higher‐order chromatin structures on the Hi‐C heatmaps (Figure 2A; Figure S4). Notably, at the chromosome level, the number of short‐range interactions (≤2 Mb) exceeded long‐range interactions (>2 Mb) (Figure 2B), indicating that short‐range chromatin interactions are more prevalent than long‐range ones in cotton anthers. During the HT response process, intra‐chromosomal interactions increased and inter‐chromosomal interactions decreased in 84021 TS, TDS, ADS and in H05 ADS, while the opposite trend was observed in H05 TS and TDS (Figure 2C). Next, the chromosome contact probability analysis demonstrated a negative correlation with increasing genomic distance, reflecting a gradual decline in interaction strength (Figure 2D). And interestingly, chromatin compaction responses to HT varied between the two lines. H05 showed increased compaction across all stages, while 84021 only experienced compaction at ADS, with decondensation at other stages (Figure 2E; Figure S5). Comparison between the At and Dt sub‐genomes showed that the Dt subgenome was more compact than the At subgenome, and its response to HT seemed to be more dramatic (Figure 2E). This dynamic 3D genome structure of cotton anthers under HT established here provides a valuable resource for understanding the higher‐order chromatin structure and transcriptional regulation during cotton anther development under stress.
High‐resolution three‐dimensional dynamic mapping of chromatin in cotton anthers under HT. (A) Multi‐resolution Hi‐C interaction heat maps of different cotton anther samples. A standardized color scale is provided to facilitate comparison across stages and conditions. The red squares in the diagram represent the max contacts. (B) Total number of Hi‐C contacts for short‐range (≤ 2 Mb) and long‐range (> 2 Mb) chromatin interactions under NT and HT in different samples. Each data point represents one of the 26 chromosomes of cotton (n = 26). Statistical significance was assessed using a two‐sided Wilcoxon signed‐rank test. (C) Proportion of the number of intra‐ and interchromosomal contacts in each sample of anther. (D) Chromatin contact probabilities calculated relative to genomic distances from 1 Kb to 100 Mb in each sample of the anther. (E) Comparative analysis of chromatin compactness under HT. Each boxplot represents the distribution of compactness values for the A and D subgenomes across different conditions. Sample sizes: A subgenome (n = 141,397), D subgenome (n = 82,014). Different lowercase letters above bars indicate statistically significant differences (p < 0.001) based on two‐sided Wilcoxon signed‐rank test was used. TS, tetrads stage; TDS, tapetal degradation stage; ADS, anther dehiscence stage; NT, normal temperature; NT, normal temperature; HT, high temperature.
Compartment Switching Reveals Functional Divergence Between At and Dt Subgenomes, with At Subgenome May Governing Heat Stress Adaptation
2.3
Chromatin compartment, which is the next level below chromatin domains in the chromatin hierarchy. At this level, we categorized the chromatin of each sample into A and B compartments using a 40 Kb Hi‐C matrix, representing active and inactive chromatin regions, respectively. A/B compartments were pervasive throughout the genome, encompassing approximately 97% of the genomic regions (Figure S6A). On this basis, we found that the genomic length occupied by the B compartment exceeded that of the A compartment in all stages by examining the proportions of genomic length occupied by the A and B compartments (Figure 3A). But the A compartment had a higher gene density and expression level compared to the B compartment in all samples (Figure S6B; Figure 3B). Furthermore, the genomic proportions of A/B compartments altered under HT, particularly at the ADS stage, where significant differences between 84021 and H05 were observed. In 84021, the proportion of the A compartment on the genome decreased under HT, while in H05, it increased obviously (Figure 3A). These results indicate dynamic changes in A/B compartments states under HT. So, we investigated the dynamic conversions of A/B compartment states in all samples (Figure 3C,D). As an illustrative example, we first visualized the chromatin compartment structures in the 11 Mb‐25 Mb region of chromosome A03 at different anther stages under HT stress, as this region showed clear dynamic changes in compartment states (Figure 3C). Further observation of the changes in A/B compartments under HT revealed that most A/B compartments remained in a relatively stable state, yet some regions underwent dynamic conversion. At TS, more regions converted from the B to the A compartment (B to A) in both cotton lines, with a higher proportion in H05. At TDS, more regions converted from the B to the A compartment in 84021 compared to H05, while at ADS, a large number of regions converted from the B to the A compartment in H05, but the opposite trend occurred in 84021, with many regions converting from the A to the B compartment (AtoB) (Figure 3D). Similarly, the number of genes contained in the dynamically converted A/B compartment regions was consistent with the size of the converted regions (Figure 3E,F). Such compartment conversions occurred on every chromosome (Figure S7), suggesting that at the chromatin A/B compartment level, H05 was in a more active state during the TS and ADS stages under HT, while 84021 was in a relatively silenced state, which is consistent with the higher proportion of upregulated genes in H05 under HT stress (Figure 1A).
*Dynamic changes in the A/B compartment in cotton anthers under HT. (A) The bar plot depicted the proportion of genomic length categorized as A/B compartment in different anther samples. (B) Gene expression levels in A compartment are significantly higher than those in B compartment. The sample size (n) for each box plot, corresponding to the number of genes in each compartment/sample group, is provided in Table S8. Two‐sided Wilcoxon signed‐rank test was used (**p < 0.001). FPKM, fragments per kilobase of transcript per million mapped reads. (C) A representative genomic region showing the dynamic state switch between A/B compartments in chromosome A03. Red color represents region A compartment and blue color represents region B compartment. Some regions showing A/B compartment switching are marked in the red box. (D) Global dynamic switching of chromatin compartment states under HT. Chromosomal regions with stable chromatin states are shown in red (A compartment) or light blue (B compartment). Heatmap showing chromatin compartments whose chromatin compartment state switches between A compartment and B compartment. “AB” indicates switching from A to B compartment and “BA” indicates switching from B to A compartment. (E) Size of the A/B compartment switching region in each sample under HT. (F) The number of genes contained in the A/B compartment switching region in each sample under HT. (G) The bar plot shows the proportion of differentially expressed genes contained in switching regions. (H) The size of the A/B compartment switching regions at the subgenome level in different samples under HT. (I) Proportion of active and inactive genomic regions contained in the A/B compartments, respectively, in different anther samples at the subgenome level. TS, tetrads stage; TDS, tapetal degradation stage; ADS, anther dehiscence stage; NT, normal temperature; HT, high temperature.
We further wondered whether the dynamic conversion of A/B compartments during the HT response of anthers would affect gene expression. The expression levels of genes in the AtoB regions and BtoA regions were analyzed, the result showed genes located in the AtoB regions were more likely to be downregulated, while genes in the BtoA regions tended to be upregulated under HT (Figure 3G), this is consistent with studies suggesting that large‐scale 3D structural changes are associated with gene expression. Comparing the A/B compartment conversion between the At and Dt subgenomes, we found that at the TS stage, the At subgenome had more BtoA conversions than the Dt subgenome in both H05 and 84021. At the TDS stage, in both lines, the At subgenome had more regions exhibited BtoA, while the Dt subgenome had more regions exhibited AtoB. At the ADS stage, the dynamic conversion patterns of A/B compartments under HT differed obviously between the two lines. In 84021, the primary conversion pattern in both subgenomes was AtoB, while in H05, the main conversion pattern, primarily occurring in the At subgenome, was BtoA (Figure 3H). These results suggest that the At subgenome may be more responsive to HT stress, with the Dt subgenome exhibiting a more stable chromatin organization under HT.
We further performed functional enrichment analysis on the genes contained in the A/B compartment conversion regions of the At and Dt subgenomes separately (Figure S8A) revealed that genes in the At subgenome conversion regions were enriched in processes such as gene expression regulation and phosphorylation, while the genes in the Dt subgenome conversion regions were enriched in processes such as vesicle fusion, hydrogen ion transmembrane transport, and acetyl‐CoA metabolism (Figure S8B). These results demonstrated that the At and Dt subgenomes may contribute differentially to anther response to HT at the compartment level.
Further, to understand how the compartment structure affects gene expression under HT, we integrated our compartment data with ChIP‐seq and RNA‐seq data. Based on the presence and type of histone modifications (H3K27me3 and H3K4me3), we categorized the compartment structure regions into four classes (unmarked; H3K4me3‐marked regions; H3K27me3‐marked regions; regions marked by both H3K4me3 and H3K27me3) (Figure S9). We found that genes in the H3K4me3‐marked regions, typically associated with active regions, had higher expression levels compared to those in the H3K27me3‐marked regions and unmarked regions (Figure S9). Furthermore, active regions (with H3K4me3 modification) tended to be located in the A compartment, while inactive regions (unmarked or only marked by H3K27me3) predominantly located in the B compartment (Figure 3I; Figure S9). Interestingly, at the TDS stage, the proportions of active and inactive regions in the A and B compartments remained relatively stable under HT in both 84021 and H05. However, at ADS under HT, the proportion of active regions in the A compartment showed a decreasing trend in 84021, while the proportion of inactive regions in the B compartment increased. In contrast, H05 exhibited the opposite trend, the proportion of active regions in the A compartment increased, and the proportion of inactive regions in the B compartment decreased (Figure 3I). This trend was more pronounced in the At subgenome (Figure 3I). Moreover, the changes in active regions under HT were consistent with the expression changes in the two cotton lines (Figure 1A). These results suggest a potential tendency of the At subgenome to respond more dynamically to HT stress, possibly through coordinated chromatin reorganization and epigenetic regulation. However, given that Hi‐C data were derived from merged replicates, further validation with additional biological replicates is required to confirm the subgenome‐specific responses.
TE‐driven Rewiring of TAD‐Like Structures Boundaries Balances Insulation Stability and Transcriptional Agility Under HT in Cotton Anthers
2.4
In the 3D genome, the next level of chromatin compartment structure is the TAD‐like structure. TAD‐like structures are critical for organizing regulatory elements and target genes within the 3D chromatin space, and their reorganization can alter chromatin interactions near TAD‐like structures boundaries, potentially affecting gene expression [40, 41]. To investigate the impact of HT on the 3D organization of TAD‐like structures within the cotton anther genome, we employed a 20 Kb resolution matrix to identify and analyze TAD‐like structures in different samples. We found that the number of TAD‐like structures increased slightly under HT in different samples, while the TAD‐like structures size did not change obviously, except for a decrease in 84021 at TS under HT (Figure 4A). After classifying the TAD‐like structure boundaries into dynamic (those that gained or lost interactions under HT) and conservative (those that maintained interactions), we found that the majority (at least 76% in each sample) of TAD‐like structure boundaries were conservative under HT (Figure 4B). This suggests that most TAD‐like structures maintained their integrity under HT (Figure S10). The proportion of dynamic boundaries was slightly higher in 84021 than in H05 at TS and TDS stages, but lower at ADS under HT. We used insulation scores to evaluate boundary strength and found that conservative boundaries exhibited stronger insulation than dynamic ones (Figure 4C), consistent with previous findings [25]. Representative examples of both boundary types are provided (Figure 4D; Figure S11). These results indicate that, although some TAD‐like boundaries are responsive to HT, the overall TAD‐like structures landscape remains largely stable, underscoring the role of conserved TAD‐like architecture in maintaining chromatin organization during stress responses.
*Global dynamic TAD‐like structure organization under HT. (A) Size and number of TAD‐like structures in each cotton sample. (B) Proportion of dynamic TAD‐like structure boundaries and conserved TAD‐like structure boundaries in each sample under HT. (C) Insulation scores at HT dynamic TAD‐like structure boundaries and at HT conserved TAD‐like structure boundaries. (D) Demonstration of a representative dynamic TAD‐like structure under HT. The insulating scores are represented below. Orange lines represent conserved TAD boundaries, while blue lines represent dynamic TAD boundaries. (E) Expression of genes located in TAD‐like structure boundary or TAD‐like structure interior. (F) Number of differentially expressed genes located at the HT dynamic TAD‐like structure boundary. (G) Expression of TEs located in TAD‐like structure boundary or TAD‐like structure interior. (H) Expression of TEs located at the HT gained, HT loss, and HT conserved TAD‐like structure boundary. In (E, G and H), two‐sided Wilcoxon signed‐rank test was used (***P < 0.001; **P < 0.01; P < 0.05; ns: no significance). Sample size (n) information in (E, G and H) is provided in Table S8. TS, tetrads stage; TDS, tapetal degradation stage; ADS, anther dehiscence stage; TAD, topologically associated domains; NT, normal temperature; HT, high temperature; TE, transposable element; FPKM, fragments per kilobase of transcript per million mapped reads.
Next, we investigated the impact of TAD‐like structures on gene transcription. Comparison of the coverage and expression levels of genes at TAD‐like structure boundaries with those at TAD‐like structures interiors in cotton anther samples, showed that gene coverage was higher at TAD‐like structure boundaries than within TAD‐like structures in anther (Figure S12), indicating that more genes were preferentially located at TAD‐like structure boundaries. Additionally, the expression levels of genes at TAD‐like structure boundaries were consistently higher than those of within TAD‐like structures (Figure 4E), indicating that these boundaries are critical for gene regulation. We further compared the expression changes of genes at dynamic TAD‐like structure boundaries under HT and found that, compared to H05, these boundary changes in 84021 were associated with the expression changes of more genes at each stage (Table S4). However, among the DEGs at dynamic TAD‐like structure boundaries under HT, the proportion of upregulated genes was higher in H05 than in 84021 (Figure 4F; Table S4), demonstrating that the changes in TAD‐like structures under HT response in H05 is associated with more widespread activation of genes at these boundaries.
The activation of transposable elements (TEs) has been reported to be related to the formation of TAD‐like boundaries [42], we sought to investigate whether the dynamics of TAD‐like boundaries in cotton anthers under HT are related to TE activity. In cotton anthers, TEs are primarily enriched within TAD‐like structures (Figure S12), and their expression of TEs is significantly higher at TAD‐like structure boundaries than at TAD‐like structures interiors (Figure 4G). Notably, the expression of TEs at newly formed TAD‐like structure boundaries is higher than that at lost boundaries under HT, but lower than at conserved TAD‐like structure boundaries (Figure 4H). These findings suggest that the establishment and maintenance of TAD‐like boundaries under HT are associated with heightened TE transcriptional activity. These analyses record the details of the reorganization of TAD‐like structures in cotton anthers in response to HT and provide clues for studying the transcriptional regulation of genes around these boundaries.
Spatial Genome Reorganization Under HT Alters Subgenome Expression Balance via TE‐Mediated TAD‐Like Structures Dynamics in Cotton Anthers
2.5
The impact of TAD‐like structure alterations on subgenome‐specific gene expression under HT was a key focus of our investigation. We observed that the conservation and specificity of TAD‐like structures in the At and Dt subgenomes under HT were similar (Figure S13A). At the same time, we noted that the size and number of TAD‐like structures in the At subgenome are generally larger than those in the Dt subgenome in each anther sample (Figure S13B), which may be related to the larger size of the At subgenome. Based on the consistency of homologous genes contained in the TAD‐like structures of the two subgenomes, we defined the TAD‐like structures as homologous when they contained consistent homologous genes across both subgenomes and as partitioned when homologous genes were distributed across different TAD‐like structures between the subgenomes (Figure 5A). A total of 7 968 homologous TAD‐like structures were identified in all samples (Table S5). The comparison between the At and Dt subgenomes showed that the number and length of partitioned TAD‐like structures in the At subgenome are more than those in the Dt subgenome (Figure 5B), this is also likely due to the larger size of the A subgenome. Furthermore, we found that the expression of TEs at the boundaries of partitioned TAD‐like structures was significantly higher than at the boundaries of homologous TAD‐like structures (Figure 5C; Figure S14). This indicates that the transcriptional activity of TEs is related to the transformation of subgenome‐specific TAD‐like structures boundaries, leading to the partitioning or merging of TAD‐like structures between the At and Dt subgenomes. This is consistent with the pattern in cotton fibers [25].
*Different organization of TAD‐like structures at the subgenome level leads to homologous gene expression bias. (A) Examples of homologous TAD‐like structures and partitioned TAD‐like structures are shown. Homologous gene pairs are connected with gray lines. Non‐homologous genes are not connected. (B) Number of partitioned TAD‐like structures on At subgenome and Dt subgenome in different samples. (C) TE expression at the boundary of homologous TAD‐like structures and at the boundary of partitioned TAD‐like structures. Sample size (n) information is provided in Supplementary Table S8. Two‐sided Wilcoxon signed‐rank test, ***P <0.001. (D) Proportion of homologous TAD‐like structures and partitioned TAD‐like structures on the genome in each sample. (E) Proportion of HT dynamic TAD‐like structure boundaries (HT‐gain and HT‐loss) occurring at homologous TAD‐like structure boundaries and partitioned TAD‐like structure boundaries. (F) The number of homologous gene pairs within homologous TAD‐like structures and partitioned TAD‐like structures in each sample of anthers. (G)The number of homologous genes at different spatial locations (located at the border or interior) of TAD‐like structures of the two subgenomes and showing expression bias. Comparison with a random distribution of homologous genes showing expression bias. Pearson chi‐square test, **P <0.001. (H) Differences in biased expression of GhUBC32 homologous genes in 84021 and H05 are associated with changes in TAD‐like structures at the subgenome level under HT. TS, tetrads stage; TDS, tapetal degradation stage; ADS, anther dehiscence stage; NT, normal temperature; HT, high temperature; TAD, topologically associated domains; TE, transposable element; UBC32, UBIQUITIN‐CONJUGATING ENZYME 32; FPKM, fragments per kilobase of transcript per million mapped reads.
Although partitioned TAD‐like structures account for the vast majority of the genome (Table S6), the proportion of homologous TAD‐like structures and partitioned TAD‐like structures both exhibited dynamic changes under HT (Figure 5D). Under HT, at TS, the proportion of homologous TAD‐like structures increased slightly in both 84021 and H05. At TDS, the proportions of the two types of TAD‐like structures remained basically unchanged in 84021 and H05. However, at ADS, the proportion of homologous TAD‐like structures increased slightly in 84021, but decreased obviously in H05 (Figure 5D). Moreover, at TS, the proportion of newly gained TAD‐like structures boundaries compared to lost TAD‐like structures boundaries under HT was higher at the boundaries of homologous TAD‐like structures in both 84021 and H05. At TDS, in 84021, the proportion of newly acquired TAD‐like structure boundaries compared to lost TAD‐like structures boundaries under HT was similarly higher at the boundaries of homologous TAD‐like structures, while in H05, the proportions were not significantly different. At ADS, in 84021, the proportion of newly acquired TAD‐like structures boundaries compared to lost TAD‐like structures boundaries under HT was similarly higher at the boundaries of homologous TAD‐like structures, but in H05, the opposite was observed (Figure 5E).
Next, we explored the effects of different TAD‐like structures on gene expression under HT. We found that the expression of genes at the boundaries of partitioned TAD‐like structures was significantly higher than that of genes at the boundaries of homologous TAD‐like structures, and this trend persisted under HT stress (Figure S14). This suggests that partitioned TAD‐like structures are associated with enhanced transcriptional activity. Then, when comparing the number of homologous genes contained within homologous and partitioned TAD‐like structures, we found the number of homologous genes contained within partitioned TAD‐like structures was far more than that within homologous TAD‐like structures (Figure 5F). This observation indicates that partitioned TAD‐like structures are enriched for homologous genes, which could contribute to their higher transcriptional activity. Furthermore, our analysis revealed that homologous genes with biased expression preferentially localized to regions where the spatial positions of TAD‐like structures changed (Figure 5G). Therefore, the expression bias of homologous genes under HT may be related to changes in the spatial reorganization of TAD‐like structures. Combining the above results, we observed that at ADS under HT, 84021 would gain more homologous TAD‐like structures compared to H05, and a higher proportion of the newly acquired TAD‐like structure boundaries in 84021 occurred at the boundaries of homologous TAD‐like structures. This may explain why 84021 has a greater proportion of co‐expressed homologous genes than H05 under HT. For instance, UBIQUITIN‐CONJUGATING ENZYME 32 (UBC32), known to be a negative regulator of drought stress in Arabidopsis [43], demonstrates a differential expression pattern in 84021 and H05 under HT. In 84021, the homologous gene of UBC32 in the Dt subgenome would transition from the boundary to the interior of the TAD‐like structure, causing the homologous gene pair of UBC32 to change from biased expression toward the Dt subgenome to unbiased expression between the At and Dt subgenomes. In contrast, in H05, under both NT and HT, the UBC32 gene in the A subgenome is always located within the TAD‐like structure interior region, while the UBC32 gene in the D subgenome remains at the TAD‐like structure boundary, so the expression of this gene is always bias toward the D subgenome (Figure 5H). This differential expression pattern may contribute to the differential tolerance of HT between 84021 and H05. In summary, these data indicate that the transcriptional activity of TEs is intricately linked to the reorganization of TAD‐like structures under HT, which in turn influences the spatial positioning of homologous genes in these TAD‐like structures. These structural changes further alter the expression bias of homologous gene pairs under HT.
Chromatin Loop Dynamics Imbalance Underlies Heat‐Induced Male Sterility: Subgenome‐Coordinated G‐ON Loop Hyperactivation Drives Transcriptional Dysregulation in Cotton Anthers
2.6
Chromatin loops are important units for gene transcription regulation, referring to the chromatin structure in which two linearly distant genomic elements (e.g., enhancers, silencers, promoters, etc.) become spatially close through the mediation of proteins and RNAs [10, 11]. Utilizing a 5 Kb resolution matrix, we identified chromatin loops at different stages under HT stress (Figure 6A). Except for the TDS in 84021, the number of chromatin loops increased under HT in other samples (Figure 6A). This suggests that chromatin loops exhibit a higher degree of dynamism compared to TAD‐like structures under HT (Figure 6B; Figure S15). For instance, at TS stage of H05, only 30.3% of the loops were conserved under HT (Figure 6C), highlighting the high dynamism of chromatin loops during HT. Moreover, we found that the number and proportion of dynamic loops between anther samples at different developmental stages under the same temperature condition were obviously lower than those observed between samples under HT and NT conditions. This indicates that the observed loop dynamics are primarily driven by heat stress rather than developmental variation (Figure S15). Additionally, under HT, in 84021 and H05, we observed that the proportion of chromatin loop anchors located in the A compartment increased at TS, decreased at TDS, and further decreased at ADS, and with a more pronounced decrease in 84021 compared to H05 (Figure 6D). Moreover, the changes in the distribution of chromatin loop anchors in the A/B compartments under HT are completely consistent with the changes in their distribution in gene‐rich and gene‐poor regions (Figure 6D). Additionally, at the chromosome level, the changes in the distribution of chromatin loop anchors along the chromosomes were not very obvious under HT at TS and TDS. However, at ADS, the density of loops at the chromosome ends decreased greatly, while the density of loops around the chromosome centers increased in 84021 under HT, but the changes were less significant in H05 (Figure 6E). These results indicate that at ADS under HT, many chromatin loops and interactions in the active genomic regions were lost in 84021, while the changes were relatively minor in H05, and H05 was in a more active state at the chromatin loop level under HT, which is basically consistent with the trends of changes in the A/B compartments under HT (Figure 2E).
*The HT‐sensitive cotton line has more active chromatin loops under HT. (A) Number of chromatin loops in each sample. (B) The heatmap shows the dynamic changes of chromatin loops under HT. Black circles represent chromatin loops. (C) Chromatin loops exhibit high dynamics under high temperature. This sample is from the TS of H05. (D) The bar plots depict the normalized proportions of chromatin loops in the A/B compartments (upper panel), and in the gene‐poor and gene‐rich regions (lower panel) for each sample. (E) The heatmap and line plot show the changes in loop density on chromosome D02 during the HT response. (F) Proportions of G‐G, G‐N, N‐N loop types in each sample. (G) Number of genes with dynamic changes of G‐G loops between 84021 and H05 under HT. (H) Chromatin accessibility levels of G‐G loops and G‐N loops show that the chromatin accessibility of highly expressed genes and their linked genes or intergenic regions is higher compared to non‐expressed genes in anthers. Highly expressed genes: FPKM (fragments per kilobase of transcript per million mapped reads)>100; non‐expressed genes: FPKM< 1. The target genes and their loop regions are located to the left and right of the loop regions in the figure, respectively. Sample size (n) information is provided in Supplementary Table S8. This sample is from the TS of H05 under normal temperature. (I) Gene expression levels of G‐G loops show that highly expressed genes and their linked genes have higher expression levels compared to non‐expressed genes in anthers. Sample size (n) information is provided in Table S8. This sample is from the TDS period of H05 under normal temperature. (J) In G‐N loops, genes connected to open chromatin regions have higher expression than those connected to non‐open chromatin regions. (K) Proportions of G‐ON and G‐CN in each sample. G‐ON loops, where the intergenic region is located within an open chromatin region, and G‐CN loops, where the intergenic region resides in a closed (non‐accessible) chromatin region. (L) Loops number, expression levels, and epigenetic states around the GhCLASP gene at the ADS stage in 84021 and H05. The red lines represent loops related to GhCLASP. (M) The larger the difference in the proportions of G‐G or G‐N type chromatin loops between homoeologous gene pairs on the subgenomes, the more likely it is for homoeologous gene expression bias to occur. This sample is from the ADS stage of H05 under high temperature. For both Figure H and I, two‐sided Wilcoxon signed‐rank test was used (**p < 0.001). TS, tetrads stage; TDS, tapetal degradation stage; ADS, anther dehiscence stage; NT, normal temperature; HT, high temperature. CLASP, CLIP‐ASSOCIATED PROTEIN. G‐G, gene‐gene loop; G‐N, gene‐non‐coding region loop.
To further investigate the relationship between changes in chromatin loops and gene transcriptional regulation under HT, we observed the correlation between the number of chromatin loop anchors on genes and their expression levels. The results showed that in cotton anthers, the number of chromatin loop anchors was initially positively correlated with gene expression levels, but as the number of chromatin loop anchors increased, the correlation became negative. This phenomenon was consistent across all cotton anther samples and was not affected by HT stress (Figure S16). To explain this phenomenon, we categorized chromatin loops into three types: G‐G (loops between genes), G‐N (loops between genes and non‐genic regions), and N‐N (loops between non‐genic regions) (Figure 6F). We found that the number of G‐G loops were positively correlated with gene expression, while the number of G‐N loops were negatively correlated (Figure S17). In all cotton anther samples, as the number of chromatin loops on genes increased, the proportion of G‐G loops increased initially but decreased as the number of chromatin loops continued to increase, while the proportion of G‐N loops increased (Figure S18). This change in the ratio of G‐G to G‐N loops on genes accounts for the initial positive correlation and subsequent negative correlation between the number of chromatin loops and gene expression. Across the genome, G‐G type loops predominate, with G‐N and N‐N type loops being relatively rare. Under HT, at TS, the proportions of G‐G and G‐N type chromatin loops increased in both 84021 and H05. At TDS, the proportions of the three types of chromatin loops remained relatively stable in both samples. At ADS, the proportion of G‐G type chromatin loops decreased in both 84021 and H05, but the decrease was more drastic in 84021 (Figure 6F). This result also indicates that at ADS, a higher proportion of active chromatin loops were lost in 84021 compared to H05. Next, we statistically analyzed the dynamic changes of G‐G loops under HT in each sample and identified that the number of genes with highly dynamic G‐G loops was significant higher in H05 than in 84021 (Figure 6G), indicating that H05 exhibited more pronounced chromatin reorganization in response to HT stress compared to 84021.
Trans‐regulatory elements and trans‐acting factors regulate target gene expression through long‐range chromatin interactions in the intergenic regions [44]. In cotton anthers, the abundance of G‐N type chromatin loops, which connect genes with non‐genic regions, and these G‐N type chromatin loops change dynamically under HT. To better elucidate the regulatory role of G‐N type chromatin loops under HT, we performed chromatin accessibility sequencing (ATAC‐seq) on the same samples used for in situ Hi‐C. The ATAC‐seq data had two replicates for each sample, with high Pearson correlation coefficients (Figure S19), confirming the high reproducibility of our ATAC‐seq results. In cotton anthers, genes with increased chromatin accessibility exhibited significantly higher expression levels compared to those with reduced chromatin accessibility, indicating that a positive correlation between chromatin accessibility and gene expression (Figure S20), in agreement with previous studies [31, 45]. We also found that genes or intergenic regions connected to highly expressed genes through chromatin loops had higher chromatin accessibility than those connected to non‐expressed genes (Figure 6H; Figure S21). Moreover, the expression of genes linked to highly expressed genes through chromatin loops was higher than that of genes connected to non‐expressed genes (Figure 6I; Figure S22). Therefore, under HT, chromatin loops may regulate gene expression by affecting chromatin accessibility. Furthermore, we further classified G‐N type chromatin loops‐those connecting a gene to an intergenic region‐into two categories: G‐ON loops, where the intergenic region is located within an open chromatin region, and G‐CN loops, where the intergenic region resides in a closed (non‐accessible) chromatin region. The expression of genes connected by G‐ON loops was higher than that of genes connected by G‐CN loops (Figure 6J; Figure S23). By calculating the proportions of G‐ON and G‐CN in the G‐N chromatin loops during different stages of 84021 and H05 under HT, we found that the proportion of G‐ON was higher in H05 than that in 84021 at three stages under HT (Figure 6K), indicating that H05 possesses a more active chromatin loops landscape under HT. Our findings indicate that H05 is in a state of heightened chromatin loop activity compared to 84021 under HT, which may lead to its male sterility through excessive chromatin activation. For example, CLIP‐ASSOCIATED PROTEIN (CLASP) gene encodes a microtubule‐associated protein involved in cell division and proliferation [46], exhibited a slight increase in the number of G‐G and G‐N type chromatin loops and gene expression in 84021 under HT at ADS, while H05 showed a significant increase in both the number of G‐G and G‐ON type chromatin loops and CLASP gene expression (Figure 6L). Similarly, at TS, the number of the chromatin remodeling ATPase CHROMATIN REMODELING 1 (CHR1) [47] displayed a decrease in G‐G type chromatin loops in 84021, while H05 exhibited a significant increase. And the number of G‐N type loops of CHR1 increased in 84021 but remained stable in H05, resulting in downregulation of CHR1 expression in 84021 and upregulation in H05 under HT (Figure S24). To investigate whether chromatin loops in cotton anthers regulate the expression of homologous genes between different subgenomes, and to assess the differential effects of G‐G and G‐N type loops on gene expression, we divided the homologous gene pairs into four types: 1) A subgenome homolog with 75%–100% G‐G loops and D subgenome homolog with 0%–25%; 2) A subgenome homolog with 50%–75% G‐G loops and D subgenome homolog with 25%–50%; 3) A subgenome homolog with 25%–50% G‐G loops and D subgenome homolog with 50%–75%; 4) A subgenome homolog with 0%–25% G‐G loops and D subgenome homolog with 75%–100%, based on the proportions of G‐G and G‐N type loops on the homologous genes from different subgenomes (Figure 6M; Figure S25). We found that a greater disparity in the distribution of G‐G or G‐N type chromatin loops between homologous genes corresponded to a higher likelihood of expression bias (Figure 6M; Figure S25). In summary, our analysis suggests that one of the reasons why H05 is prone to male sterility under HT may be the excessive activation of the chromatin state under HT, leading to an inability to effectively cope with HT stress.
Chromatin‐Regulated GhAL5 Mediated 3D Chromatin Rebalancing Drives Cross‐Species Thermotolerance by Coordinating Male Fertility Under Heat Stress
2.7
To dissect the functional significance of chromatin dynamics in HT stress adaptation, we focused on ALFIN‐LIKE 5 (AL5), a chromatin‐regulated transcription factor. GhAL5 was prioritized as a candidate through integration of 3D structural genomic dynamics (Figures 4F and 6G) and co‐expression network analysis of HT‐responsive DEGs (Figure 1B). Specifically, GhAL5 emerged as a hub gene within a key heat‐responsive co‐expression module and was concurrently located within dynamic TAD boundaries and/or associated with differential chromatin loops under heat stress. It was among the top candidates consistently supported by both datasets Figure S26A). Strikingly, GhAL5 exhibited contrasting chromatin positioning between the heat‐tolerant (84021) and heat‐sensitive (H05) lines. In line 84021, GhAL5 consistently localized at the boundaries of TAD‐like structures within ADS under both normal temperature (NT) and high‐temperature (HT) conditions, with abundant chromatin loops associated with GhAL5 observed across both regimes. This spatial configuration correlated with its constitutively high expression (Figure 7A). Conversely, in line H05, GhAL5 initially resided within TAD‐like structure interiors under NT but migrated to TAD‐like structures boundaries under HT stress. This relocation coincided with a marked increase in chromatin loop formation and significantly elevated GhAL5 expression levels under heat treatment (Figure 7A). This chromatin repositioning suggests GhAL5 accessibility is dynamically regulated in heat‐sensitive genotypes under stress. To validate whether the dynamic chromatin positioning of GhAL5 observed in Hi‐C data reflects locus‐specific chromatin interactions, we performed Capture‐C experiments targeting the GhAL5 locus using anther tissues from both 84021 and heat‐sensitive H05 lines under NT and HT conditions at the ADS stage. Strikingly, under heat stress, the GhAL5 locus exhibited a marked increase in chromatin loop interactions, particularly in the heat‐sensitive H05 line (Figure S27). This observation is consistent with the genome‐wide loop gain observed under HT and further supports that GhAL5 is embedded within a locally responsive chromatin interaction hub (Figure 7A). These results provide direct evidence that GhAL5 is involved in heat‐induced chromatin reorganization through dynamic modulation of chromatin contacts.
*GhAL5 is a potential thermotolerant gene. (A) Chromatin status, expression levels and epigenetic states near GhAL5 at the ADS under NT and HT in 84021 and H05. (B) Sequence analysis of GhAL5 mutants Ghal5‐ko at the target sites. The PAM motif is marked in green. The black dotted line represents deletion compared with WT. (C) Relative expression of GhAL5 in cotton anthers of transgenic lines. WT, wild type; Ri2 and Ri7, RNA interference lines; OE1 and OE7, Overexpression lines. The expression of GhUB7 (Gossypium hirsutum UBIQUITIN 7) was used as internal control. (D) The existence of MYC‐GhAL5 protein in GhAL5 overexpressing plants was verified by western blotting. (E) Phenotypes of different genetic materials of GhAL5 at seedling stage under HT. The pictures on the left are leaf phenotypes of different genetic materials of GhAL5 under HT. (F) Examination of male fertility of different genetic materials of GhAL5 under NT and HT. Scale bars above = 1 cm. Scale bars below = 200 µm. The red pollen grains represent fertile pollen grains and the white pollen grains represent sterile pollen grains. (G‐H) Positive detection (G) and expression level analysis (H) of the GhAL5 gene in GhAL5 overexpressing rice lines. (I) HT tolerance assay of GhAL5 overexpressing rice. Scale bars below = 4 cm. (J) Comparative analysis of A/B compartment transitions between wild‐type and various GhAL5 genetic lines revealed distinct patterns. Statistical analyses were performed using a two tailed unpaired Student’ s t test in (C and H), *P <0.01. Values are means ± SD from 3 independent biological replicates. NT, normal temperature; HT, high temperature; ADS, anther dehiscence stage; PAM, protospacer adjacent motif.
Moreover, given that GhAL5 is significantly upregulated by HT, particularly at the TDS and ADS in both H05 and 84021, with especially strong induction at the ADS stage in the heat‐sensitive line H05 (Figure 7A; Figure S26B), we hypothesized that this gene may act as a positive regulator controlling the basal thermotolerance in cotton. To further investigate the function of GhAL5 in the HT response, we generated overexpression, RNAi, and CRISPR knockout lines for GhAL5 in cotton. Quantitative analyses revealed elevated GhAL5 transcript levels and robust MYC‐tagged protein accumulation in overexpression lines (OE1/OE7) (Figure 7C,D). Under 48°C seedling‐stage heat stress, GhAL5‐OE plants exhibited significantly enhanced thermotolerance compared to wild‐type (Figure 7E). Crucially, GhAL5‐OE cotton maintained pollen viability at 78%–83% under heat treatment, starkly contrasting with the severe sterility (5%–28% viable pollen) observed in controls (Figure 7F). Conversely, RNAi lines (Ri2/Ri7) displayed a 2.5‐4 fold transcriptional suppression of GhAL5 (Figure 7C), while CRISPR‐generated Ghal5 knockout mutants (Ghal5‐ko) (Figure 7B) showed hypersensitivity to heat. These loss‐of‐function lines exhibited compromised seedling thermotolerance and dramatically elevated pollen sterility rates (86%–97% unviable pollen) under stress conditions. The results consistently showed that GhAL5 overexpression lines exhibited enhanced thermotolerance (Figure 7E,F), while the RNAi lines and knockout mutants displayed diminished thermotolerance, demonstrating that GhAL5 is an excellent positive regulator of thermotolerance in cotton.
Remarkably, heterologous GhAL5 expression in rice recapitulated this thermoprotective role. After confirming positive transformation and expression levels of GhAL5 in rice (Figure 7G,H), we subjected the transgenic rice plants (GhAL5‐OE10/OE27) to HT treatment. The results showed that overexpressing of GhAL5 significantly enhanced rice plants HT tolerance (Figure 7I). Meanwhile, virus‐induced gene silencing (VIGS) and transient overexpression of CaAL5 in pepper seedlings respectively decreased and increased thermotolerance at the seedling stage (Figure S28). These findings suggest that the chromatin‐regulated GhAL5 can act as a heat protector in various crops.
To further characterize the 3D genomic organization in anthers of GhAL5 genetic materials respond to HT, comparative analysis using the wild‐type as control revealed that GhAL5 RNAi lines exhibited a higher frequency of A/B compartment transitions (49.35 Mb) compared to overexpression lines (30.85 Mb) (Figure 7J). Notably, the majority of compartment conversions in GhAL5 RNAi lines involved transformation from B to A compartments (90.4%) (Figure 7J), this parallels our observation that excessive B to A transitions in heat‐sensitive line H05 correlate with maladaptive transcriptional activation (Figure 3E,F), suggesting GhAL5 fine‐tunes genome topology to prevent stress‐induced chromatin overactivation. Collectively, these results position GhAL5 as a central regulator that couples 3D chromatin reconfiguration to transcriptional resilience. Its conserved function across cotton and rice underscores its potential as a universal thermotolerance gene, offering a dual mechanism‐direct transcriptional regulation and genome topology modulation‐to combat heat‐induced male sterility.
Discussion
3
The genetic mechanisms underlying crop thermotolerance are complex, involving the coordinated action of numerous genes and interactions with other environmental factors [48, 49, 50]. Deepening our understanding of the thermotolerance mechanisms in crop and pinpointing key thermotolerance genes is an effective approach to address the scarcity of thermotolerance resources and to improve the efficacy of thermotolerance‐oriented breeding programs. The recent advancements and applications of phenomics, genomics, transcriptomics, and metabolomics platforms have provided a robust framework for advancing thermotolerance breeding [48, 51, 52]. In this study, we successfully mapped the dynamic 3D genomic structures of cotton anthers under high‐temperature (HT) stress, revealing a critical connection between chromatin topology and male fertility. The findings underscore the importance of chromatin remodeling in regulating gene expression and plant thermotolerance, offering valuable insights into the molecular mechanisms that enable certain cotton lines to maintain fertility under HT stress.
Polyploidization, characterized by the doubling of chromosome sets within a cell nucleus and their heritable transmission to offspring, has been a recurrent event throughout plant evolution [53, 54]. Polyploids provide opportunities for plants to adapt to environmental stresses and diversify their species, serving as a key driver of plant evolution and biodiversity. Polyploidization is also regarded as a key factor in the domestication of crops such as wheat, rapeseed, soybeans, and cotton [55, 56]. For instance, allopolyploid cotton, with its AD subgenomes, has demonstrated enhanced flexibility in responding to moderate salt stress compared to its diploid progenitors, exhibiting novel phenotypic adaptations [57]. Studies on the global transcriptional profiles in allopolyploid plants have revealed expression biases among homoeologous genes, potentially leading to functional diversification [58, 59, 60]. In hexaploid wheat, the coordinated expression of homeologs from the A, B, and D subgenomes under long‐term salt stress has been linked to the plant's robust salt tolerance through a dosage effect [61]. Similarly, in cotton fibers, the simultaneous selection and utilization of homoeologous gene pairs from both subgenomes can enhance fiber length and strength [62]. Our comparative analyses of the HT‐tolerant line 84021 and the HT‐sensitive line H05 provide compelling evidence that subgenome‐specific regulation is tightly linked to thermotolerance, not only at the bulk tissue level, but also in a cell‐type‐specific manner. The relatively unbiased expression of homoeologous genes observed in 84021‐particularly in heat‐sensitive cell types such as tapetal cells and microspores‐suggests a model in which subgenomic expression balance contributes to cellular resilience under HT stress (Figure 1E–G). These results align with previous studies showing that polyploidy provides evolutionary advantages to plants under stress by promoting regulatory flexibility between subgenomes.
In previous studies, the rate of homoeologous chromosome exchanges, transposon densities, and specific epigenetic modifications at the subgenomic level have been demonstrated to modulate the expression of homoeologous genes [63, 64]. The chromatin landscape observed in this study highlighted significant differences between HT‐tolerant and HT‐sensitive lines. Notably, the At subgenome exhibited a more pronounced response to HT stress, particularly at the chromatin compartment level, which contributed to the differential regulation of gene expression (Figure 3H). This finding reinforces the idea that the spatial reorganization of chromatin, coupled with dynamic chromatin looping, plays a crucial role in determining a plant's response to abiotic stress. One limitation of our study is that the Hi‐C compartment analysis was conducted on merged replicates, and statistical comparisons between biological replicates could not be performed. Therefore, the observed subgenome‐specific patterns should be interpreted with caution and warrant further validation. Additionally, the increased transcriptional activity at the boundaries of TAD‐like structures suggests that chromatin interactions at these loci are essential for activating or silencing key stress‐responsive genes (Figure 4E). Our identification of GhAL5 as a key thermotolerance regulator further supports the role of chromatin dynamics in stress responses. GhAL5, regulated by changes in TAD‐like structures boundaries, was shown to positively influence male fertility under HT stress. However, unlike conventional heat shock factors that act through linear promoter interactions [36], GhAL5 employs a dual regulatory strategy to regulate crop heat tolerance: (1) stabilizing TAD‐like structures boundaries to prevent maladaptive chromatin decompression (Figure 7J); and (2) facilitating enhancer‐promoter cycling through the histone modification recognition ability of its PHD structural domain. Functional analyses, including overexpression in rice, demonstrate its potential as a heat protector across different plant species, emphasizing the broader applicability of our findings. Importantly, the 3D genome‐level regulatory mechanism employed by GhAL5 opens new possibilities for chromatin‐aware crop improvement strategies. Unlike classical promoters or transcription factors, chromatin regulators can reshape higher‐order genome architecture, allowing coordinated activation of stress‐response networks. Given its dual role in chromatin topology and enhancer dynamics, GhAL5 represents a promising target for genome editing approaches aimed at enhancing thermotolerance in cotton and potentially in other crops.
Moreover, growing evidence from other species suggests that 3D genome reorganization is a conserved mechanism in plant stress responses. For instance, in tomato, the heat shock factor HSFA1a mediates heat‐induced enhancer–promoter interactions to activate defense genes [65]; in apple, the chromatin remodeler MdRAD5B modulates both chromatin accessibility and H3K27me3 to regulate drought tolerance [66]; in rice, OsHMGB1 enhances phosphate starvation responses by increasing chromatin openness [67]. Notably, in Arabidopsis, heat stress has been shown to trigger 3D genome reorganization and transposable element activation in heterochromatic regions [29]. These findings collectively support the idea that dynamic chromatin architecture is an evolutionarily conserved and versatile regulatory layer in plant adaptation to environmental stress. While certain features, such as subgenome coordination in allotetraploid cotton, may be species‐specific, the core principle of chromatin‐based regulation is widely shared and can be harnessed in breeding programs across diverse crops.
Three‐dimensional genome topological structures have provided substantial information for studying genome structures and interactions between cis‐regulatory elements and genes in both animals and plants [10, 23]. Promoters, silencers, enhancers, and other cis‐regulatory elements, as well as trans‐acting factors, can be connected through chromatin loops to regulate gene expression, which is the basis for the involvement of chromatin loops in gene expression regulation [23]. Utilizing in situ Hi‐C technology, we obtained high‐resolution chromatin interaction maps under HTs and revealed numerous chromatin loop topological structures across the genome. Importantly, the dynamic reorganization of chromatin loops under HT stress suggests that long‐range chromatin interactions are integral to gene regulation. Genes connected by G‐G loops showed enhanced expression, while G‐N loops, particularly those linking genes with non‐genic regulatory elements (Figure 6F,K), were implicated in fine‐tuning gene activity. These findings shed light on the complex regulatory networks orchestrated by chromatin loops, providing new targets for genetic manipulation aimed at improving crop resilience. This study shifts the paradigm of abiotic stress research from “stress‐responsive genes” to “stress‐responsive genome architectures”. The finding that many HT‐induced DEGs reside in dynamically reorganized chromatin regions (Table S4) underscores 3D genome engineering's potential for climate‐resilient crop development. Future efforts should explore CRISPR‐based chromatin loop engineeringand synthetic polyploid designs optimizing subgenome compartmentalization‐strategies poised to revolutionize precision breeding in the post‐genomics era.
Conclusion
4
In conclusion, our findings underscore the potential of targeting chromatin structure and dynamics as a strategy for breeding thermotolerant crops. By elucidating the regulatory mechanisms underlying homoeologous gene expression and chromatin reorganization, this study provides valuable insights for improving stress resilience in cotton and other crops. Future efforts should focus on translating these findings into practical breeding programs, leveraging chromatin engineering and genome‐editing technologies to enhance crop performance under increasingly variable environmental conditions.
Experimental Section
5
Plant Materials and HT Treatment
5.1
The cotton varieties 84021 (heat‐tolerant) and H05 (heat‐sensitive) exhibited distinct responses to high‐temperature (HT) stress. HT treatment during the flowering stage of cotton was initiated when daytime temperatures exceeded 35°C and nighttime temperatures exceeded 28°C. After seven days of exposure, cotton anther samples were collected at three key developmental stages: the tetrad stage (TS), tapetum degradation stage (TDS), and anther dehiscence stage (ADS). All samples were immediately snap‐frozen in liquid nitrogen and stored at −80°C for subsequent analyses. For seedling‐stage HT treatment, rice and cotton seedlings at the four‐leaf stage were subjected to 42°C and 48°C, respectively, for 48 h, followed by recovery under normal conditions for phenotypic documentation.
Genome Collinearity and Homoeologous Gene Analysis
5.2
The identification of collinear blocks between the At and Dt subgenomes was performed using MCScanX software [68]. This involved an all‐vs‐all blastp search across all genes, with parameters set to ‐e 1e‐10, ‐v 5, and ‐b 5. Collinear blocks were defined as genomic regions containing at least five collinear genes. Additionally, a reciprocal all‐vs‐all blastn alignment was carried out on the At and Dt subgenomes to identify homoeologous genes, using the coding sequences (CDS) of each gene (parameters: ‐e 1e‐10, ‐v 1, ‐b 1). The best reciprocal matches within collinear blocks were considered as homoeologous gene pairs.
RNA‐seq and ChIP‐seq Data Analysis
5.3
For RNA‐seq analysis, clean reads were aligned to the upland cotton reference genome using HISAT2 (version 2.1.0) [69] with default settings. Reads with mapping quality scores above 25 were used to assemble transcripts through StringTie [70] (parameters: –fr ‐e ‐G), and gene expression was quantified. In this study, genes were considered expressed if their FPKM (fragments per kilobase of transcript per million mapped reads) values exceeded 1 in all three biological replicates. For ChIP‐seq analysis, clean reads were mapped to the reference genome using Bowtie2 (version 2.2.4) [71] with parameters set to ‐N 1 and ‐L 30. After removing PCR duplicates, reads with high mapping quality (score >25) were used for peak calling with the MACS software (version 2.2.7.1) [72], using input DNA as the control and parameters set to ‐q 0.01.
In Situ Hi‐C Experiment and Sequencing
5.4
Anther samples were ground into powder in liquid nitrogen and fixed using NI buffer containing 1% formaldehyde for 15 min at room temperature. The fixation was stopped by adding glycine. The fixed materials were filtered through miracloth and a 20 µm nylon membrane to collect pure nuclei. These nuclei were treated with NEB 3.1 buffer and SDS to increase membrane permeability, followed by digestion with the DpnII restriction enzyme. Biotin was used for labeling, and DNA ligase was added for proximal ligation. After nuclei lysis and chromatin de‐crosslinking, DNA was isolated. DNA fragments of 300–500 bp were selected for library construction. Biotin‐labeled DNA fragments were purified using Dynabeads MyOne Streptavidin T1 beads (Invitrogen; 65602), and libraries were prepared with the VAHTS Universal DNA library construction kit for MGI (Vazyme, NDM607). Sequencing was performed on the MGI2000 platform, with two biological replicates. The sequencing depths estimated by raw data were 187× (H05‐HT‐ADS), 214× (H05‐HT‐TDS), 195× (H05‐HT‐TS), 169×(H05‐NT‐ADS), 194×(H05‐NT‐TDS), 181× (H05‐NT‐TS), 173× (H05‐HT‐ADS), 177×(H05‐HT‐ADS), 179×(H05‐HT‐ADS), 209×(H05‐HT‐ADS), 210×(H05‐HT‐ADS) and 227×(H05‐HT‐ADS).
Hi‐C Sequencing Data Processing
5.5
The HiC‐Pro pipeline [73] was employed to process Hi‐C data. Clean reads were aligned to the upland cotton genome using Bowtie2. Low‐quality reads, multiple hits, PCR duplicates, and other artifacts were removed. The iterative correction and eigenvector decomposition (ICE) method was used to normalize the Hi‐C contact matrices at varying resolutions (5, 10, 20, 40, and 200 Kb) for downstream analyses.
Interaction Frequency Heatmaps and Visualization
5.6
Heatmaps displaying chromatin interactions were generated using HiC‐Pro‐produced matrices, along with HiCPlotter and Juicebox tools [73, 74, 75]. Chromatin interactions were categorized into short‐range (≤2 Mb) and long‐range (>2 Mb) based on genomic distance. Chromatin compaction was assessed by counting reads that interacted with each bin within a 1 Mb surrounding window at a 10 Kb resolution.
Contact Probability Analysis
5.7
Contact intensity, representing interaction probability between two bins, was calculated as log10 (sum of observed reads between bins/average reads) at a 10 Kb resolution, with distances ranging from 10 Kb to 100 Mb.
A and B Compartment Analysis
5.8
The hicPCA program from HiCExplorer was used to identify A/B compartments at a 40 Kb resolution [74]. Positive first principal component values, corresponding to high gene density regions, indicated A compartments, while negative values marked B compartments. Transitions between A and B compartments were identified when at least two consecutive bins switched between positive and negative values.
Topologically Associated Domain (TAD‐Like) Structure Analysis
5.9
TAD‐like structures were identified using the hicFindTADs program in HiCExplorer at a 20 Kb resolution [74]. Insulation scores were calculated with the same tool. Conserved TAD‐like boundaries between two developmental stages were defined as boundaries that shifted by no more than two resolution distances (40 Kb). Homologous TAD‐like structures between subgenomes were those containing the same homoeologous genes, while partitioned TAD‐like structures contained differing homoeologous genes.
Chromatin Loop Analysis
5.10
Chromatin loops (FDR < 0.005, contacts > 10) were identified using FitHiC2 at a 5 Kb resolution [76]. To account for chromatin length differences between gene‐poor (≤20 genes in 500 Kb) and gene‐rich regions (>20 genes in 500 Kb), loop anchor ratios were normalized. Chromatin loop densities were also normalized by the total number of loops across different periods.
Statistical Analysis
5.11
Differentially expressed genes were determined using DESeq2 [77] with a negative binomial distribution model. Genes with fold changes greater than 2 and FDR values less than 0.01 were considered significant. Gene Ontology (GO) enrichment was performed using Fisher's exact test, with GO terms having FDR < 0.05 considered enriched. All statistical significance tests in the manuscript, including t‐tests and Wilcoxon signed‐rank tests, were performed using R.
ATAC Experiment and ATAC‐seq Data Analysis
5.12
ATAC‐seq libraries were generated following a modified protocol from a previous study [78]. Intact nuclei were isolated from frozen anther samples, and after treatment with Triton X‐100, 10,000 nuclei were subjected to transposition using Tn5 transposase (VAHTS, TD501). The resulting DNA fragments were purified and amplified for 10–13 PCR cycles to construct libraries, which were size‐selected using DNA clean beads (KAPA) and sequenced on the Illumina Hiseq X‐Ten platform. Raw sequencing data were processed using Trimmomatic and mapped to the genome with Bowtie2 [71, 79]. PCR duplicates and low‐quality reads were removed, and peaks were called using MACS2 with specific parameters (–shift ‐100 –extsize 200 –nomodel ‐B ‐SPMR) [63].
Plant Materials, Vector Construction, and Genetic Transformation
5.13
The cotton variety G. hirsutum cv Jin668 was used for genetic transformation experiments. RNAi vectors for GhAL5 were constructed in the pHellsgate 4 vector, while overexpression constructs were made in the pK2GW7 vector. CRISPR/Cas9 constructs were created as previously described [80]. Transgenic plants were generated via Agrobacterium‐mediated transformation [81], and GhAL5 expression levels were assessed by qRT‐PCR using GhUB7 as the internal control. Mutations in CRISPR/Cas9 lines were analyzed using the Hi‐TOM method [82]. The rice transgenic work was carried out by Wuhan Tianwen Biotechnology Co., Ltd., and the transgenic rice recipient was ZH11. All primer sequences used for vector construction are listed in Table S7.
Protein Extraction and Western Blot Analysis
5.14
Total protein was extracted from cotton anthers using a phenol‐based method as previously described [83]. Proteins were separated via SDS‐PAGE, transferred to a PVDF or nitrocellulose membrane, and blocked with 5% milk in TBST. The membrane was incubated with a primary antibody specific to MYC‐AL5, followed by a secondary HRP‐conjugated antibody. Detection was performed using an ECL substrate and a chemiluminescent imaging system.
Author Contributions
Y.L.L., J.Y., and W.W. contributed equally to this project (co‐first author). L.M. and J.K. are both co‐last authors. L.M. and J.K. designed and supervised the research, Y.L.L. carried out the experiments and wrote the main manuscript text, J.Y., W.W., L.P., C.Z., Y.M., R.Z., Y.F., Y.W.L., H.M., R.L. performed anther sampling. X.Z., L.Z., and S.J. participated in discussion. L.M., J.K., and X.Z. revised the manuscript. All authors reviewed the manuscript.
Conflicts of Interest
The authors declare no conflicts of interest.
Ethics Approval and Consent to Participate
Not applicable
Supporting information
Supporting File 1: advs73612‐sup‐0001‐SuppMat.docx.
Supporting File 2: advs73612‐sup‐0002‐Table S1.xlsx.
Supporting File 3: advs73612‐sup‐0003‐Table S2.xlsx.
Supporting File 4: advs73612‐sup‐0004‐Table S3.xlsx.
Supporting File 5: advs73612‐sup‐0005‐Table S4.xlsx.
Supporting File 6: advs73612‐sup‐0006‐Table S5.xlsx.
Supporting File 7: advs73612‐sup‐0007‐Table S6.xlsx.
Supporting File 8: advs73612‐sup‐0008‐Table S7.xlsx.
Supporting File 9: advs73612‐sup‐0009‐Table S8.xlsx.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1M. R. Anwar , B. Wang , L. Liu , and C. Waters , “Late Planting Has Great Potential to Mitigate the Effects of Future Climate Change on Australian Rain‐fed Cotton,” Science of The Total Environment 714 (2020): 136806, 10.1016/j.scitotenv.2020.136806.31982770 · doi ↗ · pubmed ↗
- 2Y. Li , Y. Li , Q. Su , et al., “High Temperature Induces Male Sterility via MYB 66–MYB 4–Casein Kinase I Signaling in Cotton,” Plant Physiology 189 (2022): 2091–2109, 10.1093/plphys/kiac 213.35522025 PMC 9342968 · doi ↗ · pubmed ↗
- 3L. Min , Y. Li , Q. Hu , et al., “Sugar and Auxin Signaling Pathways Respond to High‐temperature Stress during anther Development as Revealed by Transcript Profiling Analysis in Cotton,” Plant Physiology 164 (2014): 1293–1308, 10.1104/pp.113.232314.24481135 PMC 3938621 · doi ↗ · pubmed ↗
- 4A. H. Khan , L. Min , Y. Ma , M. Zeeshan , S. Jin , and X. Zhang , “High‐Temperature Stress in Crops: Male Sterility, Yield Loss and Potential Remedy Approaches,” Plant Biotechnology Journal 21 (2023): 680–697, 10.1111/pbi.13946.36221230 PMC 10037161 · doi ↗ · pubmed ↗
- 5A. V. Probst , E. Dunleavy , and G. Almouzni , “Epigenetic Inheritance during the Cell Cycle,” Nature Reviews Molecular Cell Biology 10 (2009): 192–206, 10.1038/nrm 2640.19234478 · doi ↗ · pubmed ↗
- 6D. Sitbon , K. Podsypanina , T. Yadav , and G. Almouzni , “Shaping Chromatin in the Nucleus: the Bricks and the Architects,” Cold Spring Harbor Symposia on Quantitative Biology 82 (2017): 1–14, 10.1101/sqb.2017.82.033753.29208640 · doi ↗ · pubmed ↗
- 7T. Yadav , J. P. Quivy , and G. Almouzni , “Chromatin Plasticity: a Versatile Landscape That Underlies Cell Fate and Identity,” Science 361 (2018): 1332–1336, 10.1126/science.aat 8950.30262494 · doi ↗ · pubmed ↗
- 8F. X. Wang , G. D. Shang , L. Y. Wu , Z. G. Xu , X. Y. Zhao , and J. W. Wang , “Chromatin Accessibility Dynamics and a Hierarchical Transcriptional Regulatory Network Structure for Plant Somatic Embryogenesis,” Developmental Cell 54 (2020): 742–757, 10.1016/j.devcel.2020.07.003.32755547 · doi ↗ · pubmed ↗