Comparative transcriptome analysis reveals evolutionary divergence and shared network of cold and salt stress response in diploid D-genome cotton

Abstract

Background: Wild species of cotton are excellent resistance to abiotic stress. Diploid D-genome cotton showed abundant phenotypic diversity and was the putative donor species of allotetraploid cotton which produce the largest textile natural fiber.

Results: A total of 41,053 genes were expressed in all samples by mapping RNA-seq Illumina reads of G. thurberi (D₁), G. klotzschianum (D_3-k), G. raimondii (D₅) and G. trilobum (D₈) to reference genome. The numbers of differently expressed genes (DEGs) were significantly higher under cold stress than salt stress. However, 34.1% DEGs under salt stress were overlapped with cold stress in four species. Notably, a potential shared network (cold and salt response, including 16 genes) was mined out by gene co-expression analysis. A total of 47,180-55,548 unique genes were identified in four diploid species by De novo assembly. Furthermore, 163, 344, 330, and 161 positively selected genes (PSGs) were detected in thurberi, G. klotzschianum, G. raimondii and G. trilobum by evolutionary analysis, respectively, and 9.5-17% PSGs of four species were DEGs in corresponding species under cold or salt stress. What's more, most of PSGs were enriched GO term related to response to stimulation. G. klotzschianum showed the best tolerance under both cold and salt stress. Interestingly, we found that a RALF-like protein coding gene not only is PSGs of G. klotzschianum, but also belongs to the potential shared network.

Conclusion: Our study provided new evidence that gene expression variations of evolution by natural selection were essential drivers of the morphological variations related to environmental adaptation during evolution. Additionally, there exist shared regulated networks under cold and salt stress, such as Ca²⁺ signal transduction and oxidation-reduction mechanisms. Our work establishes a transcriptomic selection mechanism for altering gene expression of the four diploid D-genome cotton and provides available gene resource underlying multi-abiotic resistant cotton breeding strategy.

Keywords: Co-expression; Comparative transcriptome; Diploid D-genome cotton; Evolutionary divergence; Shared network.

Fig. 1

Phenotypic variations and correlation of whole-genome expression of four species. a Phenotypic variations in flower color and leaf shape. b Tolerance divergence of four species in seedlings. c Heatmap of correlation value (R square) of 40 libraries

Fig. 2

Sequence diversity of four species of cotton. a Unrooted phylogenetic tree of four species using SNPs, which obtained from transcriptomic data. The scale bar represents the expected number of substitutions per site. b The relationship showed by the principal component cluster among four species. c Phylogenetic tree of nine genomes using the identified orthologous genes. The scale bar represents the expected number of substitutions per site. d Boxplot of the dN/dS ratio of nine genomes. Wild species, red boxes. Cultivated species, green boxes

Fig. 3

Gene expression pattern across three-time points (C0, C6, and C12) of four species under normal condition and corresponding top fifteen enrichment GO terms. a Gene expression pattern across three-time points. Eight gene clusters (profile 1–8) were identified using k-means clustering. In each cluster, the y-axis represents log₂ (FPKM+ 1) derived from RNA-seq data for each biological sample, while the x-axis represents the biological samples that are ordered as C0 (R1 and R2), C6 (R1 and R2), and C12 (R1 and R2) for each species. b Heatmap of –log₁₀ (p-value) of biological process category enrichment among the eight profiles

Fig. 4

Co-expression network analyses by WGCNA. a Hierarchical cluster tree showing co-expression modules identified by WGCNA. Each leaf in the tree represents one gene. The major tree branches constitute 33 modules labeled with different colors. b Module–sample association. Each row corresponds to a module labeled with color as in (A) Modules are distinguished by different colors which were arbitrarily assigned by the WGCNA package. Each column corresponds to a tissue type as indicated. The color of each cell at the row-column intersection indicates the correlation coefficient (R) between the module and the tissue type. *Significance at P < 0.05; **Significance at P < 0.01

Fig. 5

Co-expression network of hub gene in skyblue3 modules. Genes were related to responses to stimuli in a red background. Five function unknown genes were displayed in white background

Fig. 6

Characterization of DEGs under cold and salt stress. a Numbers of DEGs of four species under cold and salt. Red block means numbers of down-regulated genes. Blue block means numbers of up-regulated genes. b DEGs of GD3 under cold and salt stress; PSGs of GD3; and genes of skyblue3 module showed by Venn. c Common DEGs of GD3 under cold and salt stress

Fig. 7

Phenotype observed in the silenced plants with the TRV: 00 empty vectors, wild type plants, and the silenced plants at 12 days post-inoculation. a Albino’s appearance on the leaves of the PDS infused plants. b RT-qPCR analysis of the change in the expression level of the RALF and FLA genes in cotton plants treated with VIGS. c Evaluation of fresh shoot biomass, fresh root biomass, Na^+, and K⁺ ions concentration and their ratios. Letters a/b indicate statistically significant differences (two-tailed, p < 0.05). The error bars of the RALF and FLA gene expression level represent the standard deviation of three biological replicates