Transcriptome analysis of upland cotton revealed novel pathways to scavenge reactive oxygen species (ROS) responding to Na2SO4 tolerance

Salinity is an extensive and adverse environmental stress to crop plants across the globe, and a major abiotic constraint responsible for limited crop production threatening the crop security. Soil salinization is a widespread problem across the globe, threatening the crop production and food security. Salinity impairs plant growth and development via reduction in osmotic potential, cytotoxicity due to excessive uptake of ions such as sodium (Na+) and chloride (Cl−), and nutritional imbalance. Cotton, being the most cultivated crop on saline-alkaline soils, it is of great importance to elucidate the mechanisms involved in Na2SO4 tolerance which is still lacking in upland cotton. Zhong 9835, a Na2SO4 resistant cultivar was screened for transcriptomic studies through various levels of Na2SO4 treatments, which results into identification of 3329 differentially expressed genes (DEGs) in roots, stems and leave at 300 mM Na2SO4 stress till 12 h in compared to control. According to gene functional annotation analysis, genes involved in reactive oxygen species (ROS) scavenging system including osmotic stress and ion toxicity were significantly up-regulated, especially GST (glutathione transferase). In addition, analysis for sulfur metabolism, results in to identification of two rate limiting enzymes [APR (Gh_D05G1637) and OASTL (Gh_A13G0863)] during synthesis of GSH from SO42−. Furthermore, we also observed a crosstalk of the hormones and TFs (transcription factors) enriched in hormone signal transduction pathway. Genes related to IAA exceeds the rest of hormones followed by ubiquitin related genes which are greater than TFs. The analysis of the expression profiles of diverse tissues under Na2SO4 stress and identification of relevant key hub genes in a network crosstalk will provide a strong foundation and valuable clues for genetic improvements of cotton in response to various salt stresses.

. These hormones mainly increases the plant salt tolerance by scavenging active oxygen 26 . In recent years, there have been a large number of studies on the relationship between ABA and plant salt tolerance. Among them, salt stress can induce the accumulation of endogenous ABA, and then induce the expression of E3 ubiquitin ligase or others transcription factors 25,27 , ion transporter 22 , antioxidant enzyme 28 and other salt stress-related genes to enhance plant salt tolerance 29 .
For ion transport, NHX1 regulates the export and import of Na + in and out of vacuoles 10 , and HKT1 transporters have been found to reduce Na + toxicity by regulating Na + /K + balance in several species 30 . In plant cells, the proton pump, including H + -PPase (proton pump pyrophosphoric acid hydrolase), V-PPase (vacuole proton pump ATP hydrolysis enzyme) and PMH-ATPase (plasma membrane proton pump ATP hydrolysis enzyme) on the one hand provides energy for cellular metabolism, maintain normal metabolism and cell growth 31 . On the other hand, proton pump such as H + -PPases can improve the salt-resistance 32 . In Arabidopsis, 14-3-3 proteins in the SOS system are able to regulate the V-PPase (EC 3.6.1.1) gene AVP1 and mitigate the harm of Na +33 . Besides, aquaporin protein, as a transporter of water molecules, responds to salt stress mainly by increasing antioxidant activity and maintaining ion balance 34 . The aquaporin protein related to salt resistance is mainly concentrated on the PIPs (Plasma membrane intrinsic proteins) gene. In Arabidopsis thaliana, Tip2 (Tonoplast intrinsic protein) regulates MDA (malondialdehyde) and salt-tolerance related genes (SOS 1 , SOS 2 , SOS 3 , DREB1A and P5CS1) 35 . GhSIP may be involved in endoplasmic reticulum osmotic equilibrium 36 .
Soil salinization has always been a major problem in global agriculture, because of inadequate irrigation and climate change, the area of saline-alkali land is likely to increase with the passage time 37 , so it is of great significance to study the salt tolerance in plants. The ions commonly present in saline soil are Na + , Cl − and SO 4 2− 38 . There are many studies on NaCl, however, few on Na 2 SO 4 . SO 4 2− is different from the Na + , K + and Cl − 39 , and Na 2 SO 4 is more toxic than NaCl 40,41 . Prosopis Strombulifera (Lam.) Benth, a kind of halophytic shrub with high NaCl tolerance, is found in high saline-alkali soil in Argentina, but affected with decreased Fv/Fm, ETR and Y(II) photosynthetic parameters significantly under Na 2 SO 4 condition 40 . While, in Kalidium foliatum the activity of Rubisco (Ribulose-1,5-bisphosphate Carboxylase/Oxygenase) treated with NaCl was higher than that of treated with 100 mM Na 2 SO 4 , and there was no significant change under NaCl + Na 2 SO 4 mixture treatment 42 . In Brassica rapa, the result of gene analysis of GSLs (Glucosinolates) biosynthetic pathway and transcription factor showed that sulfate had the strongest inhibition on growth under different treatment with NaCl, Na 2 SO 4 , KCl and K 2 SO 4 , respectively 41 .
So far, there is a lack of research on Na 2 SO 4 salt stress, especially in cotton, and the mechanism of its toxicity is still unclear. In this study, through the analysis of the phenotypic, physiological and biochemical indexes and morphological analysis of Zhong 9835 in response to salt stress during 0 h, 6 h, 12 h and 24 h, we found that it is best to analysis Zhong 9835 transcriptome at 12 h under 300 mM Na 2 SO 4 treatment. Sulfur metabolism was enriched in 3329 differentially expressed genes (DEGs) among roots, stems and leaves, especially GST, followed by Ubiquitin transcription factors. This study not only provides complement data for regulatory network at early stage under Na 2 SO 4 stress but also provides a strong foundation and valuable clues for genetic improvements of cotton in response to various salt stresses.

Results
Phenotypic responses of Zhong 9835 to Na 2 SO 4 stress. Previous studies have reported that cotton is more sensitive to abiotic stresses at three-leaf stage 43 . Different morphological data has been observed in G. hirsutum Zhong 9835 during its three-leaf stage under various concentrations of Na 2 SO 4 stress. We observed the significant phenotypic difference among roots, leaves and shoots of cotton seedlings under 300 mM Na 2 SO 4 stress after 12 h of treatment as compared to control conditions without any salt treatment (Fig. 1A). After 12 h of treatment, leaves start wilting, and roots gradually truing into black (Fig. 1A), while after 24 h, leaves and stems were seriously wilted (Fig. 1A). The result showed that under 300 mM Na 2 SO 4 solution the pH value ranges from 6.93 at 0 h to 5.6 at 24 h, among which it becomes more acidic after 6 h (Fig. 1B). The minimum pH value was noticed at 12 h to 24 h under the treatment of Na 2 SO 4 and control, especially in control from 7 to 6.5 (Fig. 1B). Meanwhile, SOD, POD, Proline and MDA contents in roots, stems and leaves all increased gradually from 0 to 24 h ( Fig. 1 C/D/E/F). Among which, SOD up to a same content among root, stem and leaf, especially largest increase from 6 to 12 h in root and stem (Fig. 1C). POD of root increased widely from 6 to 12 h (Fig. 1D). Substantial increase of proline content of root and stem was observed at 6 h to 12 h (Fig. 1E), MDA of leaf increased greatly from 6 to 12 h, while root amplitude increased from 0 to 6 h and stem from 12 to 24 h (Fig. 1F).
Transcriptome sequencing and alignment. Two groups (treated versus control) with three biological replications were conducted among the samples of roots, stems, and leaves, respectively. Totally 18 qualified libraries were established from the tissues of the roots, stems and leaves at 12 h with salt stress and control conditions. On average, of 42.5 million raw reads for the 18 libraries were obtained by using an Illumina Novaseq 6000 sequencing platform (Table 1). 111.34 Gb (Gigabyte) of sequence data and over 96% of the clean reads with a Q30 level were done through approximately 742 million clean valid reads. With the process of adaptor deletion, junk filtering and low copy filtering, > 95% of the sequences were confirmed as clean data, which then mapped to cotton whole genome (G. hirsutum) by using His at software 44 . In final, > 93.94% of the total reads were mapped to the reference genome, while the unique mapped reads were 63.46%-70.27% by the String Tie software 45 .
Exploration of DEGs in roots, stems and leaves in response to Na 2 SO 4 stress. Gene expression levels were estimated by fragments per kilo base of transcript per million fragments mapped (FPKM). Differential expression analysis of treatments and control group was performed using the DESeq among CK_R, SS_R www.nature.com/scientificreports/ CK_S, SS_S, CK_L, SS_L. According to root, stem, and leaf, transcriptome data between control and treatment with Na 2 SO 4 were divided into 3 groups (CK_R vs SS_R, CK_S vs SS_S, CK_L vs SS_L) ( Fig. 2A,B). The samples from Root (CK_R vs SS_R) showed 15,524 DEGs, among which 10,787 genes were up-regulated and 4,737 genes were down-regulated ( Fig. 2A). A total of 20,409 genes were differentially expressed in the samples of Stem (CK_S vs SS_S), among which 6,426 genes were up-regulated and 13,983 genes were down-regulated ( Fig. 2A). There were 12,146 DEGs identified from the Leaf sample (CK_L vs SS_L), among which 6,521 genes were upregulated and 5,625 genes were down-regulated ( Fig. 2A).

Validation of RNA-Seq data by quantitative real-time PCR.
In order to validate the differential expression analysis of RNA-seq data, we performed the quantitative real-time PCR (qRT-PCR) of 20 genes to confirm the reliability of RNA-seq data by using the same RNA samples (Table S1). To corroborate the expression levels measured by RNA-Seq data, the ratio of expression levels among root, stem and leaf under Na 2 SO 4 stress and control using RNA-Seq was compared to the ratios of expression measured by qRT-PCR. The results showed that there was a good correlation between RNA-Seq and real-time PCR results among three tissues (coefficient of determination R 2 = 0.86, 0.82 and 0.91) (Fig. 3A-C). The validation experiments support the accuracy of the RNA-Seq quantification of gene expression by relative values provided by the qRT-PCR analysis.
To study the clustering profiles of these 3329 DEGs, six samples with three biological repeats were carried out (Fig. 4). These 3329 DEGs were divided into six clusters, in which cluster had the same expression profile (Fig. 4). There were 630 genes in cluster 1, which had higher FPKM values especially in CK_L group with three biological repeats. In cluster 2, the samples SS_R possess 691 genes with three replication and harbor higher FPKM values. The FPKM values of cluster 3 were found higher than others in samples CK_S with three biological repeats, which includes 622 genes. Cluster 4 containing 320 genes was the smaller one in the clustering profile, in which the FPKM is greater in samples CK_R with three biological repeats particularly. Cluster 5 with 760 genes was the largest one in the clustering profile, which the FPKM of the samples SS_S with three biological repeats is higher. Cluster 6 of 306 genes was the smallest one, in which the FPKM of the samples SS_L with three biological repeats is greater.
Functional enrichment of differentially-expressed genes. To further understand the molecular mechanism of 3329 DEGs, we mapped all of DEGs to the GO database and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways 46 . In final, the identified 3,329 DEGs were classified into 100 pathways (1877 DEGs, account for 56.90%) and 383 Gene Ontology (GO) annotations (3144 DEGs, accounts for 95.30%) which includes the biological processes, cellular components and molecular functions. Then the first 30 GO terms (Fig. 5) and a threshold of top 30 set for KEGG pathways analysis (Fig. 6) were chosen.

Discussion
Na 2 SO 4 , being a neutral salt, is made up of Na + and SO 4 2− . High concentrations of Na + and SO 4 2− in the soil not only causes the salt toxicity to plants 50 , but also hindered the uptake of other minerals 51 . In our present study, we used Na 2 SO 4 solution of 300 mM as a salt stress to Zhong 9835 at three leaf stage, in which the pH value ranges from 6.93 at 0 h to 6.2 at 12 h (Fig. 1B). Under the mentioned treatment plants shows significant phenotypic differences among tissues of roots, shoots and leaves. As the roots get start blacking cotyledon turn wilted and leaves lost the turgidity due to loss of water and becomes weak and thin along with stem browning. (Fig. 1A). The pH value at 24 h was 5.6 ( Fig. 1B), the root blackened more significantly, the stem browned, the cotyledon shriveled, and the true leaves are wilting and turning brown (Fig. 1A). Previous studies report that the main site of Na + toxicity for most plants is the leaf blade, where Na + accumulates after being deposited in the transpiration stream, rather than in the roots 52 . And an important oxidative damage only to be induced by SO 4 2− anion with an increase in H 2 O 2 and MDA content 40 , although sulfur could be taken up by the roots and stored in the vacuoles of root and xylem parenchymal cells 53 . In other words, at 12 h root blackening is most likely the toxic phenomenon caused by the oxidation of SO 4 2− in weak acidic solution. Osmotic stress and ionic toxicity can cause oxidative damage 54 . In response to osmotic stress ( Figure S3), we found that starch and sucrose metabolism enriched (Fig. 6) and some organic material such as LEA (Fig. 9), HSFs (Fig. 8), proline and its biosynthesis key enzymes P5CS (Fig. 9) is up-regulated consistent with previous research report 49 . As ion transport for a Na + detoxification way ( Figure S5), SOS system, on the one hand, can potentially ejected Na + by Na + /H + exchangers located in the plasma membrane: on the other hand, sequestered it into the vacuole by Na + /H + exchangers (e.g. NHX proteins) located in the tonoplast 55 . As previous studies, HKT1, as a Na + /K + transporters, regulated the equilibrium of the Na + /K + decreasing of Na + toxicity 43 . Arabidopsis K + transporter1 (AKT1) activity is repressed by SCaBP8 (CALCINEURIN B-LIKE10 or CBL10, known as SOS3-LIKE CALCIUM-BINDING PROTEIN8, SCaBP8), which interacted with and activated by SOS 3 -SOS 2 complex (Fig. 11) 56 .
ROS scavenging system includes the enzymatic antioxidants (SOD, APX, PRX, GPX, CAT, GRX and TRX) and non-enzymatic scavengers (GSH-ASA system). It is well known that in cotton the enzymatic antioxidants and non-enzymatic scavengers participated in salt stress 57 . Among root, stem and leaf, not only the enzymatic www.nature.com/scientificreports/ antioxidants POD and CAT (Fig. 9), but also the non-enzymatic scavengers GPX, GST, DHAR and GST (Fig. 10B) genes were found to be up-expressed under Na 2 SO 4 condition. And it is interesting that the number of the non-enzymatic scavengers, as a role in the reaction of peroxide detoxification by catalyzing GSH to GSSG 58 , especially GST is more than the enzymatic antioxidants (Fig. 10B). These results suggest that GST may be the major scavengers of ROS under Na 2 SO 4 stress. In response to ROS ( Figure S1), both the ROS scavenging system and others signaling pathways regulated by transcription factors ( Figure S4) work synergistically with efforts in avoiding PCD ( Figure S2) 59 . Hormones including JA, ETH, BR, IAA and in particular ABA ( Figure S4) although just one gene with two sub-genomes (Fig. 7), under salt stress a mass of researches had reported the crosstalk of ABA and some TFs belonging to others class such as Ubiquitin, MYB, NAC, bZIP and AP2/ERF 60 . Among these transcription factors, ubiquitin is most than others (Fig. 8). In apple, MdbHLH3 gene (an anthocyanin-related basic helix-loop-helix transcription factor (bHLH TF) gene) promotes ethylene production involved in ethylene biosynthesis including MdACO1, MdACS1, MdACS5A MdACS1, and MdACS5A 61 . U-box E3 ubiquitin ligase gene MdPUB29, highly homologous with AtPUB29, direct ubiquitination of the MdbHLH3 protein, positively regulating salt tolerance 54 . In addition, VTC1-CSN5B associated with the COP9 signalosome complex promotes ubiquitination-dependent VTC1 degradation through the 26S proteasome pathway, affecting the response to salt stress by regulating ASA synthesis 62 . In hormones, ETH and ABA induced cell senesce or cell death 62,63 , while auxin, BR and JA as roles related to abiotic stress can prevent the accumulation of ROS from that 64 . For example, IAA/auxin can reverse the hypersensitive response stimulated by purified harpin protein to extent 65 . AUX/IAA-mediated auxin signaling contributes to ethylene-dependent inducible aerenchyma formation in rice roots 66 . To Arabidopsis cell suspension cultures, auxin had an effect on the control of cell wall composition and rigidity, preventing the cell death 67 . In present study, the largest number of hormones is IAA/auxin, followed by JA and BR (Fig. 7). These results indicate that in Zhong 9835 the root, stem, leaf cell of treatment of 12 h under Na 2 SO 4 stress is becoming balance by removing ROS toxicity (Fig. 11).   Mapping and differential expression genes analysis. Based on the quality results of the paired-end reads, we removed the low-quality reads by per script, which included only adaptor, unknown nucleotides > 5% and Q20 < 20% (percentage of sequences with sequencing error rates < 1%). After filtering from the raw reads, the clean reads were mapped to cotton genome (G.hirsutum) by Hisat software 44 . According to the mapped reads from the reference cotton genome, String Tie software 44 was used to estimate quantification of the gene expression levels with fragments per kilobase of transcript per million fragments mapped (FKPM) 68 . And an edger package, one of R packages, was applied for differential expression analysis between two groups with three tissues respectively. Fold Change ≥ 2 and FDR (false discovery rate) < 0.01 were taken as the threshold of the P-value in multiple tests for computing the significance of the differences. analysis implemented by GO seq R packages 69 , were divided into 3 classes, molecular function, cellular process, and biological process. KEGG enrichment analysis of the DEGs was applied to summary the statistical enrichment of differential expression of genes in KEGG pathways (http:// www. genome. jp/ kegg/) by KOBAS software 70 .

Validation of RNA-Seq by qRT-PCR.
Each sample with 3 biological replicates was performed by Realtime RT-PCR (qRT-PCR) 75,76 . A set of 20 genes was chosen randomly by the FPKM. Specific primers for the chosen genes were designed through Primer 3 software. cDNA was synthesized by using an EASY spin Plus Plant RNA Kit (TIANGEN) with gDNA Eraser (TaKaRa) 77 . The qRT-PCR reactions were conducted using a SYBR Green I Master mixture (Bio-Rad, CFX96, Switzerland) according to the manufacturer's protocol on a Light Cycler 480II system (Bio-Rad, CFX96, Switzerland). The results of qRT-PCR were analysed via the ΔΔCt method 78 , the cotton histone (His) 3 gene (GenBank accession no. AF024716) was used as a standardcontrol 79   . The leftmost part of the networks is related to SO 4 2− stress. The others part of the networks is related to ROS, which mainly include Na + stress and ABA signal. The vacuole is the crosstalk of SO 4 2− stress and Na + stress. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.