Additionally, only three (Gohir.A02G026700,Gohir.D02G035100, and Gohir.D13G131100) of the 21 structural genes were correlated to specific regulatory genes, but the other 86% (18) shared at least one regulatory gene with other genes (Table S5). For example, a structural gene of GSTU(Gohir.D02G035100) was only correlated with a regulatory gene ofNAC017 (Gohir.D02G067300), whereas the other structural gene of GSTU (Gohir.D13G158600) was simultaneously correlated with three regulatory genes, NAC017(Gohir.D02G067300), C2H2 (Gohir.A11G051900), andWRKY28 (Gohir.D10G039400) (Fig. S2, Table S5).
Among the 98 regulatory genes, we identified a total of 26 homoeologous pairs (52 genes) belonging to 13 gene families. Although the At- or Dt-homoeolog in each pair encode the same enzyme, they were not always co-expressed with the same set of target genes (Fig. S2, Table S5). In some cases, At- or Dt-homoeolog was correlated with different target genes in the same or different co-expression networks and there was no subgenomic selection preference, i.e., At- or Dt-homoeolog don’t only select genes with the same A or D subgenome for highly correlation, respectively. For example, the At-homoeolog of a WRKY24 gene (Gohir.A03G025100) was highly correlated with three structural genes, including a GSTU (Gohir.D02G034300), aPGDC3 (Gohir.A08G052200), and a GSTL3(Gohir.A09G160500) genes. However, the Dt-homoeolog of the WRKY24(Gohir.D03G143400) was correlated with only a GSTU gene (Gohir.D02G034300). In addition, c o-expression network analysis further showed that all but one gene (Gohir.D07G048100 ) of the 98 regulators positively regulated to their corresponding target structural genes (Table S5). These results indicate that the evolution of regulatory system for the duplicated GSH metabolic pathway in polyploid is complicated, reticulately correlated, with no independent parallel evolution of subgenomic components.
3.4 Expression characteristics of the structural genes in the GSH metabolism pathway under UV-B stress
In order to further elucidate the composition and express patterns of structural genes in the GSH metabolic pathway under UV-B stress, we performed a RNA-sequencing on cotton leaves under both UV-B stress and control conditions. Six high-quality transcriptomes (with mapping rates ≥ 94%, GC contents ≥ 43%, and Q30 values ≥ 93%), including approximately 41–45 million raw reads from each transcriptomes, were obtained (PRJNA893188). The clean reads ratio of the transcriptomes ranges from 98% to 99% (Table S6).
Of the 205 structural genes characterized at the genome-wide level, only 98 were DEGs under UV-B stress (Fig. 3, Table S3). In addition, these DEGs occurred in most structural gene families (G6PDH ,GGCT , GGT , GPX , GR , GST , IDH ,LAP , OXP , and PGDC ) in the GSH metabolic pathway except for the GCL and GS families. The expression direction of these DEGs was not completely consistent, with 68 up-regulated and 30 down-regulated. Within gene families, two expression patterns have been identified. One pattern was that all gene members of the same family showed similar expression changes under UV-B stress. For example, gene members in the GGCT , GR , LAP ,OXP , and PGDC families were totally up-regulated, whereas members in the GGT family were completely down-regulated. Another pattern was that the expression of different gene members of the same family (G6PDH , GPX , GST, and IDH ) was alternatively regulated in an upward or downward manner, respectively.
Furthermore, we found that only 19 of these 98 DEGs fell into the list of 21 structural genes in the co-expression networks (Table S7). AGSTU gene, Gohir.D02G034300 , was correlated with the largest number of regulatory genes (53). In contrast, two otherGSTU genes, Gohir.A02G026700 and Gohir.D02G035100 , were correlated with a single regulatory gene each.
Among the 205 structural genes, 80 homoeologous pairs (160 genes) were identified (Table S3), of which, 72 pairs of homoeologs (144 genes) were simultaneously expressed under both CK and UV-B stress conditions. The relative expression pattern of At- to Dt-homoeolog of each homoeologous pair could be classified into unbiased (At/Dt=1) and biased (At/Dt≠1) expression. Unbiased expression included no bias (33 pairs, 66 genes) and lost bias (14 pairs, 28 genes); and biased expression included four patterns of At bias (six pairs, 12 genes), Dt bias (nine pairs, 18 genes), get At bias (six pairs, 12 genes), and get Dt bias (four pairs, eight genes). In addition, in the no biased pattern (At/Dt=1 under both CK and UV-B stress), only 13 homoeologous pairs showed significant changes in the expression of both At- and Dt- homoeologs under UV-B stress (Table S8). Thus, the effects of UV-B stimulation on At- and Dt- homoeologs at the same structural gene loci in the GSH metabolic pathway tend to more diverse and independently driven.
3.5 Potential regulatory genes involved in the GSH metabolic pathway under UV-B stress
Through weighted co-expression network analysis, we found that only 21 structural DEGs in the GSH metabolic pathway could establish a high degree of co-correlation with regulatory genes, based on which we therefore uncovered 98 regulatory genes that may be involved in the GSH metabolic process. Comparative transcriptome analysis showed that only 96 of these genes were differentially expressed under UV-B stress, of which 95 were up-regulated and one was down-regulated (Table S7). Therefore, these 96 genes were used for the subsequent analyses (Table S7). These regulatory genes can be classified into 13 gene families, theWRKY family had the most abundant gene members (24), and theB3 , HD-Zip , NF-X1 , and RAV families had the least gene members (2). In addition, from the co-expression networks, aWRKY28 gene, Gohir.D10G039400 , was correlated with the largest number of structural genes (6), but 18 regulatory genes frombHLH , C3H , ERF , GRAS , HSF NAC ,WRKY gene families were each co-expressed with a corresponding structural gene (Table S7).
Twenty-five homoeologous pairs (50 genes) in the 96 regulatory DEGs were identified (Table S8). The homoeologous pairs belonging to theB3 , bHLH , C2H2 , C3H , ERF ,GRAS , HSF , MYB , NAC , NF-X1 ,RAV , and WRKY gene families exhibited consistent up-regulation of their expressions under UV-B stress. Furthermore, the patterns of relative expression between At- and Dt-homoeologs of these regulatory genes could be categorized into the unbiased (20 no bias, and two lost bias) and biased expression (one At bias, one get At bias, and one get Dt bias) (Table S8). These results show that, in most cases, the effects of UV-B stress on At- and Dt homoeologs at each regulatory locus are often unbiased.