Parental legacy versus regulatory innovation in salt stress responsiveness of allopolyploid cotton (Gossypium) species.

Polyploidy provides an opportunity for evolutionary innovation and species diversification, especially under stressful conditions. In allopolyploids, the conditional dynamics of homoeologous gene expression can be either inherited from ancestral states pre-existing in the parental diploids or novel upon polyploidization, the latter potentially permitting a wider range of phenotypic responses to stresses. To gain insight into regulatory mechanisms underlying the diversity of salt resistance in Gossypium species, we compared global transcriptomic responses to modest salinity stress in two allotetraploid (AD-genome) cotton species, Gossypium hirsutum and G. mustelinum, relative to their model diploid progenitors (A-genome and D-genome). Multivariate and pairwise analyses of salt-responsive changes revealed a profound alteration of gene expression for about one third of the transcriptome. Transcriptional responses and associated functional implications of salt acclimation varied across species, as did species-specific coexpression modules among species and ploidy levels. Salt responsiveness in both allopolyploids was strongly biased toward the D-genome progenitor. A much lower level of transgressive downregulation was observed in the more salt-tolerant G. mustelinum than in the less tolerant G. hirsutum. By disentangling inherited effects from evolved responses, we show that expression biases that are not conditional upon salt stress approximately equally reflect parental legacy and regulatory novelty upon allopolyploidization, whereas stress-responsive biases are predominantly novel, or evolved, in allopolyploids. Overall, our work suggests that allopolyploid cottons acquired a wide range of stress response flexibility relative to their diploid ancestors, most likely mediated by complex suites of duplicated genes and regulatory factors.

INTRODUCTION

Polyploidy is a prevalent phenomenon among eukaryotes, especially in plants, and is recognized as a driving force in evolution and diversification (Blanc & Wolfe, 2004; Bowman et al., 2007; Wood et al., 2009; Jiao et al., 2011; Wendel, 2015; Van de Peer et al., 2017; Fox et al., 2020; Van de Peer et al., 2020). Polyploidy is also considered a key factor in domestication of crops such as wheat (Triticum aestivum), rapeseed (Brassica napus), cotton (Gossypium spp.), and soybean (Glycine max) (among others; Udall & Wendel, 2006; Salman‐Minkov et al., 2016), as genome doubling is thought to promote diversification and adaptation via various genetic and epigenetic processes (Udall & Wendel, 2006; Renny‐Byfield & Wendel, 2014; Salman‐Minkov et al., 2016; Ding & Chen, 2018). Allopolyploids, which combine two divergent parental genomes (Adams & Wendel, 2005; Doyle et al., 2008; Doyle & Sherman‐Broyles, 2017), often show enhanced characteristics, e.g., better growth vigor (Solhaug et al., 2016), broader ecological adaptation (Shimizu‐Inatsugi et al., 2017), and better stress tolerance (Liu et al., 2012; Yang et al., 2014), compared to their diploid ancestors (Comai, 2005; Fawcett & Maere, 2009; Ni et al., 2009; Chen, 2010; Chen & Birchler, 2013). The increased competitiveness of allopolyploid species in environments that are challenging or unsuitable for the parental species is thought to be attributable to processes including fixed heterozygosity, intergenomic interactions, and expression dosage effects that accompany hybridization and genome doubling (Solhaug et al., 2016).

Allopolyploid success involves accommodating divergent genomes and transcriptomes within the same nucleus (Comai, 2005; Chen, 2007; Leitch & Leitch, 2008; Jackson & Chen, 2010), generating enormous combinatorial complexity and thus the potential for evolutionary novelty. Although additive expression of duplicated genes (homoeologs) may represent an expected outcome for polyploid species (where joint expression of homoeologs is equal to the arithmetic average or summed expression of parental species, i.e., the mid‐parent value), non‐additive expression has been widely documented in a variety of polyploid plants (Flagel & Wendel, 2010; Bardil et al., 2011; Grover et al., 2012; Chelaifa et al., 2013; Yoo et al., 2014; Wendel, 2015; Powell & Doyle, 2017), representing a deviation from vertical transmission of the pre‐existing parental state, also known as parental legacy (Buggs et al., 2014). Examination of non‐additive expression patterns at the whole transcriptome level (Pumphrey et al., 2009; Yoo et al., 2013; Powell et al., 2017; Takahagi et al., 2018; Yu et al., 2018; Peng et al., 2020) led to the recognition of the widely observed phenomena of homoeolog expression bias (HEB), expression level dominance, and transgressive expression of homoeologs (Grover et al., 2012; Yoo et al., 2014; Wu et al., 2018; Sigel et al., 2019; Li et al., 2020). These phenomena collectively represent regulatory novelty that either accompanied or evolved following allopolyploidization, often leading to phenotypic innovation in polyploids relative to their diploid parents. For example, conditional HEB during development and in response to environmental stress, respectively, have been implicated in causing delayed leaf maturation and impaired chlorophyll accumulation in Cucumis × hytivus (Yu et al., 2018) and influencing the response to Fusarium infection (Powell et al., 2017) in allohexaploid wheat. These and other studies support the notion that polyploid success is intimately related to stress responses (Van de Peer et al., 2020), although at present, the patterns and mechanisms by which polyploids alter homoeolog expression to accommodate stressful environments are not yet understood.

Cotton is a major crop with global economic significance as a natural source of textile fiber. Notably, two of the major cotton crop species, Gossypium hirsutum and Gossypium barbadense, are descendants of a single, natural allopolyploidization event that occurred approximately 1–2 million years ago (mya), leading to the evolution of seven modern polyploid species. The evolutionary history of Gossypium L. (the cotton genus) is well known (Wendel & Grover, 2015; Hu et al., 2021) and model diploid progenitors for the allopolyploid clade have been identified (Gossypium arboreum and Gossypium herbaceum, maternal; Gossypium raimondii, paternal; Wendel & Cronn, 2003). Gossypium contains more than 50 recognized species that are extensively distributed worldwide in arid and semi‐arid regions of the tropics and subtropics (Endrizzi et al., 1985; Wendel et al., 2010; Wendel & Grover, 2015; Wang et al., 2018). In general, there exists a geographic distinction in natural distribution of diploid cotton species (mostly found inland) versus polyploids (largely species of littoral environments), which suggests an association between polyploidization and ecological adaptation. However, the underlying genetic basis and molecular processes still remain to be discovered. In terms of the consequences of polyploidy for gene expression, HEB and other phenomena have been globally characterized in various tissues and/or developmental stages of allopolyploid cotton (Rapp et al., 2009; Flagel & Wendel, 2010; Hu et al., 2013; Yoo et al., 2013), but the conditional changes under environmental stress have been minimally surveyed (Dong & Adams, 2011).

As a major abiotic stress, soil salinity limits plant growth and causes substantial yield loss worldwide (Munns & Tester, 2008). Although cotton is classified as a moderately salt‐tolerant crop, soil salinity can still inhibit cotton growth and fiber production, depending upon the genotype and growth stage (Abdelraheem et al., 2019). In a recent survey of salt tolerance among Gossypium species/accessions, we found highly variable responses to a 2‐week treatment of hydroponic salinity (50 and 100 mm NaCl) in allopolyploids relative to the diploid progenitors (Dong et al., 2020). Although the ranking of allopolyploids with respect to salt tolerance is generally intermediate between the more tolerant A‐genome and the less tolerant D‐genome diploids, the allopolyploid species Gossypium mustelinum was most tolerant among all genotypes surveyed. Another interesting finding was that under control conditions, the morphological and physiological traits of allopolyploids surveyed were more like those of the A‐genome diploid progenitor, whereas under salt stress the opposite pattern (D‐like traits) was observed. To further understand the underlying mechanisms of diversity and divergence in salt tolerance following polyploidy, here we investigated the transcriptomic response to salt stress for two allopolyploid cotton species (G. hirsutum and G. mustelinum) relative to their model diploid progenitors, G. arboreum (A‐genome) and G. raimondii (D‐genome). Based on the total expression of homoeologous gene pairs, we conducted differential expression and coexpression analyses in conjunction with our previously generated morpho‐physiological data (Dong et al., 2020). We also studied partitioned homoeolog expression, investigating stress responsiveness in allopolyploid cottons. By comparing both the total and partitioned expression patterns between the salt‐sensitive G. hirsutum and salt‐tolerant G. mustelinum, our study provides new evidence regarding how polyploids gain better tolerance to salt stress.

RESULTS

RNA‐seq analysis of salt tolerance in four Gossypium species

To gain insight into the underlying mechanisms of divergence in salt tolerance among Gossypium species, we implemented RNA‐seq to survey the global transcriptomic response to modest salinity stress (50 mm NaCl) in four cotton species, including two allopolyploids, i.e., G. hirsutum (AD1) and G. mustelinum (AD4), as well as representatives of their diploid progenitors, A‐genome G. arboreum (A2) and D‐genome G. raimondii (D5). All surveyed species were able to survive the salinity treatment while displaying different levels of stress symptoms (e.g., inhibited growth and smaller, thicker leaves) in treated versus control plants. Consistent with our previous findings (Dong et al., 2020), G. mustelinum was most tolerant, followed sequentially by G. arboreum, G. raimondii, and G. hirsutum.

A total of 24 RNA‐seq libraries were obtained with an average of 48 million paired‐end reads sequenced per library, of which 95.9–98.8% were mapped onto the reference genome of G. raimondii (Table S1). For the allopolyploid species, approximately equal proportions of A‐ and D‐homoeolog‐specific reads were assigned (Table S1); in the following analyses, the total expression for each homoeologous pair is referred to as gene expression unless further specified. Principal component analysis (PCA) of the total expression revealed clear clustering of replicated samples, except for one replicate from the G. raimondii control group, which instead clustered with the allopolyploids (Figure 1a); a consistent pattern was also observed for the subgenome partitioned expression (Figure S1b). Hence, this mislabeled G. raimondii sample was excluded from further analysis. The first two principal components (PCs) accounted for around 62% of the total variance: PC1 mainly separated G. raimondii samples from the others, while PC2 was able to separate the stress‐treated and control samples per species except for G. arboreum.

Figure 1

Transcriptomic diversity and salt‐responsive dynamics in cotton. (a) Principal component analysis (PCA) of 24 RNA‐seq samples across diploid and allopolyploid Gossypium species. The mislabeled D5 sample is circled in red. (b) Barplot of differentially expressed genes (DEGs) in response to salt stress treatment in each species. (c, d) UpSet plots of upregulated (c) and downregulated (d) DEGs and their overlap between species. The bottom left horizontal barplot shows the total number of DEGs identified in each species. The top barplot summarizes the number of DEGs for each unique or overlapping combination as indicated by the connecting circles below. (e, f) Gene Ontology (GO) enrichment analysis of upregulated and downregulated DEGs in each species. The top 15 significant, non‐redundant GO terms of each species are shown. The color intensity in both panels denotes −log10 transformed ‘adjusted P‐values’, where darker red/blue indicates a higher degree of enrichment.

Transcriptomic changes in response to salt treatment in both diploid and allopolyploid cottons

A primary goal of this study was to explore the transcriptomic diversity of salt stress responses among the four Gossypium taxa. To characterize the interspecific divergence of salt stress responses, we first applied a multifactor design (species + treatment + species:treatment) to all samples (Table 1). A significant species effect was identified for 29 054 genes, or approximately 80% of the transcriptome, reflecting massive interspecific differences regardless of the salt stress treatment. A significant treatment effect common to all four species was identified for 13 378 (35.9%) genes. Gene Ontology (GO) enrichment analysis revealed that the upregulated genes were primarily involved in oxidoreductase, peptidase, ATPase, and GTPase activity, peptidase activity, and various cellular biosynthetic and catabolic processes, whereas the downregulated genes were primarily enriched in photosynthesis, ion binding, and response to endogenous stimulus (Table S2 and Figure S2). The species:treatment interaction effect reflects species‐specific responses to salt treatment, which was significant for 12 979 (34.9%) genes enriched in GO terms related to ubiquitin‐protein transferase activity, the oxidoreductase complex, and photosynthesis (Table S2).

Table 1

Multivariate analysis of differential gene expression between species and under salt stress treatment

Effect	Differentially expressed	Upregulated	Downregulated
Species	29 054	–	–
Treatment	13 378	6852	6526
Species × treatment	12 979	–	–

We next conducted pairwise comparisons between control and treated samples for each species to further examine the effect of salt treatment. The numbers of differentially expressed genes (DEGs) in response to salt stress per species ranged from 6077 to 11 544 genes, corresponding to 16.2–30.8% of the leaf transcriptome (Figure 1b). The diploid species G. arboreum and G. raimondii exhibited the lowest and highest number of DEGs, respectively, while the allotetraploid species G. hirsutum and G. mustelinum showed intermediate amounts of change between their diploid parents. Intersections between these DEG sets showed that G. arboreum and G. raimondii exhibited the largest number of species‐specific changes for both up‐ and downregulation (Figure 1c,d, respectively). The common set of DEGs shared by the two allopolyploids (1372 upregulated and 756 downregulated) was much larger than the common set shared by two diploids (352 upregulated and 473 downregulated), perhaps reflecting the approximately fivefold greater time since divergence of the diploids than since divergence of the two polyploid species (Wendel & Grover, 2015; Hu et al., 2021). Interestingly, the number of shared DEGs between both allopolyploid cottons and G. raimondii (1427 upregulated and 1316 downregulated) was remarkably higher than the number of DEGs shared with G. arboreum (133 upregulated and 156 downregulated). Only a small set of DEGs (219 upregulated and 273 downregulated) overlapped in all four species. Enriched GO terms of pairwise DEG sets were cataloged (Table S3; Figures S3 and S4) and compared to assess whether the associated biological functions are common or species‐specific (Figure 1e,f). For example, photosynthesis was upregulated in G. arboreum yet downregulated in the other species, while the response to hormones was uniformly enriched in the downregulated DEGs. Notably, no GO term was enriched in upregulated DEGs across species. These pairwise comparisons, together with the multivariate analysis results, suggest that the transcriptional response to salt stress and the involved cellular processes differed significantly among cotton species.

Coexpression gene network analysis uncovers divergent responsiveness to salt stress

To further investigate global transcriptional organization underlying salt responses and interspecific divergence, we applied weighted gene coexpression network analysis (WGCNA) to all samples across the different cotton species. A total of 35 network modules were identified, containing from 102 to 4946 coexpressed member genes (Figure S5). By relating the representative coexpression pattern (i.e., module eigengene) of each module to the 29 phenotypic traits previously measured (Dong et al., 2020), we found that nearly all modules (except ME23 and ME32) were significantly correlated with at least one trait (Figure S5d). More than half of these modules showed significant correlations with either the categorical ‘species’ (i.e., samples from the species of interest were labeled 1; other sample were labeled 0) or ‘condition’ traits (i.e., the treated samples were labeled 1; other samples were labeled 0) (Figure 2a). Eight modules exhibiting a significant ‘condition’ correlation, including ME4, ME6, ME10, ME16, ME17, ME20, ME22, and ME30 (Figure 2b), were further examined to characterize the salt‐responsive coexpression dynamics and functional association.

Figure 2

WGCNA of four cotton species in response to salt stress. (a) Heatmap showing Pearson correlation between module eigengenes (by row) and categorical traits of ‘condition’ and ‘species’ (by column). The ‘condition’ trait is labeled 1 for the salt stress‐treated samples and 0 for the untreated samples. The ‘species_A2’ trait, for example, is labeled 1 for A2 and 0 for the other samples. Each cell shows the correlation coefficient with significance as follows: *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. (b) Eight modules of interest with significant correlations. Red circles indicate positive correlations, blue circles indicate negative correlations, and gray circles indicate non‐significant correlations. (c) Barplots of module eigengenes, with error bars denoting the standard deviation calculated from replicates.

According to the salt tolerance ranking (AD4 > A2 > D5 > AD1) reported previously, we were particularly interested in modules exhibiting expression profiles correlated with the ranking; however, no such modules were identified. Among the eight salt‐responsive modules, ME6 and ME17 were significantly correlated with the two more tolerant species G. arboreum and G. mustelinum (Figure 2b). Interestingly, both modules were positively correlated with G. mustelinum but negatively correlated with G. arboreum; that is, the highest and the lowest expression levels were observed in G. mustelinum and G. arboreum, respectively. No GO terms were significantly enriched for ME17, but the functional enrichment of ME6 pointed to biological regulation, oxidation–reduction, reproduction, and protein phosphorylation (Table S5). ME4 showed the lowest expression in G. arboreum, whereas ME16 showed the highest expression in G. arboreum (Figure 2c and Table S5), both exhibiting increased expression by salt treatment and enriched for ‘translation’. ME20 and ME30 showed opposite AD1 and A2 ‘species’ correlation patterns, with both exhibiting decreased expression by salt treatment; no GO terms were significantly enriched, probably due to their small module sizes (281 and 134 member genes, respectively).

Without significant ‘species’ correlations, ME10 exhibited reduced expression by salt treatment in all four cotton species (condition‐specific); this module was significantly enriched with GO terms related to photosynthesis (Figure 2b). These results are also consistent with the positive correlation between ME10 and photosynthetic‐related traits (Figure S5d). ME22 exhibited the opposite condition‐specific pattern to that of ME10, where the genes whose expression increased during salt treatment were strongly enriched for GO terms associated with ‘endoplasmic reticulum membrane’. Together, our module‐specific GO analysis further indicated that salt resistance is mediated by divergent strategies in the four closely related Gossypium species.

Gossypium mustelinum shows lower level of transgressive downregulation upon salt stress than does G. hirsutum

Despite their shared ancestry, the two allopolyploid cottons represent the most (G. mustelinum) and least (G. hirsutum) salt‐tolerant species surveyed here, while their model progenitors were intermediate to the two, prompting the question of how these two allotetraploid cottons respond to salt treatment relative to each other and to their parental G. arboreum and G. raimondii diploids. For each trio of allopolyploid + diploids (i.e., AD1 versus A2 versus D5 or AD4 versus A2 versus D5), DEGs with significant species:treatment interaction effects from our initial analyses were extracted to compare their expression changes by salt treatment (i.e., log2(fold change) values). Heatmap and clustering analysis revealed that gene regulation patterns in both polyploids were first clustered together with G. raimondii rather than G. arboreum (Figure 3), consistent with the higher number of shared pairwise DEGs between both allopolyploids and G. raimondii (Figure 1c,d).

Figure 3

Heatmap of expression changes by salt treatment in both allopolyploids relative to their parental parents. (a) The trio of AD1, A2, and D5. (b) The trio of AD4, A2, and D5. For 12 979 genes with significant species:treatment interaction effects, the log2 treated/control fold change values were normalized and scaled by z‐score in each trio. Beside each gene cluster, the cluster size and representative enriched functional classification terms are noted. Clusters displaying transgressive upregulation in allopolyploids are boxed in red, and clusters displaying transgressive downregulation are boxed in blue.

Additionally, both comparisons exhibited allopolyploid‐specific expression changes outside the range of the parental diploids (Figure 3), either upregulated (red box, termed as ‘transgressive upregulation’) or downregulated (blue box, termed as ‘transgressive downregulation’) in the allopolyploid. Interestingly, the number of genes that showed transgressive upregulation in G. mustelinum (3205 genes, clusters 4 and 5) was higher than that in G. hirsutum (2948 genes, clusters 1 and 2). In contrast, there were 2629 genes (clusters 3 and 7) transgressively downregulated in G. hirsutum, and only 1034 genes were transgressively downregulated in G. mustelinum (cluster 1). These results suggest that salt exposure provokes a less dramatic transcriptional response in G. mustelinum than in G. hirsutum and that the lower level of transgressive downregulation in G. mustelinum may account for its higher resistance to stress. GO enrichment results of the transgressively upregulated genes were generally consistent between G. hirsutum and G. mustelinum; i.e., GO terms linked to ‘protein modification’, ‘gene expression regulation’, and ‘macromolecule metabolism’ were significantly overrepresented in both. On the contrary, the transgressively downregulated genes in G. hirsutum were associated with categories such as photosynthesis, DNA replication, and amide metabolism, whereas transgressively downregulated genes in G. mustelinum were specifically involved in protein polymerization (full GO lists are provided in Table S6). In addition to the transgressive regulatory patterns, clusters displaying expression changes intermediate to those of diploids were also observed, with G. mustelinum (clusters 3 and 7, 3744 genes) containing more genes than those of G. hirsutum (cluster 5, 1992 genes).

Stress‐responsive homoeolog expression bias is novel in allopolyploids

In addition to the above analyses of total homoeolog expression in allopolyploids, we next asked whether salt stress treatment also affects relative homoeologous contributions to the total expression. In other words, we tested whether patterns of HEB were different between treated and control samples. PCA of subgenomic expression profiles (Figure S1b) showed that the subgenomic expression in allopolyploid homoeologs is, in general, closer to that of their corresponding diploids than to each other. This observation suggested that expression divergence between homoeologs in the allopolyploids mostly reflects divergence between parental diploids, or ‘parental legacy’ (Buggs et al., 2014). To further characterize the homoeolog expression divergence that is novel in each allopolyploid and/or in response to salt stress treatment, we applied a multifactor design (homoeolog + ploidy + treatment + homoeolog:ploidy + homoeolog:treatment + ploidy:treatment + homoeolog:ploidy:treatment) to each quartet of (sub)genome‐specific gene expression profiles (e.g., AD1:At versus AD1:Dt versus A2 versus D5). Four types of HEB patterns were categorized to dissect the static and conditional HEB for each allopolyploid relative to the parental diploids. In G. hirsutum for example (Figure 4), a total of 4738 homoeolog pairs exhibited type I HEB characterized by a significant homoeolog effect, which represent A versus D expression divergence at both the diploid (A2 versus D5) and allopolyploid (At versus Dt) levels regardless of salt treatment; in other words, this category of HEB was inherited from the diploid parental state and is constitutive. Type II HEB, characterized by a significant homoeolog:polyploid interaction effect, was detected for 4787 homoeolog pairs and represents novel HEB patterns not previously detected in diploids. In contrast to type I and type II patterns, which are not condition‐specific, type III and type IV HEB are responsive to salt treatment, as characterized by the significant homoeolog:treatment and homoeolog:treatment:polyploidy interaction effects, respectively. Like type I, type III effects represent those inherited from the parental state and type IV are salt‐responsive effects novel in allopolyploids. Whereas inherited and novel categories for both polyploid species are nearly equivalent under normal conditions (types I and II, Figure 4a), both polyploids exhibit a similar excess of novel HEB during salt stress (average 152 versus 1732 homoeologous pairs; Figure 4b), suggesting that the salt‐responsive HEB was predominantly novel upon allopolyploidization rather than inherited from diploids.

Figure 4

Homoeolog expression bias (HEB) categorization in G. hirsutum (AD1) and G. mustelinum (AD4). (a) Constitutive, not affected by salt stress. (b) Interaction, affected by salt treatment.

For genes exhibiting stress‐responsive HEB (types III and IV), we further examined the magnitude and direction of changes in bias between control and treated conditions (Table 2). Under each condition, pairwise At versus Dt comparisons revealed more biases under salt treatment than under control conditions in both allopolyploid species (G. hirsutum, 24.2% versus 10.5%; G. mustelinum, 18.3% versus 12.9%). In response to the salt stress treatment, 387 and 282 homoeologous pairs gained new HEB patterns in G. hirsutum and G. mustelinum, respectively, while an additional 125 and 179 pairs lost the HEB exhibited under control conditions; fewer pairs switched direction of bias. GO analyses revealed that genes with conserved patterns of HEB were enriched in photosynthesis‐related terms in both species (Table S7). Genes that gained bias in G. hirsutum under salt conditions were associated with catalytic activity and transcription factors while no statistically significant GO terms were enriched for those that gained bias in G. mustelinum, which might suggest diverse functions. Because no significant difference was detected between the numbers of A‐ and D‐biased genes, the overall HEB was considered balanced, as previously defined (Grover et al., 2012), in all scenarios (Table 2). For example, of the 191 homoeolog pairs exhibiting HEB in G. hirsutum under control conditions, 89 homoeologous pairs exhibited At bias while 102 exhibited Dt bias (Chi‐squared test, P > 0.05).

Table 2

Stress‐responsive homoeolog expression bias (HEB) in allopolyploid cottons

Bias pattern	Control	Salt	AD1	AD4
			Number (%)	Number (%)
Direction conserved	At = Dt	At = Dt	1259 (68.9%)	1397 (72.4%)
	At > Dt	At > Dt	28 (1.5%)	28 (1.5%)
	At < Dt	At < Dt	28 (1.5%)	35 (1.8%)
Bias gain by stress	At = Dt	At > Dt	176 (9.6%)	150 (7.8%)
		At < Dt	201 (11.0%)	132 (6.8%)
Direction switched	At > Dt	At < Dt	5 (0.3%)	3 (0.2%)
	At < Dt	At > Dt	5 (0.3%)	6 (0.3%)
Bias loss by stress	At > Dt	At = Dt	56 (3.1%)	96 (5.0%)
	At < Dt		69 (3.1%)	83 (4.3%)
Total At bias under control conditions			89 (4.9%)	127 (6.5%)
Total Dt bias under control conditions			102 (5.6%)	124 (6.4%)
Total At bias under salt conditions			209 (11.4%)	184 (9.5%)
Total Dt bias under salt conditions			234 (12.8%)	170 (8.8%)

At = A‐homoeolog; Dt = D‐homoeolog; At = Dt, equal expression of homoeologs; At > Dt, biased expression of the At homoeolog; At < Dt, biased expression of the Dt homoeolog; AD1 = Gossypium hirsutum; AD4 = Gossypium mustelinum.

Stress‐responsive HEB can alter expression in one or both homoeologs

Although all species studied generally exhibited balanced HEB, we consider the possibility that salt stress may disproportionately alter homoeolog usage from one or the other of the two divergent and coresident genomes (‘At’ and ‘Dt’) in allopolyploid cotton, a phenomenon known as ‘homoeolog induction bias’, previously proposed in hexaploid wheat (Powell et al., 2017). For homoeologous pairs exhibiting stress‐responsive HEB in G. hirsutum, 283 pairs exhibited significant expression alteration of only the At homoeolog and 167 pairs exhibited significant expression alteration of only Dt (Figure 5a, first four black sets in top barplot); 595 pairs exhibited significant changes of both At and Dt in the same direction (middle two red sets); and only five pairs exhibited changes in opposite directions (last two blue sets). A similar pattern was observed in G. mustelinum, where 245 and 162 pairs exhibited significant expression alteration of At and Dt, respectively; 648 pairs exhibited significant changes in both At and Dt in the same direction; and only six pairs exhibited changes in opposite directions (Figure 5b). These results indicate that the salt treatment was more likely to affect both homoeologs than only one homoeolog (42.9% versus 57.1% in G. hirsutum; 38.4% versus 61.6% in G. mustelinum); however, when only one homoeolog was affected, it was more frequently an expression change in the At homoeolog, possibly indicating that the A‐subgenome is more responsive to salt stress than the D‐subgenome in allopolyploid cottons. Notably, G. hirsutum showed more down‐ than upregulation of both homoeologs, whereas G. mustelinum showed more up‐ than downregulation of both homoeologs, congruent with the aforementioned results.

Figure 5

Identification of the homoeolog regulation patterns in response to salt stress in G. hirsutum (AD1) (a) and G. mustelinum (AD4) (b). For homoeologous pairs exhibiting stress‐responsive HEBs, further differential expression analyses were applied to identify regulation patterns after salt treatment. For example, this analysis identified 389 At homoeologs that were upregulated and 494 that were downregulated in response to salinity stress in AD1. Also shown are UpSet plots with counts of differentially expressed homoeologs in each pattern. Counts represent homoeologous genes where one (no intersection) or both (intersection of two circles) homoeologs were differentially expressed.

DISCUSSION

Gossypium species exhibit divergent global transcriptional responses to salt stress

Salt tolerance is a complex quantitative trait that involves many genes and pathways. Global analysis of the stress‐responsive gene expression profiles have been widely applied to understanding the molecular mechanisms and evolution of salt tolerance in plants (Edelist et al., 2009; Peng et al., 2014; Guo et al., 2015; Zhang et al., 2016; Lei et al., 2018; Yuan et al., 2019; Wang et al., 2021). In spite of many previous RNA‐seq studies in cotton, the transcriptional landscape of salt tolerance diversity between allopolyploid and diploid cottons has not been explored. Based on our previous understanding of the diverse morphological and physiological responses of 17 Gossypium species and genotypes (Dong et al., 2020), we selected four genotypes to examine their transcriptomic responses to relatively long‐term (14 days treatment) and mild salinity (50 mm NaCl) stress. The two tetraploid species G. hirsutum (AD1) and G. mustelinum (AD4) represent the least and the most salt‐tolerant allopolyploid genotypes, respectively, outside the range of two parental diploid species G. arboreum (A2) and G. raimondii (D5). Aiming to dissect the transcriptional changes underlying the diversity of salt resistance in Gossypium, we specifically asked how the tetraploid species acquired a wider range of tolerance than their parental diploids and why the two allopolyploid species exhibit such high differences in salt tolerance.

Our transcriptomic results demonstrated that salt stress has a profound global effect that alters gene expression for about one third of the entire transcriptome across four species. Phenotypically, we previously observed that nearly all of the phenotypic and physiological traits related to the salt response were altered after the same treatment (Dong et al., 2020). These observations together lend support to the notion that plant salt tolerance is a complex emergent phenotype with a broad genomic underpinning (Gugger et al., 2017). Across the four species, multivariate and pairwise analyses of salt‐responsive changes revealed diversified transcriptional programs and functional implications of salt acclimation (Figure 1), and subsequent coexpression network analysis further identified species‐specific coexpression profiles that distinguish the functional responsiveness between species and ploidy levels (Figure 2). For example, the relatively low number of DEGs found in G. arboreum suggests that this species requires less transcriptional adjustment responding to salt stress compared to the other cottons; this pattern was also evident in the coexpression network profiles (Figure 2c). One potential explanation is that G. arboreum presents a more ‘stress‐ready’ state that is pre‐adapted for potential salt stress, in a manner similar to many halophytes or extremophytes (Kazachkova et al., 2013). Unfortunately, no wild form of G. arboreum has ever been discovered, so nothing is known about the conditions under which it evolved.

Interestingly, we found a D‐dominant salt‐responsive pattern in both allopolyploids, evidenced by the high number of salt stress DEGs shared between both tetraploids and G. raimondii (Figures 1c,d), as well as by the observation that the salt‐induced expression changes in both tetraploids were more similar to those of G. raimondii than to those of G. arboreum (Figure 3). These results provide a transcriptional basis for the phenotypic observation that the allopolyploid traits were more like those of G. raimondii under salt stress; interestingly, the allopolyploid traits under control conditions were more like those of G. arboreum instead (Dong et al., 2020). Hence, we conclude that the D‐genome progenitor, modeled here by G. raimondii, acts as the ‘dominant’ parent of allopolyploid cotton with respect to salt responsiveness. This conditional asymmetry in transcriptional profiles between allopolyploids and their diploid progenitors has previously been reported in allopolyploid coffee (Coffea arabica) under two temperature conditions (Bardil et al., 2011; Combes et al., 2013).

Two allopolyploid species, G. hirsutum and G. mustelinum, exhibit the lowest and highest tolerance among the four surveyed species, respectively, outside the range of parental diploids. From the perspective of functional genomics, it is of particular interest to characterize the molecular basis of higher tolerance in G. mustelinum, which can be the key for improving salt tolerance in cultivated G. hirsutum. Coexpression network analysis revealed that the module ME6 specifically showed the highest level of upregulation by salt treatment in G. mustelinum, which links its salt tolerance phenotype to effective scavenging of ROS, since genes in this module were significantly enriched in the GO term ‘oxidation–reduction process’ (Table S5), as reported in some halophytes (Flowers & Colmer, 2008). This is in agreement with our previous physiological results that G. mustelinum showed highest peroxidase (POD) activities under mild salinity stress (Dong et al., 2020). Sodium can enter plant cells through non‐selective cation channels that are regulated by salinity‐induced signals such as ROS and Ca2+ (Deinlein et al., 2014; van Zelm et al., 2020). In addition, the enrichment of phosphorylation‐related processes in this module (Table S5) implied that protein kinases involved in calcium flux mediation might have been induced during sodium signaling in G. mustelinum. For instance, the calcineurin B‐like proteins (CBLs) bind calcium and cause protein phosphorylation through their interaction with CBL‐interacting protein kinases (CIPKs), and such CBL–CIPK pathways have been widely studied (reviewed in Manishankar et al., 2018). As the best characterized CBL–CIPK pathway, the salt overly sensitive (SOS) pathway has been proposed to mediate cellular signaling during salt stress (Ji et al., 2013). By phosphorylating the H+/cation antiporters, SOS complexes can transport sodium out of the cell, thus maintaining ion homoeostasis (Halfter et al., 2000; van Zelm et al., 2020). To some extent, this can be reflected by our previous physiological data that G. mustelinum showed a higher K+/Na+ ratio under salt treatment (Dong et al., 2020).

Balanced homoeolog expression bias in response to salt stress in allopolyploid cottons

In addition to providing insights into the genetic underpinnings of salt responsiveness among Gossypium species, our transcriptome analysis also has relevance with respect to long‐standing questions regarding the functional and adaptive significance of polyploidy with respect to stress (Van de Peer et al., 2017; Van de Peer et al., 2020). Central to the connection between polyploidy and stress is the increased genomic complexity following whole genome duplication, which has been hypothesized to enable the rapid evolution of stress adaptation, hence enabling the establishing of polyploids under stressful environments. Proposed mechanisms responsible for increased biotic or abiotic stress tolerance include heterosis, expression dosage due to increased gene copy number, and functional plasticity of duplicate genes in polyploid species (Chen, 2007; Liu & Adams, 2007; Dong & Adams, 2011; Madlung, 2013; Liu et al., 2015; Takahagi et al., 2018).

At the transcriptome level, HEB is a widely observed phenomenon in plant allopolyploids. Earlier empirical examples suggest that allopolyploids selectively express homoeologs relevant for an appropriate response to environmental stresses. For instance, the allotetraploid grass Brachypodium hybridum shows long‐term heat stress tolerance similar to its diploid progenitor Brachypodium stacei, unlike the other progenitor Brachypodium distachyon, and this physiological trait can be attributed to the preferential usage of homoeologs originating from the progenitor B. stacei to maintain cellular homeostasis under stress (Takahagi et al., 2018). Powell et al. (2017) demonstrated that allohexaploid wheat displayed HEB in response to pathogen infection and that such biases were mainly contributed by the B‐ and D‐homoeologs. In a more recent study of B. napus, more AT‐ than CT‐biases were found across control and three abiotic stress conditions (i.e., cold, drought, and heat) (Lee and Adams, 2020). Here, in tetraploid cottons, the salt‐responsive HEB displayed an equivalent level of differentially altered At and Dt homoeologs (Figure 5; 52–54% versus 46–48%, respectively). In addition, the A‐ and D‐biases were statistically balanced under both control and treatment conditions (Table 2). This is consistent with previous studies in other polyploids such as coffee (Combes et al., 2012; Combes et al., 2013) and wheat (Bhanbhro et al., 2020), where the contribution of different subgenomes appeared to be balanced during varied growing conditions. These results collectively imply that allopolyploids might respond to environmental stresses through differential or subgenome‐specific usage of homoeologs, but that these responses will vary across species and conditions. We note that our study only investigated the cotton leaf transcriptional response to salinity stress at the seedling stage and that there is a clear need for further research using additional tissues and including different developmental stages.

Disentangling parental legacy from evolutionary responses in allopolyploid cottons

The analytical approach used in the present study was expressly designed to elucidate and distinguish evolutionary changes in stress responsiveness that arose subsequent to polyploidization from those that pre‐existed in the diploid progenitors and thus were vertically inherited. This is important in that a leading idea about the adaptive relevance of polyploidy is that polyploids might be able to respond to changing or stressful conditions better than do their diploid relatives (Van de Peer et al., 2017; Van de Peer et al., 2020). To gain insight into this we jointly modeled stress responsiveness and ploidy level, with respect to a partitioned analysis of homoeologs. This multivariate HEB analysis identified more than 9000 homoeologous gene pairs that showed biased expression patterns in both polyploids, irrespective of stress treatment (Figure 4a). Notably, about half of these HEB patterns were inherited from the diploid parental state (type I reflecting ‘parental legacy’), whereas the other half were novel in allopolyploids (type II reflecting ‘regulatory novelty’). Perhaps the most important result, though, from the general standpoint of understanding the significance of polyploidy in evolution, is that among the stress‐responsive HEB patterns (Figure 4b), over 90% were novel in allopolyploids (type IV). This sharp contrast between control and stress‐responsive HEB patterns suggests that allopolyploid‐specific regulatory novelty likely plays a prominent role in salt stress responsiveness. We speculate that this reflects the numerous complex and novel interactions among both homoeologous regulatory factors and their target genes in allopolyploids, relative to their diploid progenitors, that arise immediately with polyploidization (Hu and Wendel, 2019). This regulatory robustness of gene regulatory networks conferred by polyploidy in adaptations to stresses has been testified in yeast through evolution experiments (Keane et al., 2014; Selmecki et al., 2015), as well as by recent computational simulation work (Yao et al., 2019). That is, deleterious mutations can often be buffered by the more complex polyploid gene networks or make the polyploids slightly less fit under stable environments, whereas the newly formed interactions or disturbances of the network structure may actually allow a wider range of evolvability of polyploids (reviewed in Carretero‐Paulet and Van de Peer, 2020). In this light, we further note that many homoeologs displayed novel expression patterns under stress in polyploid cottons. Thus, it is reasonable to assume that the involved cis‐ and trans‐regulatory elements act as interaction partners to coevolve the new phenotypic space that is specific to polyploids. It is worth mentioning that the diploids studied here are not the actual parents of polyploid cottons but are instead the best living models of diploid donors. Thus, the expression divergence between allopolyploids and their model diploids parents can also result from the independent evolutionary trajectories of diploid sister genomes and their corresponding subgenomes in allopolyploids. Nevertheless, this study still allows us to evaluate the role of homoeolog transcriptional regulation in salt responsiveness of allopolyploids. Further analyses carried out on artificial, synthetic allopolyploids systems would help shed light on diploids parental legacy versus polyploids regulatory novelty in responding to environmental stresses.

One additional result relevant to this polyploid‐specific stress responsiveness becomes evident when one considers total homoeolog expression changes. As shown in Figure 1b, 10 690 and 8806 genes in G. hirsutum and in G. mustelinum, respectively, exhibited significant transcriptional alteration in response to modest salt stress. By examining total expression changes in allopolyploids relative to the parental diploids, we found that transgressively up‐ and downregulatory changes accounted for about one third of these changes (Figure 3). This large stress‐responsive alteration in polyploids also likely reflects the massive combinatorial increase in regulatory complexity that is inherent in an allopolyploid relative to a diploid nucleus. Furthermore, these transgressive regulatory changes distinguished the two allopolyploid species studied, G. hirsutum and G. mustelinum, in terms of global level of transgressive downregulation. That is, fewer salt‐responsive genes were transgressively downregulated in the more salt‐tolerant G. mustelinum than in the less tolerant G. hirsutum, which in turn attenuated the salt responses. Because the two allopolyploid species originate from a single allopolyploidy event (Grover et al., 2015), this divergence may lead to new understanding of the post‐polyploidization evolution of diversity in adaptation to salinity stress.

In conclusion, our study revealed divergent responses to salt stress among four closely related diploid and polyploid cottons. We (i) show that the transcriptomic landscape associated with the salinity response reflects both parental legacy and regulatory novelty in polyploids and (ii) disentangle these effects through statistical categorization of four distinct HEB patterns. Our findings highlight the transcriptional regulatory consequences of duplicated genes in response to salt stress of polyploids and provide new insights into the putative adaptive role of allopolyploidization during evolution under stressful environments.

EXPERIMENTAL PROCEDURES

Plant materials and salinity stress treatment

Four cotton species were selected based on our previous study (Dong et al., 2020), including two natural allopolyploid cottons, G. hirsutum (AD1) cv. TM1 and G. mustelinum (AD4), and their model diploid progenitors G. arboreum (A2) cv. 101 and G. raimondii (D5). Mature seeds were delinted with sulfuric acid and germinated at 30°C on moistened filter paper placed in Petri dishes. Germinated seeds were next transferred to a hydroponic system containing half‐strength Murashige and Skoog (MS) medium and grown in the Pohl Conservatory greenhouse at Iowa State University (Ames, IA). Uniformly developed 2‐week‐old seedlings were collected for each species and randomly divided into control and treatment groups, each with three biological replicates (three to four plants per replicate). Salt treatment of 50 mm NaCl for 2 weeks was applied as previously described (Dong et al., 2020).

RNA extraction, sequencing, and data processing

Fully expanded true leaves (seventh to ninth) were collected and pooled from three to four individual plants per replicate. Harvested tissues were flash‐frozen in liquid nitrogen and stored at −80°C until extraction. Total RNA was isolated using the Qiagen Plant RNeasy Kit (Qiagen, Stanford, CA, USA) following the manufacturer's instructions. The concentration and integrity of extracted RNA were determined using an Agilent Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). RNA samples were sent to the Beijing Genomics Institute (BGI, Hong Kong) for library preparation and 150‐bp paired‐end sequencing on a BGISEQ‐500 platform. Raw reads were deposited at the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) under project number PRJNA601953.

RNA‐seq reads were mapped against the G. raimondii reference genome (D5; 37 505 protein‐coding genes) (Paterson et al., 2012) using the single nucleotide polymorphism (SNP)‐tolerant mode of the Genomic Short‐Read Nucleotide Alignment Program (GSNAP) (Wu and Nacu, 2010). A set of pre‐built, species‐specific SNP indices (Page and Udall, 2015) were used to distinguish A‐ and D‐genome orthologs and homoeologs. For allopolyploids, mapped reads were partitioned by polyCat (Page et al., 2013; Page et al., 2014) as homoeolog A‐specific, D‐specific, unclassified, or chimeric reads; because the latter two lack homoeologous SNPs or contain both A‐ and D‐SNPs, respectively, reads in these categories were excluded from subsequent analyses as per (Hu et al., 2020). Read counts were generated using featureCounts from the subread package v.1.6.0 (Liao et al., 2014). The total expression of homoeologous genes in the allopolyploids was calculated as the sum of the homoeolog A‐ and D‐specific expressions.

Differential gene expression analysis

Raw read counts were used for differential expression analysis with the DESeq2 package (Love et al., 2014) in R v. 3.6.3 (R Core Team, 2017). First, a multifactor design (species + treatment + species:treatment) was used to differentiate gene expression changes due to species, salt treatment, or the interaction between the two, using a likelihood ratio test‐based approach to compare the full model against the corresponding reduced models. Next, pairwise comparisons were conducted with Wald tests to determine DEG sets between control and salt‐treated conditions, between any two genotypes, and due to species‐specific responses to salinity (interaction). The resulting P‐values were corrected using a false discovery rate (FDR) of α = 0.05 by the Benjamini–Hochberg method (Benjamini and Hochberg, 1995). The numbers of DEGs were divided by the total number of genes in the reference genome (i.e., 37 505) to obtain DEG percentages for direct comparisons (Rapp et al., 2009; Bardil et al., 2011). Clustering analysis of DEGs with interaction effect was performed by the Complexheatmap package (Gu et al., 2016) with the Spearman distance matrix of z‐score scaled log2(fold change) values.

Homoeolog expression bias

HEB refers to differential expression between homoeologous genes in allopolyploid species (Grover et al., 2012); in tetraploid cotton, an A‐bias indicates higher expression of the At than the Dt homoeolog, and a D‐bias indicates higher Dt than At expression. In addition to this classic definition of HEB under static conditions (e.g., in either treated or control plants), we are also interested in the conditional HEB in response to stress treatment. Therefore, we applied a multifactor design to each allopolyploid species relative to the parental diploids: expression homoeolog + ploidy + treatment + homoeolog:ploidy + homoeolog:treatment + ploidy:treatment + homoeolog:ploidy:treatment, where (sub)genome‐specific gene expression is subject to the effects of genome origin (homoeolog: A2/At versus D5/Dt), ploidy level (ploidy: allopolyploid versus diploid), and experimental condition (treatment: control versus treated). Based on the significant effects of homoeolog‐related terms, four HEB patterns were inferred as follows: homoeolog, type I HEB that is inherited from the parental state and not affected by treatment; homoeolog:ploidy, type II HEB that is novel in allopolyploids and not affected by treatment; homoeolog:treatment, type III stress‐responsive HEB that is inherited from the parental state; and homoeolog:ploidy:treatment, type IV stress‐responsive HEB that is novel in allopolyploids. An FDR‐adjusted P‐value of 5% was used to identify significant effects.

Functional enrichment analysis

GO enrichment analysis was performed using the R package topGO (Alexa and Rahnenführer, 2016). Briefly, Fisher's exact test was applied to the gene sets of interest to detect the overrepresented GO terms against the universal set of genes. The P‐values were adjusted for multiple testing using the Benjamini–Hochberg method (Benjamini and Hochberg, 1995) to retain significant GO terms with an adjusted P < 0.05. GO annotations of the G. raimondii reference genes were downloaded from CottonGen (Paterson et al., 2012; Yu et al., 2014). GO enrichment results were summarized using the REVIGO web server (Supek et al., 2011).

Weighted gene coexpression network analysis

Coexpression network analysis was conducted using the R package WGCNA (Langfelder and Horvath, 2008). The raw read count data were normalized using the regularized log transformation method built in DESeq2 (Love et al., 2014). After removing genes with zero expression or without variation across samples, coexpression modules were obtained using the automatic network construction function blockwiseModules. Briefly, a similarity matrix was calculated based on Pearson correlation indices between each pair of genes, which was then applied to generate a ‘signed’ adjacency matrix with an optimized soft threshold of β = 26. Network interconnectedness was next measured by calculating the ‘signed’ topological overlaps using the TOMdist function. By performing average linkage hierarchical clustering on the topological overlap dissimilarity matrix (DistTOM = 1 − TOM), genes that were highly connected were grouped together into coexpression modules with a dynamic tree cutting algorithm. The parameter mergeCutHeight = 0.25 was used to merge highly similar modules.

For each coexpression module, a module eigengene that characterizes the representative gene expression pattern was calculated as the first principal component of the scaled module member gene expression profiles. Using the 29 phenotypic and physiological traits data obtained from our previous study (Dong et al., 2020) (the same accessions as those used in the present study), module–trait relationships were examined by Pearson correlation indices between module eigengenes and trait measurements using the cor function in WGCNA. The resulting P‐values were corrected for multiple testing by the Benjamini–Hochberg method (Benjamini and Hochberg, 1995).

AUTHOR CONTRIBUTIONS

GH, YD, and JFW conceived the research; YD and ERM carried out the experiments; YD and CEG performed data analysis; YD wrote the manuscript; GH, JFW, SZ, and CEG contributed substantially to revisions.

CONFLICT OF INTERESTS

The authors have declared that no competing interests exist.

OPEN RESEARCH BADGES

This article has earned Open Data and Open Materials badges. Data and materials are available at https://www-ncbi‐nlm‐nih‐gov.ezproxy.u‐pec.fr/bioproject/601953 and https://github.com/Wendellab/SaltStressTranscriptome.

Figure S1. (a) Correlation heatmap of RNA‐seq samples. CK and Salt refer to control and salt treatment samples, respectively; numbers 1, 2, and 3 represent three biological replicates. (b) Principal component analysis of A‐ and D‐subgenomic expression in allopolyploid cottons (G. hirsutum AD1 and G. mustelinum AD4) and their diploid progenitors G. arboreum (A2) and G. raimondii (D5). Diploid samples are circled in black, and the one mislabeled D5 sample is circled in red.

Click here for additional data file.

Figure S2. GO enrichment of differentially expressed genes with significant species, treatment, and species:treatment effects. Enriched GO terms of biological processes were summarized and visualized by REVIGO, which uses multidimensional scaling (MDS) to reduce the dimensionality of a matrix of the GO terms pairwise semantic similarities. Circle color indicates the log10 adjusted P‐value and circle size indicates the GO term frequency.

Click here for additional data file.

Figure S3. GO enrichment of up‐ and downregulated genes in response to salt stress of each species.

Click here for additional data file.

Figure S4. GO enrichment of species‐specific up‐ and downregulated genes in response to salt stress, corresponding to the first four columns in Figure 1c,d. Results are shown for A2, D5, and AD4 genes, because no significant GO terms were enriched in AD1‐specific differentially expressed genes.

Click here for additional data file.

Figure S5. WGCNA network construction. (a) Selection of the soft‐thresholding powers. The scale‐free fit index and the mean connectivity on the y‐axis are plotted against soft‐thresholding power on the x‐axis. (b) Coexpression modules were identified by cutting the hierarchical clustering tree. (c) Eigengene dendrogram and correlation heatmap showing the relationships among identified modules. (d) Correlation heatmap between module eigengenes and traits. In each cell, color and text indicate correlation and P‐value.

Click here for additional data file.

Figure S6. GO enrichment of module member genes.

Click here for additional data file.

Table S1. Summary of RNA‐seq data generated in this study.

Click here for additional data file.

Table S2. GO enrichment analysis of differentially expressed genes exhibiting species, treatment, and interaction effects.

Click here for additional data file.

Table S3. GO enrichment analysis of differentially expressed genes in response to salt stress in four cotton species.

Click here for additional data file.

Table S4. GO enrichment analysis of species‐specific salt stress‐responsive genes.

Click here for additional data file.

Table S5. GO enrichment analysis of WGCNA coexpression modules.

Click here for additional data file.

Table S6. GO enrichment analysis of genes from each cluster.

Click here for additional data file.

Table S7. GO enrichment analysis of genes showing responsiveness HEB patterns in polyploids.

Click here for additional data file.