A Comprehensive Gene Co-Expression Network Analysis Reveals a Role of GhWRKY46 in Responding to Drought and Salt Stresses

Abiotic stress, such as drought and salinity stress, seriously inhibit the growth and development of plants. Therefore, it is vital to understand the drought and salinity resistance mechanisms to enable cotton to provide more production under drought and salt conditions. In this study, we identified 8806 and 9108 differentially expressed genes (DEGs) through a comprehensive analysis of transcriptomic data related to the PEG-induced osmotic and salt stress in cotton. By performing weighted gene co-expression network analysis (WGCNA), we identified four co-expression modules in PEG treatment and five co-expression modules in salinity stress, which included 346 and 324 predicted transcription factors (TFs) in these modules, respectively. Correspondingly, whole genome duplication (WGD) events mainly contribute to the expansion of those TFs. Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) analyses revealed those different modules were associated with stress resistance, including regulating macromolecule metabolic process, peptidase activity, transporter activity, lipid metabolic process, and responses to stimulus. Quantitative RT-PCR analysis was used to confirm the expression levels of 15 hub TFs in PEG6000 and salinity treatments. We found that the hub gene GhWRKY46 could alter salt and PEG-induced drought resistance in cotton through the virus-induced gene silencing (VIGS) method. Our results provide a preliminary framework for further investigation of the cotton response to salt and drought stress, which is significant to breeding salt- and drought-tolerant cotton varieties.


Introduction
With climate change and abnormal weather events, abiotic environmental factors, such as drought and salinity, restricted the growth and yield of crops worldwide [1,2]. The hormonal pathways, signal transcription, polysaccharide content, and lipid content in the plant are easily affected by abiotic stresses [3,4]. A series of reports demonstrated that drought and salinity stresses could affect the plant's secondary metabolites and the gene expression level by disrupting the homeostasis of the cell [5,6]. In these processes, the large amount of plant stress-resistant TFs, such as WRKY, AP2, VQ, MYB, NAC, MAPK, bZIP, and more, were detected and proved to be involved in the main pathway of the abiotic stress in Arabidopsis, rice, soybean, maize, cotton, and more [7][8][9][10]. The microarray analysis and RT-PCR results provide WRKY, ERF, and JAZ genes as potential markers of tolerance to salt stress in cotton [11]. By comparing the transcriptomics of the two cultivars, GhERF12 was identified as involving salinity tolerance during the early development of cotton [12]. By analyzing the long-reads RNA sequencing in cotton, the TFs were found to widely participate in the complex nature of salt stress tolerance mechanisms [13]. Moreover, the expression of hormone-related genes and salt stress-related genes in GaJAZ1 transgenic cotton were reprogrammed [14]. Furthermore, other research also provided the hormones and TFs important to the adaptability of cotton to abiotic stress [15][16][17]. Therefore, TFs played a central role in plant tolerance and adaptability.
In cotton, WRKY proteins also could play an important role in abiotic resistance. In G.aridum, GarWRKY5 was found involved in salt stress response related to activating the hormone signaling pathway [25] Similarly, GhWRKY6-like was reported as a negative regulator in response to salt stress via the ABA signaling pathway in cotton, and GhWRKY21 plays a role in the drought-induced ABA signaling pathway [26,27]. GhWRKY33 was found as a negative regulator involved in the drought stress in cotton [28]. GhWRKY41 might be a positive regulator of stomatal closure, and by regulating reactive oxygen species, enhance the plant tolerance to stress [29]. GhWRKY91 was also identified as a hub factor in the drought stress response in cotton [30]. WRKY proteins could also be involved in the process of leaf senescence. Overexpression of GhWRKY17 in Arabidopsis enhanced the plant's susceptibility to leaf senescence [31]. GhWRKY27 could positively regulate leaf senescence via interaction with GhTT2 and binding to the promoters of GhCYP94C1 and GhRipen2-2 [32]. In addition, WRKY TFs were also involved in fiber initiation and elongation [33].
WGCNA is proven to be an effective method for identifying cluster gene modules and hub genes in various crops [34,35]. In maize, the WGCNA method was used to exploit multiple traits to determine core modules and hub genes [36]. In other research on cadmium resistance, 22 regulatory modules were identified in maize [37]. Furthermore, there were lots of studies on the identification of hub genes related to fiber quality and resistance to stress response in cotton. In a recent report, Zou et al. reported five specific modules pertaining to fiber development at different growth stages [38]. Cheng et al. indicated 574 TFs and 936 hub genes related to cotton seedling cold resistance [16]. Moreover, a meta-analysis of cotton transcriptomics identified some hub genes, including RETICULONlike 5 (RTNLB5) and PRA1, involved in regulating stress responses [39]. At the seedling stage, a comprehensive analysis of two cotton genotype transcriptomics recovered the plant MAPK signaling pathway and diterpenoid biosynthesis involved in response to salt stress [40]. Additionally, a detailed analysis of the evolution and abiotic stresses in G. thurberi, G. klotzschianum, G. raimondii, and G. trilobum might provide available gene resources underlying a multi-abiotic-resistant cotton breeding strategy [41]. In G. arboretum, various tissue and stress-related transcriptomics were used to construct co-expression networks with over 500,000 pairs of edges and 33,413 nodes [42]. In addition, ccNET (http:// structuralbiology.cau.edu.cn/gossypium/) (accessed on 20 August 2020), COTTONOMICS (http://cotton.zju.edu.cn/index.htm) (accessed on 20 August 2020), MaGenDB (http:// magen.whu.edu.cn) (accessed on 20 August 2020) et.al websites were useful to the study of the co-expression functional analysis in cotton [43,44]. Therefore, the WGCNA can be used as a reliable method to estimate the gene function and further apply it to cotton breeding.
Cotton is one of the most important industrial crops, and it is grown for its elite fibers and oil for industries worldwide. However, the drought and salinity seriously limited the production and quality of cotton [45]. Thus, it is indispensable to mine the genes related to cotton's drought and salt stress. The completion and continuous update of cotton genome sequencing and transcriptome data made it possible to identify and exploit key genes related to cotton breeding [46][47][48]. In this study, the DEGs and WGCNA methods were performed to analyze the expression profile of previous transcriptomics data related to salinity and PEG-induced drought tolerance. Four and five co-expression modules related to drought and salt stress were constructed. Transcription factors, including WRKY, MYB, bHLH, and ERF proteins, were widely identified in these modules, distributed in every chromosome, and the WGD events mainly contributed to their expansion. In addition, a hub GhWRKY46 (GH_D07G1505) was isolated, which is located in the nucleus, and qRT-PCR assays indicate that GhWRKY46 responded to salt and PEG6000 stress resistance in cotton. Furthermore, silencing GhWRKY46 decreased the salinity and PEG6000 tolerance in cotton, suggesting that GhWRKY46 played a critical role in regulating the salinity and PEGinduced drought tolerance. This study will benefit the development of abiotic-resistant varieties in cotton breeding.

Transcriptome Sequences and Differential Expression Analysis
The transcriptome sequence related to the drought and salt treatment was downloaded from the SRA database to identify genes related to abiotic stress in cotton. A total of 1754.18 million raw reads and 1678.94 clean reads were obtained, with an average of 26.23 million reads per sample. For each sample, the GC content of the clean reads was from 43% to 45%, and the high-quality reads that mapped to the G.hirsutum reference genome were ranging from 97.11% to 98.93% (Table S1).
In this study, the data sets include the transcript with expressed levels, with FPKM > 1 in at least three samples. A total of 25,655 and 24,542 genes related to drought (PEG) and salinity (NaCl) stresses were obtained, respectively (Figure 1a,b; Table S3). The largest number of DEG (1899) was identified in PEG at 3 h, and the smallest number of DEGs (1543) was identified in PEG at 24 h (Figure 1a,c; Table S3). However, under salt stress, it was estimated that the amount of DEG was the highest (1911) at 1 h, and the lowest amount of DEG (1711) was obtained at 6 h (Figure 1b,d; Table S3). As is depicted in the Venn diagram, 3494 DEGs were detected in the two stresses (Table S3; Figure S1).

Gene Co-Expression Construction and Analysis
To reveal the potential regulatory pathways for resistance to drought and salt stress in cotton, we constructed the co-expression modules through the WGCNA method. In this study, we selected the weight value β = 18 to construct the scale-free networks, describing different modules with different colors and merging similar modules. These modules were defined as clusters of highly interconnected genes, and genes within the same cluster have high correlation coefficients and potential functional relations. Four co-expression modules related to the PEG treatment were constructed. The turquoise module (2690) was with the maximum count, and the yellow module was with the minimum (119) (Figure 1e; Table 1). Under the PEG treatment, the range of the up-regulated DEGs in different modules was from 9 (7.56%, the yellow module) to 461 (17.14%, the turquoise module), and the downregulated number varied from 15 (12.61%, the yellow module) to 196 (7.29%, the turquoise module). Furthermore, clustering analysis suggested five co-expression modules related to salinity treatment. The largest module was the turquoise module, which contained 1900 genes, and the smallest module was the green module, including only 64 genes ( Figure 1f and Table 1). The up-regulated DEGs in the different modules varied from 4 (6.2%, the green module) to 265 (13.95%, the turquoise module), while the down-regulated varied from 14 (4.55%, the green module) to 266 (19.57%, the blue module) ( Table 1).
We further analyze duplication events of those TFs to explore their expansion mechanism. In the modules related to PEG and salt treatment, the duplication type of these TFs was identified (Table 2). For the PEG-treatment modules, there were 1, 3, and 132 TFs related to dispersed, tandem, and WGD events in the blue module, respectively; 1, 3, and 128 TFs were detected relating to dispersed, tandem, and WGD events in the turquoise module, respectively (Table S6). Meanwhile, in the salt stress modules, 2, 2, and 66 TFs were related to dispersed, tandem, and WGD events in the brown module, respectively. In the blue module, 1, 3, and 93 TFs might expand by dispersed, tandem, and WGD events, respectively; 2 dispersed and 95 TFs were detected relating to dispersed, tandem, and WGD events in the turquoise module, respectively (Table S6). Therefore, WGD events mainly contributed to the expansion of identified TF in the modules related to salt and PEG-induced drought stress.

Candidate Module Identification and Functional Analysis
With the selected correlation value of the |r| > 0.7, four modules were found in the PEG stresses, including the yellow module (r = 0.8, p = 0.1) at 6 h, the blue module (r = 0.85, p = 0.07) and the brown module (r = 0.83, p = 0.09) at 12 h, and the turquoise modules with a negative correlation (r = −0.72, p = 0.2) at 24 h ( Figure S2   To uncover the potential function in the above modules, we performed KEGG enrichment analysis in different modules in this study. KEGG analysis of the four modules related to PEG treatments indicated that the potential pathways were enriched in the metabolism of terpenoids and polyketides, protein export, ubiquitin mediated proteolysis, and metabolism of cofactors and vitamins (Figure 3a-d). Among the pathways related to salt treatment, the metabolism pathway is the vast majority pathway, such as peroxisome, valine, leucine and isoleucine degradation, pentose phosphate, and thiamine metabolism (Figure 3e-i). Additionally, GO enrichment analysis was also performed in this study (Table S8). In the PEG treatment analysis, the GO enrichment in the blue module suggested those genes were mainly involved in transferase activity, transferring hexosyl groups, oxidoreductase complex, and oxidoreductase complex. The genes related to cell metabolism, macromolecular metabolic process, heat shock protein, and water absorption were found in the brown module and are essential to cotton drought adaptation. At the same time, the turquoise module was enriched in peptidase activity, aspartic peptidase activity, oxidoreductase, ion transport, defense mechanisms, and other multi-biological processes. Correspondingly, GO analysis in the salt treatment of the blue module was related to material transport, transmembrane transport, ion steady state, and cell walls. Additionally, in the turquoise module, the electronic signal transmission, abscisic acid reaction, reactive oxygen response, and secondary metabolites were enriched. The gene in the brown module is mainly related to cell stimulation and toxicity response, such as peroxides and superoxide. Furthermore, in the yellow module, the genes were enriched in the organic acid metabolic process, carboxylic acid metabolic process, monocarboxylic acid biosynthetic process, and fatty acid biosynthetic process; and the genes in the green module were enriched in the chloroplast, plastid, anatomical structure development, and developmental process (Table S8). Detailed function annotations make it better to understand the function of distinguishing gene lists in different modules through WGCNA and multiple gene function annotation analysis.

Identification of Hub Genes and Gene Expression Assays
The correlation network and hub genes were further constructed and identified by preforming the CytoHubba package in Cytoscape software ( Figure 4). In the gene expression regulation network, hub genes might interact with other genes through interaction or regulation. In modules related to PEG stress, three hub TFs were identified in the blue module, including WRKY (GH_D07G1505), HAT (GH_D11G0270), and SCL (GH_D12G1100) (Figure 4a); seven hub TFs were also identified in the turquoise module, including MYBS (GH_A09G0633), MYB (GH_A09G1143), WHIRLY (GH_D08G2612), bHLH (GH_A11G1316), BBX (GH_A03G1868), and COL (GH_A12G0567 and GH_A01G2011) (Figure 4b). Moreover, in the module related salt stress, three hub TFs were identified in the blue module, which contained WRKY (GH_A02G0035), MYBS (GH_A09G0633), and REM16 (GH_D01G1221) (Figure 4c). In the brown module, three hub genes were identified, including KAN2 (GH_D08G1819), HSFA8 (GH_A12G1737), and RAP2 (GH_D06G0186) (Figure 4d). In the turquoise modules, a total of 5 hub TFs, including COL genes (GH_A01G2011 and GH_A09G0650), a GRAS gene (GH_D12G1100), P450 (GH_D01G0390), and MADS-box gene (GH_D13G2046) were identified (Figure 4e). These mined modules might provide a potential functional connection between the hub genes and other functional genes.
Next, 15 hub TFs were selected for further qRT-PCR analysis to confirm their potential function. The qRT-PCR results suggest that TFs were involved in response to stress resistance at different stages in cotton ( Figure 5; Table S7). Some selected genes had a high expression at 12 h and 24 h under the abiotic stress, especially the genes including WRKY TFs (GH_A02G0035 and GH_D07G1505), MADS-box TF (GH_D13G2046), and COL TFs (GH_A01G2011 and GH_D01G2107) (Figure 5a; Table S7). Almost all of the related genes' promoter sequences contain more than one WRKY binding site, which might be regulated by the WRKY proteins (Table S9). Therefore, we speculate that these hub genes might play a role in responding to salt and PEG-induced drought stress through participating in the pathway related to WRKY proteins.

Identification of Hub Genes and Gene Expression Assays
The correlation network and hub genes were further constructed and identified by preforming the CytoHubba package in Cytoscape software (Figure 4). In the gene expression regulation network, hub genes might interact with other genes through interaction or regulation. In modules related to PEG stress, three hub TFs were identified in the blue module, including WRKY (GH_D07G1505), HAT (GH_D11G0270), and SCL (GH_D12G1100) (Figure 4a); seven hub TFs were also identified in the turquoise module, including MYBS (GH_A09G0633), MYB (GH_A09G1143), WHIRLY (GH_D08G2612), bHLH (GH_A11G1316), BBX (GH_A03G1868), and COL (GH_A12G0567 and GH_A01G2011) (Figure 4b). Moreover, in the module related salt stress, three hub TFs were identified in the blue module, which contained WRKY (GH_A02G0035), MYBS (GH_A09G0633), and REM16 (GH_D01G1221) (Figure 4c). In the brown module, three hub genes were identified, including KAN2 (GH_D08G1819), HSFA8 (GH_A12G1737), and RAP2 (GH_D06G0186) (Figure 4d). In the turquoise modules, a total of 5 hub TFs, including COL genes (GH_A01G2011 and GH_A09G0650), a GRAS gene (GH_D12G1100), P450 (GH_D01G0390), and MADS-box gene (GH_D13G2046) were identified (Figure 4e). These mined modules might provide a potential functional connection between the hub genes and other functional genes.

Silencing of GhWRKY46 Decreased Salt and PEG-Induced Drought Tolerance
In order to verify the accuracy of data analysis and the function of the hub gene, we selected GhWRKY46 for further functional analysis. The qRT-PCR results show that GhWRKY46 participated in the salt and PEG-induced drought responses. As a homologous gene of AtWRKY46 (AT2G46400) in Arabidopsis ( Figure S3), GhWRKY46 contained a conserved WRKY domain and had three exons and two introns ( Figure S3). Next, by performing the subcellular localization assay, we found that GhWRKY46 can be transported into the nucleus of N. benthamiana cells, which suggests that GhWRKY46 might perform its functions in the nucleus ( Figure 6).
The cotton lines transformed with the pTRV2::CLA1 show an albino phenotype, indicating that the VIGS experiment was successful (Figures 7b and 8b). The expression level of GhWRKY46 in the silent plants was significantly lower than in the pTRV2::00 plants (Figures 7c and 8c). Three weeks later, the silenced and control plants were irrigated with a 400 mM NaCl solution. The results indicate that the silenced plants show a salt-sensitive phenotype after treatment for two days compared with the control (Figure 7a). The MDA content in pTRV2::GhWRKY46 plants was significantly higher than that in pTRV2::00 plants (Figure 8d). For the PEG6000 treatment, we also found that the wilting was more apparent in the leaves of pTRV2::GhWRKY46 plants than in the pTRV2::00 plants under the PEG6000 treatment ( Figure 8a). Additionally, the MDA content was higher than that in the control plants (Figure 8d). Our results prove that silencing of the GhWRKY46 can reduce cotton tolerance to salt and PEG-induced drought stresses.  Next, 15 hub TFs were selected for further qRT-PCR analysis to confirm their potential function. The qRT-PCR results suggest that TFs were involved in response to stress resistance at different stages in cotton ( Figure 5; Table S7). Some selected genes had a high expression at 12 h and 24 h under the abiotic stress, especially the genes including WRKY TFs (GH_A02G0035 and GH_D07G1505), MADS-box TF (GH_D13G2046), and COL TFs (GH_A01G2011 and GH_D01G2107) (Figure 5a; Table S7). Almost all of the related genes' promoter sequences contain more than one WRKY binding site, which might be regulated by the WRKY proteins (Table S9). Therefore, we speculate that these hub genes might play a role in responding to salt and PEG-induced drought stress through participating in the pathway related to WRKY proteins.

Silencing of GhWRKY46 Decreased Salt and PEG-Induced Drought Tolerance
In order to verify the accuracy of data analysis and the function of the hub gene, we selected GhWRKY46 for further functional analysis. The qRT-PCR results show that GhWRKY46 participated in the salt and PEG-induced drought responses. As a homologous gene of AtWRKY46 (AT2G46400) in Arabidopsis ( Figure S3), GhWRKY46 contained a conserved WRKY domain and had three exons and two introns ( Figure S3). Next, by performing the subcellular localization assay, we found that GhWRKY46 can be transported into the nucleus of N. benthamiana cells, which suggests that GhWRKY46 might perform its functions in the nucleus ( Figure 6). The cotton lines transformed with the pTRV2::CLA1 show an albino phenotype, indicating that the VIGS experiment was successful (Figures 7b and 8b). The expression level of GhWRKY46 in the silent plants was significantly lower than in the pTRV2::00 plants (Figures 7c and 8c). Three weeks later, the silenced and control plants were irrigated

Discussion
Cotton is an important economic crop with natural abiotic resistance and is widely planted in the world. Abiotic stresses, including drought and salinity stress, affected the cotton's growth and restricted the planting area in the world [1,2]. Therefore, understanding the complicated potential mechanisms and serious pathways, including genes related to hormone signal transcription, gene expression, peptide, and physiological indicators, will better explore the resistance mechanism and improve the existing crop varieties [3,5,6]. Studies in Arabidopsis, rice, cotton, and more, show that transcriptomics sequencing is a fast and effective method to obtain candidate genes and predict functional regulation pathways [44,[49][50][51][52][53]. Here, we performed a comprehensive analysis of the transcriptomics data under the PEG and salt stresses in cotton and explored the potential network in upland cotton.

DEGs, Co-Expression Network, and Polyploidization Event Analysis
The R/edgeR and R/WGCNA packages were widely used to explore DEGs and core gene-related traits or stress-tolerant genes and related mechanisms in plant transcriptomics. In our study, a large number of drought-responsive and salinity-responsive genes were identified, and their up-regulated or down-regulated DEGs were clearly displayed (Figure 1c,d and Figure S1; Table S3). These results indicate that induced gene expression, at different periods, is various, and different gene sets were activated in response to salt and drought stress. In the current study, we found 25,655 and 24,542 DEGs in the drought and salinity stresses, and 3494 DEGs were detected in both DEGs ( Figure S1; Table S3). Followed by the WGCNA analysis, four modules were found in PEG stress and five modules in salinity stresses (Figure 1, Table 1). By performing the TFs prediction, WRKY, MYB, bHLH, ERF, and more, domain-containing genes were found in various modules.
The polyploidization event significantly contributed to the plant's adaption to environmental changes and led to the expansion of plant genomes and gene numbers [54,55]. Given that the WGD event is the main factor that doubles the plant genome and promotes stress resistance adaptation, research on those genes was essential to uncover the potential regulatory mechanism [56]. We further found that the TFs expanded by the WGD event were counted at 97.36% and 96.21% in modules related to PEG and salinity stresses, respectively ( Figure S1; Table S3). Additionally, the WRKY domain-containing genes in this study were displayed as 28 (10.56%) and 27 (10.23%) in the PEG and salinity stress modules, respectively. Therefore, we speculated that WRKY proteins might be the main factor in abiotic stress resistance in cotton.

Gene Enrichment Analysis and Candidate Gene Identification
Antioxidant and transporter activity are essential protective mechanisms that protect the plant from abiotic stress. A study on the salt-tolerant genotype, Zhong9807, showed that GO terms were mainly enriched in catalytic activity, transporter activity, and antioxidant activity, and the KEGG were mainly enriched in hormone synthesis related, ROS related, and hormone signal transduction related pathways [13]. In addition, genes associated with "response to oxidative stress," "antioxidant activity," and "peroxidase activity" were significantly enriched in salt-tolerant and sensitive cotton genotypes [12]. Other abiotic stress research also indicated that "signal transduction" and "secondary metabolite biosynthesis", and more, pathways are essential to plant growth and adaptable development [13,57]. In the current study, in order to understand the function of PEG-related and salinity-related genes in cotton, the different modules' gene list was also analyzed by KEGG and GO enrichment. For both drought and salinity stresses, our study found enrichment of the KEGG pathway associated with "terpenoids and polyketides", "protein export", "ubiquitin mediated proteolysis", "metabolism of cofactors and vitamins", "fatty acid degradation", "phagosome", "carbohydrate metabolism", "circadian rhythm" and "thiamine metabolism" (Figure 3). Moreover, "transferase activity", "cell wall", "inorganic ion homeostasis", "photosystem", "biological regulation", "signal transduction", "hormone signaling (abscisic acid and ethylene)", and more, were also widely identified in GO terms (Table S8). Taken together, signal transport and hormonal pathway responses were related to the PEG and salinity stress in cotton, which corresponded to previous research.
The above analysis indicates that the hub gene might potentially regulate the abiotic stresses, especially the TFs, as the candidate stress-related genes set (Tables S3 and S4). As key regulators of abiotic stresses, WRKY TFs play a critical role in broad stress adaptation. Previous studies prove that GhWRKY17, GhWRKY21, GhWRKY27, GhWRKY33, GhWRKY41, GhWRKY70, and other GhWRKYs could regulate the resistance to salt, drought, verticillium wilt, and other abiotic stresses, respectively [27][28][29]31,32,58]. Moreover, other essential regulators of stress tolerance TFs, including MYB, ERF, bHLH, NAC, and bZIP, were also identified and selected to construct the co-expression network (Figure 4; Table S4). In addition, the analysis of their promoter sequence indicated a complicated regulation in stress tolerance in cotton (Table S9). Altogether, identification of the TFs in stresses and duplication analysis shed a new light on the significant contribution to cotton adaption in multi-abiotic stresses.

Silencing GhWRKY46 Enhanced the Sensitivity to Salinity and Drought in Cotton
Previous studies show that AtWRKY46, the homolog gene of GhWRKY46, plays dual roles in regulating plant response to drought and salt stresses. It interacts with AtWRKY50/70 as a signaling component involved in BR-regulated growth and drought responses [23]. Here, we made a comprehensive co-expression analysis and observed a hub gene, GhWRKY46, and we further found that GhWRKY46 could respond to salt and PEG6000 treatment. In addition, the silencing of GhWRKY46 enhanced sensitivity to salinity and drought in cotton (Figures 7a and 8a). These results indicate that GhWRKY46 might participate in regulating salinity and drought stress in cotton. Previous studies show that the MDA content was related to oxidative stress and redox signaling, particularly in plant abiotic stresses, and an indicator of ROS-dependent cell damage [59]. Our results here show that MDA content experienced a significant change in the GhWRKY46 silencing plants under the treatment of salt and PEG6000, indicating that GhWRKY46 may contribute to salt and drought stress response through regulating ROS scavenging (Figures 7d and 8d).

Acquisition and Comparison Analysis of Cotton Transcriptome Data
A transcriptome project (PRJNA490626) that contained 32 transcriptomes from previous research from the NCBI SRA (Sequence Read Archive) database was downloaded, which related to the cotton seedling treated with sodium chloride and PEG stress (Table S1) [46]. FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) (accessed on 16 July 2020) and Trimmomatic (version 0.3.9) were used to perform sequencing quality evaluation and low-quality read filtration, respectively [60]. HISAT2 was used to build the G.hirsutum genome (ZJU2.1 version) index file [46,61]. Samtools (version 1.9) and featureCounts (version 1.5.3) were used for data format conversion and calculating the gene FPKM (Fragments per Kilobase per Million) value, respectively [62,63].

DEGs Analysis and Gene Co-Expression Construction
The DEGs of PEG and salinity treatments post of 1 h (h), 3 h, 6 h, 12 h, and 24 h were identified by edgeR (R version 3.10) [64]. Genes with | logFC | > 1 and p-value < 0.05 were selected as DEGs in this study. The WGCNA (version 1.69) package was used to construct the weighted gene co-expression network, divide the relevant modules, and select hub genes [65]. The weight value was calculated by pickSoftThreshold in the WGCNA package, and β = 18 was selected to perform power processing to obtain a scale-free adjacency matrix on the original scaled relationship matrix. The topological disparity matrix (dissTOM = 1-TOM) and the dynamic shearing algorithm were used to classify gene clustering and the module division. The minimum number of genes in the module is 30 (min-ModuleSize = 30), and the merge threshold of similar modules is 0.25 (cutHeight = 0.25).
Furthermore, the CytoHubba package in the Cytoscape software (version 3.7.2) was used to visualize the network in the modules [66,67].

RNA Extraction and qRT-PCR Analysis
Upland cotton cultivar TM-1 seeds were grown in the growth room. Seedlings with uniform growth at the three-leaf stage were treated with 400 mM PEG and 400 mM NaCl, respectively. Leaf samples were collected at 0 h, 1 h, 3 h, 6 h, 12 h, and 24 h after treatment and rapidly saved in liquid nitrogen and stored at −80 • C in the refrigerator. The RNA extraction kit (Polysaccharides & Polyphenolics-rich, DP441) and the Takara reverse transcription kit (Mir-X TM MIRNA First-Strand Synthesis Kit) were used for RNA extraction and RNA reverse transcription. The Roche LightCycler 480 System (Roche, Germany) with the Cowin Bio and UltraSYBR One-Step Fluorescence Quantitative PCR Kit (UltraSYBR One-Step RT-qPCR Kit) were used to perform qRT-PCR. The GhUBQ7 gene was selected as the internal reference gene. Primers for qRT-PCR are shown in Table S7. The reaction procedure is: 95 • C for 10 min preheat denaturation, 95 • C for 5 s, 60 • C for 15 s, 72 • C for 10 s, 40 cycles; melting curve programs are: 95 • C for 15 s, 60 • C for 1 min, 95 • C for 15 s; reaction system is: 2xUltra SYBR Mixture 10 µL, forward primer (10 µmol L −1 ) 0.6 µL, reverse primer (10 µmol L −1 ) 0.6 µL, cDNA 0.8 µL, and ddH2O 8 µL. Three biological replicates were taken for each sample, and three independent experiments were performed. The results were calculated using the relative quantitative method 2 −∆∆Ct [73].

Subcellular Localization
The full-length CDS of GhWRKY46 was amplified from the upland cotton cultivar TM-1 and cloned into the pBI121-GFP vector. The leaves of six-week-old N. benthamiana leaves were used to inject pBI121-GFP, DAPI, and GhWRKY46-GFP, respectively. After the injection, the N.benthamiana plants were treated with dark for 24 h, then exposed to light treatment for 48 h. Observations under the laser confocal microscope were recorded.

Virus-Induced Gene Silencing of the GhWRKY46 in Cotton
A 300-bp fragment of GhWRKY46 was amplified and cloned into the pTRV2 (pYL156) vector to produce pTRV2::GhWRKY46 constructs, and the primers were listed in Table S2. The recombinant construction vector was transformed into the Agrobacterium tumefaciens strain LBA4404. The cotyledons of TM-1 cotton seedlings were used to inject an equal amount of Agrobacterium expressing the vectors, including pTRV2::00 (empty vector), pTRV2::GhWRKY46, pTRV2::CLA1 (positive control), and pTRV1 (pYL192, helper vector), as previously described [12,25,31]. Three weeks later, the plants with pTRV2::00 and pTRV2::GhWRKY46 were subjected to salt and PEG6000 stress. The malondialdehyde (MDA) contents were measured to determine the degree of damage to cotton leaves according to the standard methods (Solarbio, Beijing, China).

Conclusions
In the present study, a comprehensive analysis of transcriptomic of PEG and salinity stresses was performed in cotton. The DEGs and co-expression analysis showed differences in the number of genes contained in each module. Most of the TFs belonged to the WGD events in the gene expansion analysis. Moreover, KEGG and GO analysis proved that the peptidase activity, transporter activity, and lipid metabolic processes were critical to cotton abiotic stresses. Several hub genes contained within network modules were associated with abiotic stresses. Moreover, qRT-PCR assays demonstrated that numerous hub genes were further proved to respond to the salt and PEG-induced drought stress, including GhWRKY46. As a TF, GhWRKY46 plays its role at the nucleus. In addition, the VIGS assays and the measurement of MDA proved that the GhWRKY46 plays a positive role in the salt and PEG-induced drought stress. These results provide valuable information for further research investigating the salt and drought tolerance in cotton and provide a new gene resource for future breeding.

Conflicts of Interest:
The authors declare no conflict of interest.