Oxidant and antioxidant evaluation of the three upland cotton accessions CRI12, LAT40 and MAR85 under salt stress conditions
The three upland cotton genotypes morphologically are identical and exhibit no significant variations when exposed to any form of stress, thus no consideration of various morphological traits were considered in this study. However, we observed variation on root growth and biomass accumulation, CRI12 and LAT40 had a relatively higher root and overall higher biomass compared to MAR85 (Figure 1A). Moreover, in carrying out deep analysis on the concentration levels of oxidant and antioxidant enzymes on the leaf tissues of the three cultivars under salt stress conditions, MAR85 and CRI12 showed significantly higher concentrations of proline and superoxide dismutase (SOD) compared to LAT40 (Figure 1B), an indication that CRI12 and MAR85 were less affected under salt stress compared to LAT40. Furthermore, Malondialdehyde (MDA) concentrations levels within the plants is an indication that the plants are suffering from oxidative stress, being MDA is a byproduct of lipid peroxidation . The results obtained were in agreement to previous findings in which the knockdown of trihelix transcription factor in cotton reduced drought and salt stress tolerance and in turn accelerated the accumulation of various oxidants and MDA levels under drought and salt stress exposure .
Sequencing and Transcript Identification
The transcription analyses of three accessions of G. hirsutum differed significantly in their morphological characteristics as well as their crucial systems-level insights into molecular mechanisms that underlie the alkali-salt stress response [32,33]. In this study, one upland cotton cultivar (CRI12, China) and two other wild upland cotton cultivars (LAT40, G. hirsutum race latifolium40 and MAR85, G. hirsutum race mari-garant85) were used for transcription analysis by carrying out RNA-Sequencing on the their tissues when exposed to salinity stress. Two organs were considered, the leaf and root of the three accessions, at four different treatment stages; the samples were coded as Rt0h, Rt3h, Rt12h, Rt48h for the root samples and leaf samples were coded as Lf0h, Lf3h, Lf12h and Lf48h, all the samples were replicated three times. The performance of the three upland cotton accessions was categorized into four distinct groups, as normal growth stage (CK0), early alkali-salt stress response stage (SS3), seedlings significantly damaged stage (SS12) and seedling recovered stage (SS48). In order to investigate the transcriptome dynamics throughout alkali-salt stress treatment, total RNA was isolated from the leaf and root of three accessions of G. hirsutum in four different treatment stages. In total, 1,106,485,712 numbers of the raw reads were produced from 96 cDNA libraries, after cleaning, 1,097,617,134 (99.19%) clean reads were obtained. The percentage of mapped reads among clean reads in each library ranged from 83.88 to 89.61%, moreover the percentage of clean reads with a Phred quality value of 30% ranged from 94.37% to 97.05%. Furthermore, the clean reads were aligned to the reference genomeof Gossypium hirsutum (AD1), in which 83.06% to 89.73% reads generated from the 96 samples were mapped to reference genome, producing 77.4%—82.3% of the uniquely mapped reads to the reference genome. hirsutum genome. What was more interesting; the uniquely mapped reads were mapped to the reference genome by use of HT-seq (Python package), with higher efficiency and accuracy (Table 1)
Transcription and Expression Analysis in Various Tissues of the Three Cotton Species under Salt Stress Condition
A total of 64,737 genes were identified, and the number of expressed genes in different samples varied from 51,586 to 57,263 as detected by RNA-sequence data analysis (Figure 2A). To comprehend the global differences in the transcriptome dynamics during dissimilar treatment stages, the expression distribution analysis, correlation analysis, principal component analysis (PCA) and hierarchical clustering have realistically accomplished based on FPKM values for all the expressed genes in at least one of the 32 tissue samples (Figure 2B). The leaf samples showed a lower expression level than root, and same treatment tissue/treatment stage of three accessions of G. hirsutum shown the similar expression distribution. The Pearson correlation coefficient analysis, among the combined 32 samples illustrated more significant correlation in the same tissue/treatment stage among three accessions of G. hirsutum (Figure 2C). As expected, leaf transcriptome of three contrasting upland cotton races clustered together and showed substantial differences with treatment points. The LAT40 and CRI12 showed a more exact correlation and similar phenotype in the leaf after alkali-salt treatment indicated high similarity in their transcriptional programs. Compared with CK0 stages, SS3, SS12 and SS48 demonstrated closer in the root and leaf of three accessions of upland cotton (Figure 2D). It suggested a significant difference of transcriptional programs between the standard and alkali-salt stress condition. Overall, tissues/treatment stages exhibited a higher correlation in these analyses that expected to have more similar transcriptomes and functions/activities.
Differential gene expression after alkali-salt treatment
To investigate the differential expression pattern of the various transcriptional factors in the three upland cotton accessions, two tissues were profiled, the leaf and root tissues under salt stress condition. The total number of the differentially expressed genes (DEGs) varied, ranging from 8663 (MAR85L_CK0 vs. MAR85L_SS12) to 22068 (CRI12R_CK0 vs. CRI12R_SS12) furthermore, the DEGs were obtained from pairwise comparison (MAR85L: MAR85L_CK0 vs. MAR85L_SS3; MAR85L_CK0 vs. MAR85L_SS12; MAR85L_CK0 vs. MAR85L_SS48; MAR85R, LAT40L/R and CRI12L/R pairwise comparison was similar with MAR85L) by using DEGseq software (Figure 3A). A greater number of DEGs observed in roots (a total of 41,132 DEGs) than it in leaves (a total of 35,724 DEGs) which indicated that the roots could be the main tissue affected by salt stress and thus have a more dynamic and complex gene regulation to reduce the toxicity of salts to the root cells. The results obtained are in agreement to previous findings in gene expression pattern tends to be tissue specific Magwanga et al , found that high number of the LEA genes were highly upregulated in the leaf tissues compared to stem and root tissues under drought stress condition. Compared with SS3 (21,738 and 30,525, in leaves and roots, respectively) and SS48 (23,418 and 26,533, in leaves and roots, respectively) stages, the SS12 stage (28,521 and 32,560, in leaves and roots, respectively) exhibited the highest number of DEGs in root and leaf, indicating SS12 stage was more activated for response to alkali-salt stress. Interestingly, MAR85 showed the lowest number of the DEGs in all the stages compared to LAT40 (highest number of DEGs) and CRI12. It indicated that the inconsistent responses to alkali-salt stress of three accessions of upland cotton. Additionally, 3,509 (SS3 VS. CK0) 8,138 (SS12 VS. CK0), and 5,955 (SS48 VS. CK0) common DEGs were identified at each stage subsequently alkali-salt treatment in the leaves of three accessions of upland cotton. Then, 8,870 (SS3 VS. CK0), (Figure 3E), 10,428 (SS12 VS. CK0), and 7,281 (SS48 VS. CK0), DEGs were identified in the roots (Figure 3B). The variable number of differentially expressed genes suggested that each tissue/alkali-salt treatment stage retained their own independent developmental programs. Alternatively, the transcriptional complexity may merely reflect the intricacy of the captured seed stages, which contained more than one cell type.
Functional annotation of DEGs sets in the three accessions of upland cotton at different tissues (leaf and root)
A total of 22,359 DEGs (11,818 and 15,674 in leaf and root, respectively) was used for functional enrichment and weighted co-expression network analysis (WGCNA). The cluster analysis, consequence of the 27,492 DEGs was similar to the cluster analysis of all the genes. This empirical observation supported those 22,359 DEGs as obtained from the leaf and root tissues in the three different cotton accessions could precisely illustrate the variations in the samples analyzed. The common DEGs of each tissue/stage in three accessions of upland cotton were further used for carrying out the gene ontology (GO) enrichment and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis. A total of 17,985 genes were annotated, with 7,836 and 10,149 being found to be differentially expressed in the leaf and roots, respectively. Among the 7,836 annotated DEGs, 105 GO terms were significantly enriched at P-value≤ 0.05, FDR≤ 0.05 (Table S2). In roots, more of the DEGs were annotated, however, only 96 GO terms were found to be significantly enriched (Table S3). Moreover, 54 GO terms were commonly enriched between the leaf and root tissues, including oxidation-reduction process (GO: 0055114), carbohydrate metabolic process (GO: 0005975), oxidoreductase activity (GO: 0016491), protein serine/threonine kinase activity (GO: 0004713) and transcription factor activity (GO: 0003700), which were known to be involved abiotic stress response.
Tissue-specific enrichment GO terms were observed where the DEGs obtained from the leaf were significantly enriched in photosynthesis (GO: 0015979) and photosynthesis-related GO terms (GO: 0019684, GO: 0009765). Moreover, all the 4163 DEGs (2426 and 2667 DEGs in leaves and roots, respectively) were enriched by using the KEGG database. The DEGs of the leaf and root were significantly enriched in 30 and 34 KEGG terms, respectively, with the Q value of ≤ 0.0001, and gene number for each analyzed term set at ≥ 3 (Table S4/S5). Plant hormone signal transduction (ko04075) and the biosynthesis of secondary metabolites (ko01110) were the most common KEGG pathways detected (Figure S2)., and the two pathways have been previously found to be significantly involved in abiotic stress response [35–38]. Moreover, the leaf-specifically enriched KEGG terms detected were those related to the photosynthetic apparatus of the plant leaf, such as Photosynthesis (ko00195) and Photosynthesis - antenna proteins (ko00196). Even though alkali-salt stress response do involved in a more complex pathways, signal transduction and oxidation-reduction processes, but the DEGS detected could be playing a core role in the leaf and root tissues of cotton under salt stress conditions. Apart from the DEGs, we analyzed the various plant transcription factors with putative function in enhancing salt stress tolerance among the three cotton cultivars used. A total of 2,080 TFs were identified, in which 991 TFs were detected among the DEGs in the leaf tissues while1, 577 TFs were from the DEGs analyzed from the root tissues roots (Table S6). The results obtained for the TFs correlated positively to the distribution of the DEGs in which more were found to be expressed in the roots as opposed to the leaf tissues, an indication that the roots could be the primary plant tissue, which is highly affected under salt stress conditions. In all the TFs identified, 53 different gene families were found to be linked to the TFs such as ZF-HD, SRS, S1Fa-like, SAP and NF-X1 TFs families were identified for the DEGs for the root tissues while YABBY TF family, was detected for the DEGs found for the leaf tissues.
Co-expression network and module construction
To investigate the gene regulatory network of alkali-salt stress response in upland cotton, we identified co-expressed gene sets via weighted gene co-expression network analysis (WGCNA) which can be used to find networks (modules) of highly correlated genes (Langfelder and Horvath, 2008). In the present study, 22,359 DEGs were used for WGCNA. In order to ensure a scale-free network, the power of β = 12 (scale free R2 = 0.8740) was accurately selected by determination of soft-thresholding power (Figure S3A/B) and evaluation of scale free topology analysis (Figure S3C/D). Representing interaction among genes with similar expression profiles, several networks that were referred to as co-expression modules, were identified. A total of 20 modules was identified via the Dynamic Tree Cut method (core parameter: MEDissThres = 0.25), varying from 135 genes in the mistyrose module to 4793 genes in the darkmagenta module. Among 20 modules, 17 genes were not useful, which displayed in the grey module, could not be selected by any other modules (Figure 4A-B). A thousand genes were randomly selected for visualization of gene networks. Genes in the same module shown a higher topological overlap (Figure S4), indicating networks (modules) which were constructed by WGCNA have biological useful functions.
Identification and visualization of tissue/stage-specific modules in leaf and root of upland cotton
A summary profile (eigengene) for each module obtained by WGCNA. For example, clustering heat map and eigengene bar plot of modules that highly correlated with traits, which were shown in (Figure 5A), it represented the gene expression levels of each module. In order to determine the correlation between module and trait, we associated eigengenes with alkali-salt treatment stages and upland cotton races via Pearson correlation coefficient analysis, respectively. Interestingly, none of the module significantly associated with three upland cotton races (r ≥ 0.60, p ≤ 0.01);; however, 11 co-expression modules showed a relatively higher correlation (r ≥ 0.60) with SS3, SS12, SS48 stages in leaf and root of upland cotton (Figure 5B-C). It suggested that the co-expression networks were significantly different at different stages. On the contrary, none of the co-expression networks were found to be associated with three accessions of upland cotton, indicating that a particular regulatory network of alkali-salt stress response in different upland cotton races. SS12 and SS48 stages were correlated with more than one module. SS12 stage was correlated with the brown 2 (Gene number: 1163, r = 0.79),, dark olive green (Gene number: 416, r = 0.76) and dark sea green 4 (Gene number: 237, r = 0.68) modules in leaf, and correlated with medium orchid (Gene number: 241, r = 0.66) and blue violet (Gene number: 2632, r = 0.71) modules in the root, respectively. SS48 stage was associated with misty rose (Gene number: 135, r = 0.67) and light slate blue (Gene number: 1300, r = 0.87) modules in leaf, and associated with light steel blue1 (Gene number: 1218, 0.87) and orange red 4 (Gene number: 536, r = 0.85) modules in the root, respectively. While, only one module was correlated with SS3 stage (navajo white (Gene number: 326, r = 0.88) moduleSS3 stage in leaf; orange red1 (Gene number: 1836, r = 0.95) module-SS3 stage in root). Interestingly, the higher and significant correlated modules were all positively correlated with different stages, indicated that different expression genes in different tissues/stages were predominantly up regulated to response to alkali-salt stress.
GO enrichment analysis
We performed GO enrichment analysis of genes in tissue/stage-specific modules. The enriched GO terms are presented in (Table S8). GO analysis categorizes the genes into three possible groups, namely cellular component (CC), molecular functions (MF) and biological process (BP). The molecular function detected for the DEGs were signal transducer activity (GO:0004871) and phosphorelay response regulator activity (GO:0000156) which were significantly enriched at the SS3 stage, protein tyrosine kinase activity (GO:0004713), protein kinase activity (GO:0004672) and protein serine/threonine kinase activity (GO:0004674) were top three highly enriched at the SS12 stage in leaf; whereas regulatory region nucleic acid binding (GO:0001067), transcription regulatory region DNA binding (GO:0044212), regulatory region DNA binding (GO:0000975) were top three significantly enriched at the SS48 stage (P < 0.001, FDR < 0.05) in leaf. In the GO enrichment analysis for the biological process (BP), no significant pathway was identified at the SS3 stages in leaf; protein phosphorylation (GO:0006468), macromolecule modification (GO:0043412) and protein modification process (GO:0036211) were the top three pathways at the SS12 stage in leaf; whereas single-organism metabolic process (GO:0044710), carbohydrate metabolic process (GO:0005975) and oxidation-reduction process (GO:0055114) were the top three pathways at the SS48 stage in leaf. For GO analysis of molecular function, nucleic acid binding transcription factor activity (GO:0001071), transcription factor activity, sequence-specific DNA binding (GO:0003700) and oxidoreductase activity acting on paired donors, with incorporation or reduction of molecular oxygen (GO:0016705) were top three significantly enriched at the SS3 stage in root; Additionally, UDP-N-acetylmuramate dehydrogenase activity (GO:0008762), chromatin binding(GO:0003682), macromolecular complex binding (GO:0044877) were top three highly enriched at the SS12 stage in root; whereas oxidoreductase activity (GO:0016491), Acyl-CoA dehydrogenase activity (GO:0003995) and catalytic activity (GO:0003824) were top three significantly enriched at the SS48 stage in root. Interestingly, oxidoreductase activity was found at different stages of differentiation. In the GO enrichment analysis of biological process, regulation of RNA biosynthetic process (GO:2001141), regulation of transcription, DNA-templated (GO:0006355), regulation of nucleic acid-templated transcription(GO:1903506) were top three significantly enriched GO terms in SS3 stage in root; regulation of RNA metabolic process(GO:0051252), regulation of RNA biosynthetic process GO:2001141and regulation of transcription, DNA-templated (GO:0006355)were the top three pathways at the SS12 stage in root. Similarly with the SS12 stage in leaf, single-organism metabolic process (GO: 0044710), single-organism process (GO: 0044699) and oxidation-reduction process (GO: 0055114) were the top three pathways at the SS48 stage in root. Overall, signal transduction was performed at SS3 and SS12 stages and regulates transcription was performed at SS48 stage in leaf. Although it is complicated to regulate alkali-salt stress in root, regulation of oxidation-reduction was performed in all stages. the results obtained were in agreement to previous findings in which similar GO terms have been identified for various stress responsive genes such as LEA genes , MATE genes  among others.
3.9 Identification of central and highly connected genes
To identify hub genes under the salt-alkali stress, two different methods were used for screening genes. Each gene possibility weight (p. weight) value was obtained by the network screening function of the WGCNA package based on gene significance (GS), modular membership (MM), and lower p. The weight value indicated that the gene was a higher correlation with traits (treatment stages). The top 30 hub genes were identified based on p. Weight. In the second method, the top 150 and top 300 highly connected genes for higher correlation modules of different stages were selected for analysis based on the topological overlap matrix (TOM) of all differential expression genes. In addition, the top 30 genes that were central and highly connected were visualized using the Cytoscape 3.3.0 software (Table S7, Figure 6A-B). In total, 180 hub genes of higher correlation with different treatment stages were found and 180 higher connectivity genes were also selected in six different stages, in which 39 common genes were identified by two different methods that indicated the core role of these genes in response to alkali-salt stress. Interestingly, 19 and 18 common genes were obtained in SS3L and SS48L stages, but none of the common genes was found in other stages. Transmembrane proteins (Pfam: DUF3082), Gh_D11G2953 and Gh_A11G2587, were the most correlated genes with other genes and at SS3L stages under alkali-salt stress. Galactosyltransferase family protein (Gh_D05G1401 and Gh_A05G1229) and lysine decarboxylase family protein (Gh_D05G1724 and Gh_A03G0267) were also observed an indication that these genes were significantly important factors at the SS3L stage in relation to salt stress tolerance.
Three different starch branching enzymes (Gh_A02G1739, Gh_D02G0995 and Gh_Sca006745G03) and two RING/U-box superfamily proteins (Gh_D05G0963 and Gh_D07G1649) were identified which could be playing an important role in SS48L stages. Homolog gene pairs were identified at the same time by WGCNA, suggested their potential key roles in the regulation of alkali-salt stress in cotton. A total of 35 TFs was found among hub genes, out of 26 TFs showed higher connectivity to various DEGs. Moreover, among the 35 TFs, 4 of the common genes had a higher correlation with other genes and treatment stages under alkali-salt stress, were members of the double B-box zinc finger (DBB, Gh_A10G0877 and Gh_A11G2610) and nuclear factor YA (NF-YA, Gh_D03G0606 and Gh_Sca006219G01) TF family. Overall, 321 hub genes were identified, including 35 transcription factors.
Validation of the hub genes by RT-qPCR
In order to detect the expression level and function of the hub gene, 12 genes were selected and their expression analyses. In each of the selected genes were randomly examined in five samples by qRT-PCR. The gene primer information of 12 hub genes was designed and details are contained in (Table S9). We found that the gene expression level data from RT-qPCR were significantly highly correlated (cor = 0.7397, p-value = 1.458e–11) with the data from RNA-seq (Figure 7A). This result proved the RNA-seq data were reliable in this research. Moreover, compared with the expression level of control seedlings, the hub genes were differentially expressed under salt stress conditions, which suggested that the hub genes identified by WGCNA had a putative role in salt-alkali stress response.
Enhanced salt sensitivity in GhSOS3 and GhCBL10 Virus-induced Gene Silencing (VIGS) seedlings
To further investigate the functions of hub genes, GhSOS3 and GhCBL10, the VIGS pYL156-GhPDS, pYL156-Ctrl, pYL156-GhSOS3 and pYL156-GhCBL10 plants were observed under the salt-alkali stress. Albino leaves were observed in pYL156-PDS inoculated seedlings after 7 days of inoculation. Compared with infected seedlings, we found that control seedlings had rapid growth after 20 days of inoculation. Moreover, the no differences were noted between infected seedlings (Figure 8A). The expression levels of GhSOS3 and GhCBL10 were checked by RT-qPCR Compared with pYL156-Ctrl seedlings, expression levels of GhSOS3 and GhCBL10 were down-regulated in the corresponding gene silencing seedlings after 20 days of inoculation (Figure 8B). The leave of pYL156-GhSOS3 and pYL156-GhCBL10 seedlings withered and wilted compared with the control and the pYL156-Ctrl seedlings after 20 days of salt-alkali stress treatment (Figure 8C). Additionally, the PRO and SOD content was lower in pYL156-GhSOS3 and pYL156-GhCBL10 seedlings compared with control and pYL156-Ctrl seedlings after 20 days of salt-alkali stress treatment. On the contrary, the MDA content was higher in pYL156-GhSOS3 and pYL156-GhCBL10 seedlings compared with control and pYL156-Ctrl seedlings (Figure 8D). This result suggested pYL156-GhSOS3 and pYL156-GhCBL10 seedlings sensitivity was enhanced.