All G. hirsutum plant materials used in the research were TM-1 (Texas Marker-1, the upland cotton genetic standard line) background. For subsequent quantitative reverse transcription (qRT)-PCR experiments, seeds of TM-1 were germinated and planted in soil under the following conditions: 12000 Lux light 16 h at 25 °C/dark 8 h at 23 °C, 80% humidity for 14 days (the first true leaf appeared). The seedlings above were divided into two groups and watered by 1/2 MS nutrient solution as the control or by 1/2 MS nutrient solution with 500 mM NaCl for 12 h.

The genome datasets of G. hirsutum (ZJU, version 2.1) and G. barbadense (ZJU, version 1.1) were downloaded from COTTONOMICS (http://cotton.zju.edu.cn/index.htm), G. arboreum (WHU, version 3.0) and G. raimondii (NSF, version 1.0) from CottonGen (https://www.cottongen.org) and Arabidopsis from TAIR 10 (http://www.arabidopsis.org/). Other plant species genome datasets were downloaded from Phytozome v12.1 (https://phytozome.jgi.doe.gov/). The databases of PlantTFDB (http://planttfdb.cbi.pku.edu.cn/) and SMART (http://smart.embl.de/) were used to confirm the conserved RING-linke zinc-finger domain (DUF702). The molecular weight (MW), isoelectric point (pI) of each GhSRS were calculated using the ExPASy (https://web.expasy.org/) compute pI/Mw tool, Plant-Ploc (http://www.csbio.sjtu.edu.cn/bioinf/plant/) was for subcellular localizational prediction (Chou and Shen, 2007). All gene names and their IDs were listed in Supplementary Table 1 .

SRS family amino acid sequences were used to perform multiple alignments by MEGA X software (Kumar et al., 2018) with MUSCLE default parameters, and then visualized using DNAMAN v7. Furthermore, the rooted and unrooted phylogenetic trees were constructed using MEGA X with Neighbor-Joining (NJ) methods, and 1000 bootstrap replicates were used to test reliability in each node and the maximum likelihood (ML) tree was constructed using “One Step Build a ML Tree” plugin in Tbtool software (Chen et al., 2020) with 5000 bootstrap replicates.

The gene structures of GhSRS genes were inferred by corresponding coding sequences. The MEME (https://meme-suite.org/meme/) program was used to identify conserved motifs in GhSRS proteins (Bailey et al., 2009). Tbtools Gene Structure View was used to draw the exon-intron structure, conserved motifs, and DUF702 domain distribution (Chen et al., 2020). The MCScanX software was used to analyze SRS protein sequences synteny and collinearity relationship between G. hirsutum, G. arboreum, and G. raimondii.

Chromosomal positions of GhSRS genes were obtained from gff annotation files for G. hirsutum (ZJU, version 2.1). The synteny and collinearity analysis between G. hirsutum, G. arboreum and G. raimondii were employed by the MCScanX software (Wang et al., 2012). Tbtools Gene Location Visualize and Advanced Circos were used to draw the distribution of GhSRS genes chromosomal mapping and synteny relationships (Chen et al., 2020).

To analysis the expression patterns of SRS genes under abiotic stress, a high-through RNA-seq datasets of leaf tissue under control and two types of stress (NaCl and PEG) treatments were obtained from COTTONOMICS (http://cotton.zju.edu.cn/) (Hu et al., 2019). Fragments per kilobase of exon per million fragments mapped (FPKM) was used for the quantification of gene expression. The clustered heatmap was drawn and normalized by the average expression levels (log2) based on FPKM values.

349 bp fragment of coding DNA sequence (CDS) of GhSRS21 was amplified from G. hirsutum (TM-1) cDNA and constructed into pTRV2. Then, the recombinant vector above, pTRV1 vector and CLA-pTRV2 were transferred into Agrobacterium strain GV3101, respectively, for subsequent experiments.

Seeds of TM-1 were germinated and planted in soil under the following conditions: 12000 Lux light 16 h at 25 °C/dark 8 h at 23 °C, 80% humidity, until the cotyledons were fully opened (about 10 d after seeds germinated). The Agrobacterium strains containing the GhSRS21-pTRV2 and pTRV1 or CLA-pTRV2 and pTRV1 or pTRV2 and pTRV1 plasmids were mixed and injected into TM-1 leaves as experimental groups, positive and negative controls, respectively. Cotton infected by Agrobacterium tumefaciens was cultured under the conditions above until the white striped leaf or albino phenotype appeared in the positive control (about 10 d after infected). Half of the experimental and negative control groups was watered by 1/2 MS nutrient solution as the control while the other half was watered by 500 mM NaCl regularly after every 3 days until the phenotypes appeared (about 14 d after treatment).

Total RNA was isolated from leaf tissues of TM-1 under salt stress treatment at 0 h, 1 h, 3 h, 6 h and 12 h using RNAprep pure Plant Kit (code: DP432, Tiangen, Beijing, China). First-strand cDNA was synthesized using Hifair 1st Strand cDNA Synthesis SuperMix for qPCR (YEASEN Biotech, Shanghai, China). The qRT-PCR was performed using Hieff qPCR SYBR GreenMaster Mix (YEASEN Biotech, Shanghai, China). The cotton ubiquitin gene UBQ7 was used as the internal control for the relative expression calculation. RNA isolation and qRT-PCR was manipulated based on the manufacturer’s instructions. The primers used for qRT-PCR were listed in Supplementary Table 2 .

The samples of cotton leaves obtained from negative control groups and GhSRS21 VIGS lines (experimental groups) treated with or without salt and PEG6000 in method 2.6 were used for measurement of chlorophyll content and determination of malondialdehyde, H₂O₂, O^2·- and antioxidant enzymes activities. The detail method for determination of malondialdehyde, H₂O₂, O^2·- and antioxidant enzymes activities was performed as described by Zhan et al. (Zhan et al., 2021), and method for measurement of chlorophyll contents were performed according to our previously described methods (Liu et al., 2019).

The CDS of GhSRS21 fused with GFP (green fluorescent protein) tag at the 3′-end was ligated into the pCAMBIA1300 vector. The recombinant plasmid above was mixed with the nuclear marker NLS-mCherry and co-transformed into Arabidopsis mesophyll protoplasts as described (Yoo et al., 2007). The protoplasts were observed and photographed by a fluorescence microscopy (Zeiss Imager.A2, Germany). The primers used above were listed in Supplementary Table 2 .

Total RNA was isolated from leaf tissues of negative control groups and GhSRS21 VIGS lines when the albino-like appearance on the leaves of positive control was observed. The fragments were purified by agarose gel electrophoresis and sequenced with NovaSeq 6000 Sequencer (Illumina Inc., San Diego, CA, USA) with a read length of 150 bp. Three biological replicates were performed separately. The raw data were filtered with fastp (https://github.com/OpenGene/fastp) (Chen et al., 2018b). The reads filtered above were then mapped to the cotton reference genome using HISAT2 software (http://ccb.jhu.edu/software/hisat2) (Kim et al., 2015; Zhang et al., 2015). Transcript analysis was performed using StringTie (https://ccb.jhu.edu/software/stringtie/) (Pertea et al., 2015), and differential expression genes (DEGs) analysis was performed by DEseq2 (Love et al., 2014), FPKM > 1.0 (FPKM, Fragments Per Kilobase of exon per Million mapped reads) were regarded as valid DEGs. Subsequently, KEGG annotations were performed using the online software EggNOG-Mapper (http://eggnog-mapper.embl.de/) and KEGG enrichment analysis were performed using Tbtools software (Chen et al., 2020).

Statistical analysis was performed using SPSS version 23.0 statistical software (SPSS, Inc., Chicago, IL, USA). All data were subjected to analysis of variance (One-way ANOVA) and mean comparisons were carried out by Ducan’s multiple range test (p < 0.05).

To identify all members of the SRS gene family in cotton, the conserved RING-like zinc-finger domain (DUF702) (Pfam ID: PF05142) from the Pfam databases (http://pfam.xfam.org/) were employed as queries to search against three main representative cotton species. Then, the puatative protein sequences using the PlantTFDB (http://planttfdb.cbi.pku.edu.cn/) and SMART (http://smart.embl.de/) databases to confirm the predicted functional domains contained DUF702 (PF05142) families. A total of 53 SRS members were identified in G. arboreum, G. raimondii, and G. hirsutum, of which 14 were GaSRS genes, 13 were GrSRS genes, and 26 were GhSRS genes (13 SRSs from At subgenome, 13 SRSs from Dt subgenome). Therefore, we named them GhSRS1~GhSRS26 based on their gene ID number and genomic distribution. The encoded protein length of GhSRS genes ranging from 195 (GhSRS3) aa to 434 (GhSRS22) aa, and molecular weight (MW) from 22116.24 Da to 45074.19 Da, isoelectric point (pI) varying from 5.49 (GhSRS25) to 9.22 (GhSRS26), in addition to predicting the subcellular location of all GhSRS members, 16 of GhSRS proteins were nucleus-localized. Other basic information for all SRS members in three cotton species were listed in Supplementary Table 3 .

To better understand the phylogenetic and evolution relationships of SRS genes in cotton, the unrooted phylogenetic tree constructed by MEGA X revealed the SRS family genes can be divided into three subgroups ( Figure 1A ). The number of SRS genes in G. hirsutum was almost the sum of the number of those in G. arboreum and G. raimondii, which was consistent with polyploidy and whole-genome duplication (WGD) events during hybridization.

We found that almost all members of this family contained a RING-like zinc-finger domain (CX₂CX₇CX₄CX₂C₂X₆C) through sequence alignment of amino acid residues. But lacking part of the RING domain in GhSRS3 may leading to a decrease in the binding ability to DNA, RNA, protein, and lipid substrates ( Figure 1B ). Moreover, GhSRSs also share a IXGH domain except GhSRS13/26, and this conserved region longer than IXGH domain in GmSRSs ( Figure 1C ) (Zhao et al., 2020).

To investigate the evolutionary relationships of SRS gene family, we constructed a separate rooted phylogenetic tree using 13 plant species genome datasets from lower aquatic to higher terrestrial plants. We totally identified 115 genes in different moss (2 in P. patens), fern (4 in S. moellendorffii), monocotyledons (5 in O. sativa, 10 in Z. mays), and dicotyledons (5 in T. cacao, 10 in A. thaliana, 14 in G. arboreum, 13 in G. raimondii, 26 in G. barbadense and 26 in G. hirsutum), while no SRS gene was found in picophytoplankton (M. pusilla) and algae (Os. tauri, V. carteri) ( Figure 2A, B ). Results revealed that SRS gene first appeared in moss (P. patens), and the number of SRS gene increased dramatically in G. barbadense and G. hirsutum ( Figure 2B ).

The SRS genes in multiple plant species can be divided into 7 clades and named a to g subfamilies. The SRS subfamiliy a and c-gonly existed in monocotyledon and each subfamiliy above contained only one or two members. However, SRS subfamily b was the largest subfamily and contained 108 SRS members including all SRS genes derived from the above dicotyledons and some SRS genes derived from the above monocotyledons ( Figure 2A ). Furthermore, the number of GhSRS genes in G. barbadense or G. hirsutum were almost the sum of SRS genes in G. arboreum and G. raimondii ( Figure 2B ), which confirmed the ideas that allotretraploid cotton (G. barbadense and G. hirsutum) evolved from hybridization and polyploidization between two diploid cotton species (G. arboreum and G. raimondii).

To figure out the structure similarity of the SRS family in cotton, the full-length protein sequence of 26 GhSRSs were aligned to display phylogenetic tree with conserved motifs, exon-intron, and domain structure. Conserved motifs in the SRS protein sequences were performed by MEME online service (https://meme-suite.org/meme/), 12 different motifs were identified and distributed on the GhSRS protein sequences. GhSRSs in the same cluster shared similar conserved motif composition, especially all GhSRS members contained some conserved motifs (motif 1, motif 2, and motif 4) ( Figure 3A ). And the motif number of each protein ranging from five to eleven ( Figure 3A ). Then, we performed the exon-intron gene structure analysis by comparing genomic sequence to the extended cDNA sequence (CDS) of GhSRSs. The number of coding exons of GhSRSs in Gossypium hirsutum was conserved, as they all contained two exons ( Figure 3B ). However, the length of the introns of the GhSRS genes was largely variable, ranging from 73 bp to 497 bp ( Figure 3B ). The domain structure of SRS proteins was analyzed using the SMART protein-domain search interface (http://smart.embl.de/), and the results showed all GhSRSs shared a conserve domain named DUF702, containing the RING-like zinc-finger domain (CX₂CX₇CX₄CX₂C₂X₆C) ( Figure 3C ). Interestingly, only the GhSRS22 possessed a transmembrane domain. The results above indicated that the GhSRS22 was likely to be a membrane-bound protein, which was consistent with the result listed in Supplementary Table 3 .

To investigate the chromosomal distribution and the duplication events of the SRS family in cotton, chromosomal distribution and collinearity analysis was performed. In G. arboreum, 13 SRS family members were distributed on chromosomes At02, At03, At05, At06, At07, At08, At09, At010, At11, At13. Meanwhile, a total of 13 SRSs were unevenly distributed on chromosomes Dt01, Dt02, Dt03, Dt05, Dt06, Dt07, Dt08, Dt09, Dt10, Dt11, Dt13 ( Figure 4 ). Interestingly, although the number of SRS genes in G. hirsutum were not altered, their distributions on chromosomes in G. hirsutum displayed differences compared with those in G. arboreum and G. raimondii, implying structural variation on chromosomes during the process of cotton polyploidization. For instance, GrSRS1 was located on chromosome Dt01 in G. raimondii while there was no SRS gene on chromosome At01 in G. arboreum ( Figure 4 ). However, in G. hirsutum, GhSRS1 was located in chromosome At01 and exhibited high gene collinearity with GrSRS1, probably resulting from the gene duplication and interchromosomal translocation. In addition, one SRS gene loss occurred in chromosome Dt03 ( Figure 4 ). In summary, the SRS genes in G. hirsutum were unevenly distributed on all chromosomes except for the chromosomes At04, At12, Dt04, Dt12. Besides, SRS gene duplication and loss events occurred during ancestral allopolyploidization of G. hirsutum.

To better understand the function of GhSRS genes under salt and drought stresses in Gossypium hirsutum, the expression pattern of GhSRSs in response to salt and drought stress was examined using the FPKM values of GhSRSs extracted from COTTONOMICS (http://cotton.zju.edu.cn/index.htm) ( Figure 5 ). The results showed that most GhSRSs, such as GhSRS1, GhSRS2, GhSRS3, GhSRS4, GhSRS6, GhSRS7, GhSRS10, GhSRS12, GhSRS13, GhSRS14, GhSRS15, GhSRS16, GhSRS17, GhSRS18, GhSRS19, GhSRS20, GhSRS23, GhSRS24, GhSRS25, GhSRS26, expressed at a very low level (average FPKM<5) under control, salt and drought stress ( Figure 5 ). Of the remaining 6 genes, GhSRS21 was induced by salt stress (3h, 12h and 24h after salt treatment) while the transcripts of GhSRS5, GhSRS22 showed high level accumulation after drought treatment ( Figure 5 ), indicating their potential roles in salt or drought tolerance.

To further confirm the nuclear localization of GhSRS21 and presume its potential role in regulation of the expression of eukaryotic genes, a nuclear localization sequence fused with the mCherry (NLS- mCherry) was used as nuclear localization marker and cotransformed with the plasmid (GhSRS21 ORF fused with GFP) to the Arabidopsis protoplasts. There was overlapping between green fluorescence and red fluorescence, indicating GhSRS21 located in the nucleus ( Figure 6 ).

The expression pattern of GhSRS21 under salt stress has been described previously using the COTTONOMICS public database ((http://cotton.zju.edu.cn/index.htm)) ( Figure 5 ). To further confirm the expression pattern of GhSRS21 in response to salt stress, RT-qPCR was performed and the results revealed that the GhSRS21 gene expression was highly induced by salt stress ( Figure 7A ). Subsequently, the function of GhSRS21 under salt stress conditions was determined using the VIGS (virus-induced gene silencing) system. When the albino-like appearance on the leaves of positive control was observed ( Supplementary Figure 1 ), we examined the expression of GhSRS21 in negative control and GhSRS21 VIGS lines. The data showed that the expression levels of GhSRS21 were significantly decreased in GhSRS21 VIGS lines compared with those in negative control ( Figure 7B ). To further understand the role of GhSRS21 in G. hirsutum under salt stress, salt treatment was performed and silenced GhSRS21 in TM-1 resulted in enhanced salt tolerance ( Figure 7C ). In addition, we further detected the activities of antioxidant enzymes and relevant physiological indicators of negative control and VIGS plant under control and salt treatment. The results exhibited the increased activity of catalase (CAT), peroxidase (POD) and increased content of hydrogen peroxide (H₂O₂), malondialdehyde (MDA) significantly in both GhSRS21 VIGS cottons and negative control after salt treatment ( Figures 7D–G ). However, a significantly higher CAT activity was observed in the GhSRS21 VIGS lines compared with that of the negative control after salt treatment ( Figure 7D ). Meanwhile, the content of H₂O₂ and MDA in GhSRS21 VIGS lines was significantly lower than that in the negative control after salt treatment ( Figures 7E, G ). Taken together, GhSRS21 negatively regulates salt tolerance in G. hirsutum through increased antioxidant capacity of cotton.

To explore the potential mechanism of GhSRS21 in regulation of H₂O₂ production, genes differentially expressed in the negative control and GhSRS21 silenced plants were analyzed via the transcriptome data. The results indicated that the number of down-regulated genes (306) was more than that of up-regulated genes (71) in GhSRS21 silenced lines compared with the negative control ( Figure 8A ). Furthermore, the enrichment of the differentially expressed genes (DEGs) in the KEGG pathway was analyzed. Most up-regulated genes were enriched in the pathways of genetic information processing, carbohydrate metabolism, transcription, translation, peroxisome, etc, while the down-regulated genes showed significant enrichment on metabolism, transporters, flavonoid biosynthesis, etc ( Figures 8B, C ). There were four up-regulated genes (GH_A05G0875, GH_A09G1066, GH_A13G1914, GH_D09G1018), enriched in peroxisome KEGG pathway, belong to copper/zinc superoxide dismutase (SODC), and shared potential biological functions in ROS (reactive oxygen species) scavenging. To further figure out their roles in salt resistance of the upland cotton, the expression pattern of the mentioned SODCs above were analyzed using the public COTTONOMICS Database. The results showed that the four SODCs above were significantly suppressed by salt stress while GhSRS21 exhibited the opposite expression pattern under salt treatment ( Figure 8D ). Moreover, the four SODCs were induced when the GhSRS21 was repressed in the VIGS lines ( Figure 8E ). In summary, the GhSRS21 negatively regulates salt tolerance in a manner dependent on reactive oxygen species metabolic process and probably by negative regulation of SODCs expression.