Phenotypic differences under drought resistance of two cotton varieties

We selected two varieties, ZY007 (drought-sensitive) and ZY168 (drought-tolerant), that showed significant differences in performance under drought treatment, and observed their plant phenotypes at different stages of drought. On the fourth day of drought treatment (Initial Drought, ID), the water content of the treated group gradually decreased, and the edges of the leaves of both varieties softened slightly. On the sixth day of drought treatment (Mild Drought, MD), the water content of the treated group continued to decrease, and the leaves of both varieties wilted completely. On the ninth day of drought treatment (Severe Drought, SD), the water content of the treated group further decreased, and the true leaves of both varieties gradually wilted. At this time, the drought-tolerant variety ZY168 showed better growth than the drought-sensitive ZY007, as evidenced by more severe wilting of the leaves of ZY007. On the 11th day of drought treatment, the treated group was saturated with water, and after 24 h of watering (Re-water, RW), ZY168 gradually returned to normal, while ZY007 remained in a state of continued wilting (Fig. 1a,b).

Fig. 1

Phenotypes and i-trait (PD6) of drought-tolerant and drought-sensitive varieties’ seedlings under different drought stresses conditions. a The dynamic phenotypes of drought-tolerant and drought-sensitive variety seedlings under different drought treatments. ID (Initial Drought), MD (Mild Drought), and SD (Severe Drought) represent the fourth day, sixth day, and ninth day after drought stress, respectively. RW (re-water) represents the second day after re-watering. Scale bar is 10 cm. b Average weight per pot at different time-points after drought treatment. The gray lines represent the four time-points. C represents normal water supply condition; D represent drought stress condition, respectively. c Modeling color maps of PDs of drought-tolerant and drought-sensitive varieties under MD and SD. Scale bar is 10 cm. d PD6 of drought-tolerant and drought-sensitive varieties. Higher PD6 represents more severe wilting status

At the same time, we created pseudocolor maps of the plants from both varieties under different drought conditions (Fig. 1c), and PD6 (Plant Density dimension 6) was determined to be a new indicator for plant drought severity. The higher the PD6, the more severe the wilting of the plant [21]. We found that under SD conditions, the PD6 of ZY007 was higher, indicating more severe wilting and higher sensitivity to drought. After rehydration, the PD6 of ZY168 decreased more, indicating a faster recovery from wilting and stronger post-drought resilience (Fig. 1d).

Construction of drought-induced expression landscape in cotton

In order to explore the similarities and differences of the drought resistance of the two cotton varieties at the transcriptional level, we sampled the leaves of the drought-treated plants and control plants at the above four treatment stages (ID, MD, SD, and RW). A total of 32 samples were obtained (two biological replicates), and transcriptome sequencing was performed on the samples. After preliminary analysis, we collected a total of 2.88 billion valid reads (Additional file 2: Table S1), and the correlations between biological replicates based on reads count of the bins were between 0.87 and 0.96 (Additional file 1: Fig. S1a).

Subsequently, we calculated the expression levels of all genes and constructed an expression matrix. The correlations of gene expression between two biological replicates were between 0.96 and 0.99, which allowed us to merge two biological replicates (Additional file 1: Fig. S1b). Cluster analysis showed that all samples were grouped into four clusters. The first cluster included controls and the ID stage, the second cluster included controls and RW stage of ZY168, the third cluster included the MD stage, and the fourth cluster included the SD stage and the RW stage of ZY007 (Additional file 1: Fig. S1c). From these results, it can be inferred that different sampling time points under drought can be effectively distinguished at the transcript level. Moreover, at the transcript level, the RW samples of the drought-sensitive variety ZY007 still exhibited transcriptional expression similar to that of the SD stage samples, while the drought-tolerant variety ZY168 showed significant recovery in gene expression after rewatering, which was close to that of the ID stage.

Compared to the control, the number of expressed genes decreased with prolonged and intensified drought (Additional file 1: Fig. S2a). Subsequently, we conducted pairwise comparisons of control and drought-treated samples to identify differentially expressed genes (DEGs), resulting in 2619, 10,353, 19,460, and 17,582 DEGs in the four time points, respectively (Additional file 1: Fig. S2b; Additional file 2: Table S2). Except for the ID stage, the number of up-regulated genes was lower than that of down-regulated genes, indicating an overall decrease in gene expression levels under drought. We performed GO enrichment analysis on these DEGs, and the results showed significant enrichment in pathways related to plant stress response, such as ethylene, abscisic acid, jasmonic acid, kinase activity, water transport, and redox reactions (Additional file 1: Fig. S2c). These results illustrate a cotton expression profile in response to drought stress.

Asymmetric expression of subgenomes under drought induction

To investigate the differences between two varieties in response to drought stress, we defined drought-induced genes (see “Methods”). The results showed that in ZY007, the gene expressions levels of a total of 22,606 genes were changed in response to drought, including 8601 up-regulated and 14,005 down-regulated genes. In ZY168, the gene expressions levels of a total of 19,552 genes were changed in response to drought, including 6857 up-regulated and 12,695 down-regulated genes (Fig. 2a). Analysis of the expression patterns of differentially expressed genes shared between the two varieties revealed that compared to ZY007, ZY168 exhibited stronger expression recovery ability after rehydration. Both up-regulated and down-regulated genes showed more expression-recovered trends after rehydration in ZY168 (Fig. 2b,c). After rehydration, the numbers of recovered and differentially expressed drought-induced genes were 8109 and 14,497 in ZY007 and 3326 and 16,226 in ZY168, respectively (Chi-squared test, P ≤ 2.2 × 10–26). This suggests that ZY168 has stronger recovery ability at the expression level than ZY007 consistent with phenotype.

As an allopolyploid crop, cotton exhibits a certain degree of asymmetry between subgenomes [22]. To investigate whether the response of the two subgenomes has asymmetry under drought, we performed the following three types of analysis. Firstly, we compared the expression changes of homoeologous genes in the two subgenomes relative to controls under drought treatment. Based on the direction of expression changes, the homoeologous genes were divided into three categories: common, single, and opposite (see “Methods”). The results show that the largest category was single, followed by common, and only a very small proportion (less than 0.9%) of homoeologous genes showed differential expression in opposite directions (Fig. 2d). Moreover, as the degree of drought intensifies, the proportion of single genes decreases while that of common genes increases. This trend is restored after rehydration. These findings suggest that under drought conditions, homoeologous genes exhibit convergent expression patterns.

Fig. 2

Expression patterns of drought-induced genes and expression preferences among homoeologous genes. a The number of drought-induced genes in the two varieties, which are divided into up/down-regulated and At/Dt subgenomes. b–c Expression patterns of genes in the intersection part of the two varieties at different stages. d The genes were divided into three groups according to their expression direction of homoeologous genes in the two subgenomes under drought stress treatment compared to the control. e The numbers and proportions for the three types of homoeologous genes with bias expression under four conditions. Red, blue, and gray colors represent stable At bias expression, stable Dt bias expression, and dynamic bias expression, respectively (two-sided Chi-square test, ***P < 0.001). f Number of genes with changed expression bias induced by drought. g Genome-wide bias values for two varieties at four drought stress stages. C and D represent normal water supply and drought stress condition, respectively

Secondly, we compared the bias expression patterns of homoeologous genes during different stages. We identified a total of 23,703 subgenome homoeologous gene pairs. Among them, 8121 homoeologous gene pairs with bias expression of subgenome were identified in ZY007 under control conditions, 8789 in ZY007 under drought, 8169 in ZY168 under control conditions, and 8925 in ZY168 under drought. Analysis of the dynamic bias expression patterns of these homoeologous gene pairs during four stages (see “Methods”) showed that homoeologous gene pairs with dynamic bias expression between four stages accounted for the majority (Fig. 2e). Furthermore, the proportion of homoeologous gene pairs exhibiting dynamic bias expression increased significantly under drought and obviously decreased after rehydration in both varieties (Additional file 1: Fig. S2d).

Thirdly, to investigate the effects of drought stress on bias expression of subgenome, we identified 7778 and 7516 genes with expression preference changes induced by drought in ZY007 and ZY168, respectively (Fig. 2f). Approximately 62% of these genes were shared by both varieties. To further explore the differences between the two varieties, we defined a bias value (see “Methods”). Results from the calculation of the bias values for the whole genome at four stages indicated that the genome exhibited an expression bias of At subgenome under drought (Fig. 2g). The bias value gradually increased for ZY007 and showed a trend of first increasing and then decreasing for ZY168 (Fig. 2g). To investigate the cause of the bias expression of the At subgenome, we conducted detailed analysis of the gene expression levels. The results indicated that under drought stress, the expression levels of both subgenomes decreased, but the decrease was greater for the Dt subgenome (Additional file 1: Fig. S2e). This phenomenon was particularly evident of SD in ZY007, which corresponded to a significant increase in the bias value. These results indicate that with the increasing severity of drought, expression of the Dt subgenome was more strongly inhibited than At subgenome in drought-sensitive ZY007, while the inhibition difference of both subgenomes in ZY007 was greater than in ZY168.

Construction of drought-related co-expression networks

To further screen and identify drought-related genes from the transcriptome, we constructed a co-expression network of cotton in response to drought using the union of all DEGs, comprising 26,412 genes. A threshold weight value greater than 0.2 was used to filter the co-expression relationships, resulting in 2402 genes and 34,064 co-expression relationships in the network (Fig. 3a). The co-expression network consisted of nine modules, with each module containing a different number of genes (ranging from 7 to 1291). We then identified hub genes for each module by selecting genes with the highest total weight values (Fig. 3b). We further examined the expression patterns of genes in each module in response to drought in two varieties. As a result, we found that genes in module 2 showed significant upregulation under drought in both varieties, followed by obvious recovery after rehydration, with ZY168 exhibiting stronger recovery (Fig. 3c). We applied Gene Ontology (GO) enrichment analysis for each module and found that module 2 was enriched in pathways related to ion balance, jasmonic acid response, abscisic acid response, osmotic pressure regulation, and redox reactions, which are closely related to plant stress responses (Fig. 3d; Additional file 2: Table S3). Module 2 also contained many previously reported key genes related to drought resistance, such as PP2CA [23], RAB28 [24], and HB-7 [25] (Fig. 2b; Additional file 2: Table S4). Overall, module 2 is the co-expression module closely associated with drought stress, and its genes have close connections with drought resistance.

Fig. 3

Co-expression network construction and screening. a The display of the co-expression network, different colors represent different modules, and the size of the hub gene is positively correlated with the sum of the weight values. b A detailed display of the connection of the three hub genes in module 2. c The expression patterns of the genes of module 2 in two varieties, two treatments, and four stages are shown. d GO enrichment results of genes in module 2

A/B compartment switching under drought treatment

In addition to transcriptional regulation, emerging researches underscore the impact of environmental stressors on chromatin architecture [18,19,20]. In order to investigate the changes in chromatin higher-order structure at the 3D genomic level and their potential biological functions in cotton, we performed in situ Hi-C experiments of the 32 samples described above. Ultimately, we obtained 32 Hi-C libraries containing approximately 63.7 billion raw sequences (ranging from 1.37 billion to 5.16 billion). After initial processing, we obtained a total of 12.8 billion valid interactions (ranging from 330 to 480 million; Additional file 2: Table S5). The correlation between the two biological replicates was good (correlation coefficient ranging from 0.93 to 0.97; Additional file 1: Fig. S3a). Therefore, we merged the two biological replicates for further analysis. The resolution of the merged 16 samples can reach 10 kb, with seven samples reaching 5 kb (Additional file 1: Fig. S3b). Chromatin exhibits distinct and obvious structures in the heatmaps at different resolutions, ensuring the feasibility of our subsequent analysis of A/B compartments and TAD (Fig. 4a). Overall, the high-resolution 3D genomic map of cotton in response to drought stress is reliable.

Fig. 4

Switching between A/B compartment. a Heatmap of chromatin interaction in part of regions of A04 at different resolutions at the ID of ZY007 under control treatment as an example of 3D chromatin map. b Proportions of the A/B compartment in two varieties for the four stages under the two treatments. c Global dynamic switching of chromatin compartment status. Dark yellow and light yellow represent the A and B compartment, respectively. Heatmaps show bins with status switching between A compartment and B compartment. AB means switching from A compartment to B compartment (including AAAB, ABAB, AABB, ABBB). BA means switching from B compartment to A compartment (including BAAA, BABA, BBAA, BBBA). ABA means switching from A compartment to B compartment and then from B compartment switching to A compartment (including ABBA, ABAA, AABA). BAB means switching from B compartment to A compartment, and then from A compartment to B compartment (including BAAB, BABB, BBAB). d The figure above shows a region on chromosome D03 that undergoes drought-induced switching from A compartment to B compartment. The figure below shows the changes in the expression of two genes in the switching region

The A/B compartment is an important structure about the chromatin, where the A and B compartment represents active and repressive chromatin state, respectively [26]. We identified the A/B compartments of 16 samples using a 100-Kb resolution interaction matrix. In the control, the A compartment increased and the B compartment decreased as development progressed, while an opposite trend was observed under drought (Fig. 4b). It is noteworthy that in ZY168, the downward trend of the A compartment after rehydration finished, whereas in ZY007, the A compartment continued to decrease after rehydration. Specifically, ZY007 and ZY168 had intervals of 196 Mb and 141 Mb where the A/B compartment switching occurred during drought. From the heatmap, it can also be seen that ZY007 had more regions where the irreversible (not recovered after rehydration) switching from the A compartment to the B compartment occurred (Figs. 4c, S4a). For genes, both varieties showed more genes undergoing switching from the A compartment to the B compartment during drought stress (Additional file 1: Fig. S4b), especially during ID stage (more than 30 genes per Mb underwent switching from the A compartment to the B compartment). Combined with transcriptome data, the gene expression level in the A compartment was significantly higher than that in the B compartment (Additional file 1: Fig. S4c). These genes undergoing compartment switching were more likely to have differential expression than other genes (Chi-squared test, X-squared = 72, P < 2.2 × 10–16). Among the Module 2 genes identified above, we found two genes, LTP1 and LTP3, within the regions undergoing A/B compartment switching (Fig. 4d). The corresponding gene expression levels showed a significant increase following switching from the B compartment to the A compartment under drought stress. In Arabidopsis, many members of the LTP family are induced by drought [27]. These results indicate that under drought induction, some genes are involved in the switching between the A/B compartments, and these switching are related to differential gene expression.

Dynamic TAD landscape of cotton development

TAD is a higher-order chromatin structure that occurs at mega-base resolution, playing a crucial role in maintaining gene expression stability under stress conditions [28]. To explore how TADs change under drought treatment, we identified TAD for 16 samples at 20-Kb resolution. A total of 140,376 TAD structures were detected, with the highest number occurring during SD of ZY168 control treatment (8979) and the lowest during the SD of ZY007 drought treatment (8578; Additional file 1: Fig. S5a). The number of TADs remained relatively constant across different treatments and developmental stages. The expression levels of genes located at TAD boundaries were significantly higher than those within TADs in all samples (Additional file 1: Fig. S5b). Next, we performed iterative comparisons among different stages to construct a “pan-TAD boundary” landscape for four stages in two varieties under two treatments (see “Methods”). After removing redundancy, the pan-TAD boundary sets for the four stages ranged from 11,001 to 11,192. We classified TAD boundaries into four categories: stage-specific, conserved in two stages, conserved in three stages, and conserved in all stages. The percentages of TAD boundaries conserved in four stages were the highest, while the percentages of TAD boundaries conserved in two stages was the lowest (Fig. 5a).

Fig. 5

Comparison of TAD boundaries in different stages. a The proportion of the four-type TAD boundaries. The green represents stage-specific TAD boundaries, and the yellow, pink, and orange represent TAD boundaries conserved in two, three, and four stages, respectively. b TAD separation score of the regions where the boundaries of the four types of TADs are located (200 Kb left and right). c An example of a dynamic TAD boundary between stages. The heatmap shows the interaction strength of chromosome 22.8 Mb to 24.0 Mb of A10. The darker the color, the stronger the interaction. The figure below shows the TAD-separation score in this region. d The expression dynamics of the genes contained in the four types of TAD boundary regions in the four stages, expression dynamics was expressed using the logarithm of the range of gene expression (FPKM)

To investigate the characteristics of above four types of TAD boundaries, we examined the insulating index (inferred as TAD-separation score) of each 200-Kb region to the left and right of the TAD boundary. The results showed that more conserved TAD boundaries had lower TAD-separation score (indicating greater stability), while stage-specific TAD boundaries were less stable (Fig. 5b,c). Finally, we combined transcriptome data to examine the differences in gene expression within these four categories of TAD boundary regions. The results showed that genes located within TAD boundaries conserved across all four stages had the smallest expression differences among the four stages, whereas genes within dynamic TAD boundaries showed larger expression differences over time (Fig. 5d). This indicates that the higher-order chromatin structure of plants has a certain degree of plasticity across different developmental stages and disruption of TAD boundaries are associated with dynamic changes in gene expression.

Disruption of TAD structure induced by drought

To investigate how TAD boundaries change under drought stress, we compared TAD boundaries between control and drought-treated conditions. We classified TAD boundaries into three categories: TAD boundaries conserved under drought (conserved), TAD boundaries formed under drought (dg: drought gain), and TAD boundaries lost under drought (dl: drought loss; see “Methods”). Among these three types of boundaries, conserved TAD boundaries were the most prevalent (85.9% to 90.74%), and their proportion decreased with increasing severity of drought stress, reaching the lowest point in SD (Fig. 6a). After rehydration, the proportion of conserved TAD boundaries increased from 85.9 to 89.4% in ZY168, higher than ZY007 (from 87.5 to 89.2%). These results indicate that ZY168 had a stronger drought resistance reflected at the level of chromatin structure.

Fig. 6

TAD changes under drought induction. a Proportions of the three types of TAD boundaries relative to the control in the four stages of the two varieties (dg drought gain, dl drought loss). b The ratio of the three types of TAD boundaries in dynamic level. Green represents stage-specific TAD boundaries; yellow, pink, and orange represent conserved in two, three, and four stages, respectively. c Number of two TAD disruption events (TAD fusion and neo-TAD). d The number of DEGs in the TAD boundary region related to TAD disruption events. Here, 10 Kb on each of the left and right sides of the TAD boundary are designated as the TAD boundary region. e A drought-induced TAD fusion event. The heat map shows the chromatin interaction frequency in this region, and the darker the color, the higher the interaction frequency. The yellow highlight shows a drought-induced up-regulated gene within the TAD boundary

We investigated the characteristics of the above three types of TAD boundaries. Firstly, we observed that TADs with conserved TAD boundaries were significantly larger than the TADs with dg or dl boundaries (Additional file 1: Fig. S5c). In addition, more than 85% of conserved TAD boundaries under drought were also conserved between stages (87.8% for ZY007 and 85.8% for ZY168), whereas only 21.9% to 31.1% of dg and dl boundaries were conserved between stages (Fig. 6b). This indicates that dg and dl boundaries only occur during specific stages rather than being conserved between stages. We observed the distribution of these three types of TAD boundaries throughout the whole genome and found, interestingly, that dg was preferentially located in telomere regions during the ID and MD stages, while dl was preferentially distributed in chromatin arms. This distribution preference changed during SD and RW (Additional file 1: Fig. S5d).

We focused on two important TAD change events: TAD fusion and Neo-TAD (see “Methods”) [29]. We identified a total of 6597 TAD change events, including 3407 TAD fusions events and 3190 Neo-TAD events (Additional file 2: Table S6). Notably, both varieties exhibited significantly more TAD fusion events in SD than in the MD (Fig. 6c). After rewatering, ZY168 showed significantly more Neo-TAD events than ZY007 (506 vs. 421, Fig. 6c). In conjunction with gene expression, we found that both types of TAD changes had a significantly higher number of DEGs than expected, indicating that these two important TAD changes had some impact on gene expression (Fig. 6d). Except for RW, the At subgenome had fewer TAD change events than the Dt subgenome in the other three time points, which could be associated with the stronger expression suppression of the Dt (Additional file 1: Fig. S5e). Combining with module 2 previously screened as significantly related to drought resistance, we identified a gene HB-7 that is located at a dg boundary following a TAD fusion event (Fig. 6e). This gene is also up-regulated in response to drought and is a core gene in module 2 (Fig. 3b). In Arabidopsis, this gene encodes a transcription factor of the HOMEOBOX family, which is dependent on abscisic acid for transcriptional regulation of downstream genes and plays a role in the signal transduction pathway of drought response [25]. All of these findings suggest that TAD structures undergo certain changes under drought induction, and these changes may be associated with the transcriptional regulation of drought-related genes.

We investigated the dg and al hotspots between two varieties under drought (see “Methods”). We identified a total of 165 and 212 dg and dl hotspots in the two varieties, respectively, with an average of 8.8 and 9.2 changed TAD boundaries per hotspot (Additional file 1: Fig. S5f; Additional file 2: Table S7). With the exception of the ID stage, ZY168 had more dg hotspots than ZY007 in all other stages. Notably, dl hotspot numbers showed a significant increase during the SD and decreased significantly after rehydration in both varieties. However, the reduction of dl hotspots was greater in ZY168 (from 73 to 13) than in ZY007 (from 54 to 31; Additional file 1: Fig. S5f). Analysis of the two varieties revealed that 19.4% of dg hotspots (32/165) and 33.0% of dl hotspots (70/212) were shared. These findings suggest that 3D structures of the two varieties were greatly disrupted in SD following a large number of TAD boundary loss, and after rehydration, the TAD structure of drought-tolerant ZY168 was more resilient than drought-sensitive ZY007.

Screening of drought-induced genes associated with 3D genomic variation

To facilitate the study of functional genes related to cotton drought resistance, we summarized and screened for a set of genes (referred to as target genes) that are both related to 3D chromatin structure and induced by drought stress. We employed the following filters to identify genes of interested supported by three pieces of evidence: (1) differentially expressed genes induced by drought stress in the two varieties (ZY007 and ZY168, with 22,606 and 19,552 genes, respectively, Fig. 2a); (2) genes in module 2 (506 in total, Fig. 3a); and (3) genes associated with TAD boundary changes occurring under drought induction (specifically referring to the TAD boundaries of dg and dl, with 13,641 genes in ZY007 and 15,436 in ZY168). We then took the intersection of these three gene sets in each variety separately to obtain the target genes for each variety. Finally, we compared the target genes between the two varieties and took the intersection of the target genes from both varieties.

Using our filters, we identified 92 and 98 target genes that were supported by multiple pieces of genomic evidence in ZY007 and ZY168, respectively (Fig. 7a). Among these target genes, 59 were present in both varieties, while ZY007 and ZY168 had 33 and 39 specific target genes, respectively (Fig. 7b). These target genes include several stress-responsive genes such as PP2C [30], GST8 [31], PR, and HSP20 [32]; protein kinases and phosphatase genes involved in stress response such as CIPK6 [33], PFK4 [34], and HAI2 [10]; transport-related genes such as PDR, NRT2.4 [35], and ANN5 [36]; and 17 TF encoding genes such as TGA4 [37], WRKY23 [38], and HB-7 [25] (Additional file 2: Table S8). Overall, we identified a series of genes both related to 3D chromatin structure and induced by drought stress, which are closely related to cotton response to drought stress.

Fig. 7

A filter for genes associated with 3D chromatin structure and induced by drought. a The left and right figures represent the results of genes of ZY007 and ZY168, respectively. The blue, red, and pink circles represent drought-induced genes, drought-altered TAD boundary-related genes, and module 2 genes, respectively. b The intersection of the target genes in the two varieties. The purple and yellow circles represent the target genes of ZY007 and ZY168, respectively