Genome-Wide Identification and Functional Analysis of RF2 Gene Family and the Critical Role of GhRF2-32 in Response to Drought Stress in Cotton

Cotton is an important natural fiber crop. The RF2 gene family is a member of the bZIP transcription factor superfamily, which plays an important role in plant resistance to environmental stresses. In this paper, the RF2 gene family of four cotton species was analyzed genome-wide, and the key gene RF2-32 was cloned for functional verification. A total of 113 RF2 genes were identified in the four cotton species, and the RF2 family was relatively conserved during the evolution of cotton. Chromosome mapping and collinear analysis indicated that fragment replication was the main expansion mode of RF2 gene family during evolution. Cis-element analysis showed that there were many elements related to light response, hormone response and abiotic stress response in the promoters of RF2 genes. The transcriptome and qRT-PCR analysis of RF2 family genes in upland cotton showed that RF2 family genes responded to salt stress and drought stress. GhRF2-32 protein was localized in the cell nucleus. Silencing the GhRF2-32 gene showed less leaf wilting and increased total antioxidant capacity under drought and salt stress, decreased malondialdehyde content and increased drought and salt tolerance. This study revealed the evolutionary and functional diversity of the RF2 gene family, which laid a foundation for the further study of stress-resistant genes in cotton.


Introduction
Plants are affected by a variety of abiotic stresses in the process of growth and development, which leads to a substantial reduction in crop yield [1]. Drought, salt and low temperature are the most common abiotic stresses in nature. About 1/3 of the land in the world is in an arid or semi-arid environment [2], resulting in the degradation of natural vegetation and serious damage to the structure and function of the ecosystem [3]. Saline-alkalized land slows reduced plant growth and development as well as yield and quality by affecting seed germination, root length, plant height and fruit setting rate [4,5]. Under low temperature stress, the related physiological and metabolic activities of plant cells are affected and even lead to plant death [6]. Thus, it is extremely important to study the mechanism of plant response to abiotic stress and to identify the related genes for crop improvement.
Transcription factors, also known as trans-acting elements, are a class of proteins that specifically bind to cis elements in the promoter region of eukaryotic genes, thus specifically activating or inhibiting the expression of related genes at a spatiotemporal manor [7]. The

Identification and Physicochemical Property Analysis of RF2 Gene Family in Cotton
A total of 113 RF2 genes were identified in the four cotton species via homologous search using the Markov Model (PF00170 and PF07716) and specific conserved domain (269851-bZIP_plant_RF2) of RF2 genes (Table S2). A total of 36 and 39 RF2 genes were found in G. hirsutum and G. barbadense, and 20 and 18 RF2 genes were found in G. arboretum and G. raimondii, respectively. According to the position of RF2 gene on chromosome in each cotton species, all identified genes were renamed, such as G. hirsutum (GhRF2-1).
The results of physicochemical properties analysis showed that the amino acid length of cotton RF2 gene protein ranged from 113 (GbRF2-36) to 572 (GbRF2-4), the predicted molecular weight range was 13,220.02 (GbRF2-36)-62,572.12 (GrRF2-4) and the isoelectric point range was 5.38 (GrRF2-8)-9.17 (GaRF2-5). The instability index refers to how stable the protein was in the test tube (≤40, possibly stable; >40, possibly unstable). The prediction results showed that the instability coefficient ranges from  to 73.98 (GhRF2-36). Except for GbRF2-10, GaRF2-12 and GbRF2-36, all RF2 genes may be unstable. Subcellular prediction results showed that all genes were located in the nucleus, but eight genes may also be located in chloroplast (Table S2).

Construction of Phylogenetic Tree of RF2 Genes in Cotton
In order to study the evolutionary relationship of the cotton RF2 gene family, we sequenced the RF2 proteins of four cotton species and constructed a phylogenetic tree ( Figure 1). According to the results, all family members are divided into eight subfamilies, of which subgroup A contains the largest number of genes, namely twenty-four genes, while subgroup G has the least number of genes, only five genes. The results showed that Plants 2023, 12, 2613 3 of 18 each branch contained the genes of four cotton species, and we found that the number of tetraploid cotton in each subgroup was about 2:1 compared with the number of diploid cotton. For example, D subgroup contained two genes from G. hirsutum, two genes from G. barbadense, one gene from G. arboretum and one gene from G. raimondii. This also fully showed the evolutionary relationship among the four cotton species. We also found that the RF2 genes of G. hirsutum and G. barbadense were always clustered together, which may be due to the occurrence of gene repetition events.

Construction of Phylogenetic Tree of RF2 Genes in Cotton
In order to study the evolutionary relationship of the cotton RF2 gene family, we sequenced the RF2 proteins of four cotton species and constructed a phylogenetic tree ( Figure 1). According to the results, all family members are divided into eight subfamilies, of which subgroup A contains the largest number of genes, namely twenty-four genes, while subgroup G has the least number of genes, only five genes. The results showed that each branch contained the genes of four cotton species, and we found that the number of tetraploid cotton in each subgroup was about 2:1 compared with the number of diploid cotton. For example, D subgroup contained two genes from G. hirsutum, two genes from G. barbadense, one gene from G. arboretum and one gene from G. raimondii. This also fully showed the evolutionary relationship among the four cotton species. We also found that the RF2 genes of G. hirsutum and G. barbadense were always clustered together, which may be due to the occurrence of gene repetition events.

Analysis of RF2 Gene Conserved Motif and Exon-Intron Distribution Pattern in Cotton
The information related to the evolution process of gene family can be obtained by analyzing conservative motif. Twelve conserved motif of cotton RF2 proteins were analyzed using the MEME online tool (Figure 2b). The results showed that all cotton RF2 gene encoded proteins had motif 1. All but seven genes contained motif 2. This demonstrated that its structure was highly conservative. The GbRF2-36 gene of subgroup C contained only one motif of motif 1. The GaRF2-12 gene of subfamily F contained two motifs, including motif 1 and motif 2. Other conserved motifs existed only in one or some specific subfamilies, for example, motif 5 existed in A, C, D, and E subfamilies. Motif 12 existed only in the D subfamily. There were different conserved motifs in different subfamilies. In the same subfamily, the distribution of protein conserved motifs was basically the same, indicating that the gene function of the same subfamily was similar, and the homologous relationship was closed, which also proved the accuracy of the constructed phylogenetic tree.

Analysis of RF2 Gene Conserved Motif and Exon-Intron Distribution Pattern in Cotton
The information related to the evolution process of gene family can be obtained by analyzing conservative motif. Twelve conserved motif of cotton RF2 proteins were analyzed using the MEME online tool (Figure 2b). The results showed that all cotton RF2 gene encoded proteins had motif 1. All but seven genes contained motif 2. This demonstrated that its structure was highly conservative. The GbRF2-36 gene of subgroup C contained only one motif of motif 1. The GaRF2-12 gene of subfamily F contained two motifs, including motif 1 and motif 2. Other conserved motifs existed only in one or some specific subfamilies, for example, motif 5 existed in A, C, D, and E subfamilies. Motif 12 existed only in the D subfamily. There were different conserved motifs in different subfamilies. In the same subfamily, the distribution of protein conserved motifs was basically the same, indicating that the gene function of the same subfamily was similar, and the homologous relationship was closed, which also proved the accuracy of the constructed phylogenetic tree.  G subfamilies contained four exons and three introns. This showed that the gene structure from the same subgroup was highly conservative and was closely related to evolution.

Distribution of RF2 Gene Family on Chromosomes in Cotton
We combined the chromosomes of diploid A genome G. arboretum, diploid D genome G. raimondii and tetraploid A, D genome G. hirsutum and G. barbadense and the location information of each RF2 gene on the chromosome to draw the distribution map of RF2 gene on the chromosome (Figure 3). A total of 19 genes were distributed on 10 chromosomes of A genome and 18 genes on 11 chromosomes of D genome. The A subgenome of tetraploid contained 17 RF2 genes, while in D subgenome, G. hirsutum contained 18 RF2 genes and G. barbadense contained 21 RF2 genes. One gene was located on the scaffold in the genomes of G. arboretum, G. hirsutum and G. barbadense.
The gene structure analysis of cotton RF2 showed that the number of exons and introns of cotton RF2 gene varied, with a minimum of two exons and a maximum of six exons (Figure 2c). Most RF2 genes contained four exons and three introns. RF2 genes in the same subfamily had very similar exon and intron structures. For example, all genes of C, D and G subfamilies contained four exons and three introns. This showed that the gene structure from the same subgroup was highly conservative and was closely related to evolution.

Distribution of RF2 Gene Family on Chromosomes in Cotton
We combined the chromosomes of diploid A genome G. arboretum, diploid D genome G. raimondii and tetraploid A, D genome G. hirsutum and G. barbadense and the location information of each RF2 gene on the chromosome to draw the distribution map of RF2 gene on the chromosome (Figure 3). A total of 19 genes were distributed on 10 chromosomes of A genome and 18 genes on 11 chromosomes of D genome. The A subgenome of tetraploid contained 17 RF2 genes, while in D subgenome, G. hirsutum contained 18 RF2 genes and G. barbadense contained 21 RF2 genes. One gene was located on the scaffold in the genomes of G. arboretum, G. hirsutum and G. barbadense.

Intraspecific and Interspecific Collinearity Analysis of RF2 Gene in Cotton
Gene duplication is the main force for promoting the expansion of the gene family. Gene duplication consists of tandem duplication and large fragment duplication [22]. The number of RF2 genes in tetraploid cotton was twice as much as that in diploid cotton, which indicated that the RF2 gene family expanded during cotton polyploidy. In order to study the expansion mode of RF2 gene family in cotton, we analyzed the intraspecific and interspecific collinearity of four cotton species (Figure 4). The results showed that a total

Intraspecific and Interspecific Collinearity Analysis of RF2 Gene in Cotton
Gene duplication is the main force for promoting the expansion of the gene family. Gene duplication consists of tandem duplication and large fragment duplication [22]. The number of RF2 genes in tetraploid cotton was twice as much as that in diploid cotton, which indicated that the RF2 gene family expanded during cotton polyploidy. In order to study the expansion mode of RF2 gene family in cotton, we analyzed the intraspecific and interspecific collinearity of four cotton species (Figure 4). The results showed that a total of 162 duplication gene pairs were obtained in the four cotton varieties, including 16 pairs in G. arboretum, 17 pairs in G. raimondii, 65 pairs in G. hirsutum and 64 pairs in G. barbadense. Of these, only G. barbadense contained one pair of tandem repeat genes. In addition, we comprehensively analyzed the collinearity among different cotton species and found that 556 duplication gene pairs existed in the following six groups of comparisons: Ga-Gr (49), Gb-Ga (90), Gb-Gr (84), Gh-Ga (89), Gh-Gb (157) and Gh-Gr (87). These results show that gene duplication, especially fragment duplication, plays an important role in the expansion of cotton RF2 gene family.
of 162 duplication gene pairs were obtained in the four cotton varieties, including 16 pairs in G. arboretum, 17 pairs in G. raimondii, 65 pairs in G. hirsutum and 64 pairs in G. barbadense. Of these, only G. barbadense contained one pair of tandem repeat genes. In addition, we comprehensively analyzed the collinearity among different cotton species and found that 556 duplication gene pairs existed in the following six groups of comparisons: Ga-Gr (49), Gb-Ga (90), Gb-Gr (84), Gh-Ga (89), Gh-Gb (157) and Gh-Gr (87). These results show that gene duplication, especially fragment duplication, plays an important role in the expansion of cotton RF2 gene family.

Analysis of Cis Elements of RF2 Gene Promoters in Cotton
In order to better understand the function of cotton RF2 gene, we used the Plant-CARE website to analyze the upstream 2000 bp region of the gene coding region and obtain regulatory cis elements. The results showed that the number of light response elements of cotton RF2 gene was the largest ( Figure 5). Certain hormone-related response

Analysis of Cis Elements of RF2 Gene Promoters in Cotton
In order to better understand the function of cotton RF2 gene, we used the PlantCARE website to analyze the upstream 2000 bp region of the gene coding region and obtain regulatory cis elements. The results showed that the number of light response elements of cotton RF2 gene was the largest ( Figure 5). Certain hormone-related response elements were also obtained, including abscisic acid, MeJA, auxin, salicylic acid and gibberellin. A different number of cis elements associated with the defensive stress response were identified. We also found that the RF2 gene may be related to the circadian rhythm and cell cycle. At the same time, there were a certain number of genes associated with abiotic stress, such as drought and low temperature. To sum up, different types and numbers of cis-acting elements were distributed in cotton RF2 promoters. We analyze the function of RF2 gene under abiotic stress in the next step. elements were also obtained, including abscisic acid, MeJA, auxin, salicylic acid and gibberellin. A different number of cis elements associated with the defensive stress response were identified. We also found that the RF2 gene may be related to the circadian rhythm and cell cycle. At the same time, there were a certain number of genes associated with abiotic stress, such as drought and low temperature. To sum up, different types and numbers of cis-acting elements were distributed in cotton RF2 promoters. We analyze the function of RF2 gene under abiotic stress in the next step.

RF2 Gene Expression Patterns in Cotton under Abiotic Stress
In order to explore whether cotton RF2 genes were involved in cotton response to different abiotic stresses, we further analyzed the expression of RF2 genes in upland cotton exposed to drought, salt and low temperature stresses. The results showed that the expression of RF2 genes changed under different abiotic stresses and different treatment times ( Figure 6) (Table S3). Under low temperature stress, 14 candidate genes (GhRF2-2,

RF2 Gene Expression Patterns in Cotton under Abiotic Stress
In order to explore whether cotton RF2 genes were involved in cotton response to different abiotic stresses, we further analyzed the expression of RF2 genes in upland cotton exposed to drought, salt and low temperature stresses. The results showed that the expression of RF2 genes changed under different abiotic stresses and different treatment times ( Figure 6) (Table S3). Under low temperature stress, 14 candidate genes (GhRF2-2

Analysis of the Subcellular Localization of GhRF2-32
In order to determine the subcellular localization of the GhRF2-32 gene, the fusion expression vector pCAMBIA2300-DsRED2-GhRF2-32 was constructed and transformed into GV1301 for infecting tobacco plants. As shown in Figure 8, pCAMBIA2300-DsRED2 (empty vector) displayed strong fluorescent signals in the cell membrane and nucleus of tobacco epidermal cells, but pCAMBIA2300-DsRED2-GhRF2-32 was detected only in the nucleus, which confirmed that GhRF2-32 was expressed in the nucleus.
In order to determine the subcellular localization of the GhRF2-32 gene, the fusion expression vector pCAMBIA2300-DsRED2-GhRF2-32 was constructed and transformed into GV1301 for infecting tobacco plants. As shown in Figure 8, pCAMBIA2300-DsRED2 (empty vector) displayed strong fluorescent signals in the cell membrane and nucleus of tobacco epidermal cells, but pCAMBIA2300-DsRED2-GhRF2-32 was detected only in the nucleus, which confirmed that GhRF2-32 was expressed in the nucleus.

Virus-Induced Gene Silencing (VIGS) of GhRF2-32 in Upland Cotton
Virus-induced gene silencing (VIGS) was used to study GhRF2-32 gene function in upland cotton. In this study, the pds gene served as positive control, which encodes a gene associated with chlorophyll synthesis. Silencing the pds gene, cotton leaves could not synthesize chlorophyll and the new leaves undergo albinism at a later stage. The occurrence of an albino phenotype after infection with cotton plants reveals successful gene silencing ( Figure 9A). Subsequently, randomly selected cotton leaves were subjected to fluorescence quantification, and the expression levels of TRV empty gene and GhRF2-32 gene were verified via qRT-PCR. The GhRF2-32 gene expression levels of the normal and control plants did not significantly change, while the expression of target genes in cotton plants significantly decreased after silencing, indicating that GhRF2-32 was successfully silenced ( Figure 9A). An 18% PEG6000 solution was used to treat cotton plants with genesilencing and empty vectors to simulate drought stress. Phenotypes were observed 72 h after drought stress, and the control cotton plants withered more severely than the genesilenced plants ( Figure 9B).

Virus-Induced Gene Silencing (VIGS) of GhRF2-32 in Upland Cotton
Virus-induced gene silencing (VIGS) was used to study GhRF2-32 gene function in upland cotton. In this study, the pds gene served as positive control, which encodes a gene associated with chlorophyll synthesis. Silencing the pds gene, cotton leaves could not synthesize chlorophyll and the new leaves undergo albinism at a later stage. The occurrence of an albino phenotype after infection with cotton plants reveals successful gene silencing ( Figure 9A). Subsequently, randomly selected cotton leaves were subjected to fluorescence quantification, and the expression levels of TRV empty gene and GhRF2-32 gene were verified via qRT-PCR. The GhRF2-32 gene expression levels of the normal and control plants did not significantly change, while the expression of target genes in cotton plants significantly decreased after silencing, indicating that GhRF2-32 was successfully silenced ( Figure 9A). An 18% PEG6000 solution was used to treat cotton plants with gene-silencing and empty vectors to simulate drought stress. Phenotypes were observed 72 h after drought stress, and the control cotton plants withered more severely than the gene-silenced plants ( Figure 9B).  To study the potential mechanism of GhRF2-32-mediated drought tolerance, the T-AOC and MDA contents were measured in both normal control plants and VIGS-silenced plants ( Figure 9C,D). For T-AOC, there was no significant difference between CK and TRV2:00 at 72 h after GhRF2-32 gene silencing. However, in GhRF2-32-silenced seedlings, there was a significant difference between 12 and 24 h between the CK and TRV2:00-silenced plants, while there was no significant difference between 4 and 72 h between the CK and TRV2:00 plants. The results of MDA content were different from those of T-AOC. GhRF2-32 gene showed significant differences from CK and TRV after 72 h of VIGS silencing. The MDA content in VIGS silenced plants was significantly lower than the controls. Based on the above results, it can be inferred that after the GhRF2-32 gene was silenced, T-AOC content was increased and MDA content was decreased. This suggests that GhRF2-32 functions in abiotic stress possibly through regulating active oxygen scavenging system.

Identification and Physicochemical Property Analysis of RF2 Gene Family in Cotton
The genomic sequences of Gossypium hirsutum, Gossypium barbadense [23], Gossypium arboretum [24] and Gossypium raimondii [25] were obtained from CottonMD [26]. Arabidopsis bZIP protein sequences were downloaded from Tair (https://www.arabidopsis.org/ (accessed on 12 August 2022)) and used as reference sequences to blast cotton genome To study the potential mechanism of GhRF2-32-mediated drought tolerance, the T-AOC and MDA contents were measured in both normal control plants and VIGS-silenced plants ( Figure 9C,D). For T-AOC, there was no significant difference between CK and TRV2:00 at 72 h after GhRF2-32 gene silencing. However, in GhRF2-32-silenced seedlings, there was a significant difference between 12 and 24 h between the CK and TRV2:00-silenced plants, while there was no significant difference between 4 and 72 h between the CK and TRV2:00 plants. The results of MDA content were different from those of T-AOC. GhRF2-32 gene showed significant differences from CK and TRV after 72 h of VIGS silencing. The MDA content in VIGS silenced plants was significantly lower than the controls. Based on the above results, it can be inferred that after the GhRF2-32 gene was silenced, T-AOC content was increased and MDA content was decreased. This suggests that GhRF2-32 functions in abiotic stress possibly through regulating active oxygen scavenging system.

Identification and Physicochemical Property Analysis of RF2 Gene Family in Cotton
The genomic sequences of Gossypium hirsutum, Gossypium barbadense [23], Gossypium arboretum [24] and Gossypium raimondii [25] were obtained from CottonMD [26]. Arabidopsis bZIP protein sequences were downloaded from Tair (https://www.arabidopsis.org/ (accessed on 12 August 2022)) and used as reference sequences to blast cotton genome using TBTools v1.106 [27]. Markov Model (HMM) PF00170 and PF07716 of the bZIP family were obtained from the Pfam (https://pfam.xfam.org/ (accessed on 12 August 2022)) [28]. The preliminary identification of the cotton bZIP family was conducted using the HMM Search of TBTools. The candidate proteins obtained by the two methods were merged and deduplicated. Then, NCBI Batch CD-Search was used to screen whether the candidate protein has the specific conserved domain (269851-bZIP_plant_RF2) of the RF2 family. The location information of the screened genes was extracted and further duplicated by the location information, and the final identification result was obtained. The sequence information of RF2 protein was obtained via TBTools, and the number of amino acids, relative molecular weight, isoelectric point and instability index were predicted using the ProtParam tool (https://web.expasy.org/protparam/ (accessed on 12 August 2022)) [29]. WoLF PSORT (https://wolfpsort.hgc.jp/ (accessed on 12 August 2022)) was used to predict the subcellular localization of the protein sequences of family members [30].

Phylogenetic Analysis of RF2 Gene in Cotton
Multiple sequence alignment of all cotton RF2 proteins was carried out via Clustal X, and the phylogenetic tree was constructed via the neighbor-joining method using MEGA-X [31]. The bootstrap was set to 1000 and model was selected as the p-distance. Finally, the evolution tree was modified and visualized by an online tool EvolView [32].

Analysis of Conserved Motifs and Gene Structures of RF2 Genes in Cotton
According to the location information of the extracted family members, the exonintron structure of each gene was analyzed using TBTools, and the RF2 gene structure map was drawn. The conserved motif of RF2 protein sequence was identified via the online tool MEME (http://meme-suite.org/ (accessed on 12 August 2022)) [33]. The number of motifs was set to 12, the minimum length of motif was set to 10, the maximum length was set to 150, and other parameters were set by default. The output file in XML format was downloaded and visualized through TBTools.

Mapping RF2 Genes in Cotton Chromosome
According to the genome and genome annotation information, TBTools was used to extract the position of the RF2 genes on the chromosome and the software was used to visualize the results.

Intraspecific and Interspecific Collinearity Analysis of RF2 Genes in Cotton
The collinearity analysis of RF2 gene family sequences was performed using TBTools software, including fragment repeat, tandem repeat and genome-wide duplication events. Intraspecific and interspecific BLAST was performed on the diploid A/D genome (G. arboretum/G. raimondii) and tetraploid A/D subgenome (G. hirsutum/G. barbadense), respectively, and then the intraspecific duplicated gene pairs of G. arboretum, G. raimondii, G. hirsutum and G. barbadense and the interspecific duplicated gene pairs of G. hirsutum-G. arboretum, G. hirsutum-G. raimondii, G. barbadense-G. arboretum, G. barbadense-G. raimondii and G. hirsutum-G. barbadense obtained from duplication event analysis were used to map the intraspecific and interspecific covariance circles.

Prediction of Cis Elements of RF2 Gene in Cotton
According to the location information of RF2 genes, the 2000 bp upstream DNA sequences of each RF2 gene were extracted from genome sequence files. The cis elements of the gene promoter region were predicted via PlantCARE (http://bioinformatics.psb. ugent.be/webtools/plantcare/html/ (accessed on 12 August 2022)) [34]. The predicted regulatory elements were visualized using TBTools.

Expression Pattern of RF2 under Abiotic Stress
The transcriptome data of G. hirsutum (PRJNA490626) under drought, salt and low temperature stress were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/ (accessed on 12 August 2022)). First of all, the data were filtered using Trimmomatic software [35], the G. hirsutum reference genome database was built via the hisat2 tool [36], and then the transcriptome data were compared to the reference genome of the database. The FPKM values (fragments per kilobase of transcript per million fragments) of all genes contained in the genome were calculated by Cufflinks 2.0.2 software [37]. According to the unique ID combined transcriptome data of the proteome file corresponding to upland cotton RF2 gene, the expression amount of each gene was obtained, the results were visualized and the expression heat map was drawn using TBTools.

Stress Treatments and qRT-PCR Analysis
G. hirsutum acc. TM-1 was grown in the experimental field at the Anyang Institute of Technology. In order to study the effects of drought, salt and low temperature stress on RF2 gene expression in G. hirsutum, when the cotton seedlings grew to three true leaf stage, the seedlings were treated with 18% PEG, 250 mM NaCl and in a 4 • C incubator, respectively. The leaves of various stress treatments and the control group were collected at 0 h, 1 h, 3 h, 6 h, 12 h and 24 h, respectively. The leaves were immediately frozen in liquid nitrogen and stored at −80 • C. All treatments were repeated three times.
Total RNAs were extracted from each sample by using the EASYspin Plus Plant RNA Kit (RN38, Aidlab Biotech, Beijing, China). The quality of RNAs was determined via agarose gel electrophoresis and a Nanodrop2000 nucleic acid analyzer. mRNAs were reversely transcribed to cDNAs by using TranScript All-in-One First-Strand cDNA Synthesis SuperMix for qPCR (Transgen Biotech, Beijing, China). Finally, the qRT-PCR experiment was carried out on the ABI 7500 Fast Real-Time PCR system (Applied Biosystems, Waltham, MA, USA) using the TransStart Top Green qPCR SuperMix kit (Transgen Biotech, Beijing, China). Actin was saved as reference gene identified by RefFinder [38]. The primers of internal reference genes and differentially expressed genes in G. hirsutum were designed by the Prime-Blast of NCBI and listed in Table S1. Three biological replicates were run for each treatment and the control. Three repeats were performed on each cDNA sample, and the relative expression of the genes was calculated according to the CT value via the 2 −∆∆Ct method.

Subcellular Localization of GhRF2-32
According to the plasmid sequence of pCAMBIA2300-DsRED2 and the target gene GhRF2-32 with the stop codon removed, primers were designed according to the specification, and BamHI was selected as the tangent point of line plasmid. The specific primers used were shown in Table S1. The PCR amplification product was ligated into the pCAMBIA2300-DsRED2 vector using a ClonExpress ® II One Step Cloning Kit, according to the manufacturer's instructions. The recombinant plasmid was subsequently transformed into DH5α competent cells. The recombinant results were verified via PCR amplification and sequencing. Positive clones were screened, and plasmids were extracted to obtain the pCAMBIA2300-DsRED2-GhRF2-32 vector construct. The plasmid was injected into leaves for 4 weeks under normal growth conditions. Two days later, the leaves injected with plasmid were removed and the red fluorescence of the DsRED2 fusion protein was observed using the Leica DMi8 confocal laser scanning microscope (Leica, Wetzlar, Germany).

Using Virus-Induced Gene Silencing (VIGS) to Test RF2 Gene Function
A 300 bp fragment of RF2-32 gene was cloned into the TRV-RNA vector to develop the TRV2:GhRF2-32 using the Kpn1 and Xba1 restriction sites. The primers were designed by SnapGene 4.1.9 software (Table S1). Similarly, TRV2:PDS was constructed as a visual marker to monitor silencing efficiency. As a negative control, TRV2:00 (empty vector) was used. All the vectors were inserted into the GV3101. We injected the cotyledons of 10 days old G. hirsutum seedlings at 25 • C. Th cotton plants with injected TRV2:PDS show an albino phenotype after 2 weeks of infiltration. When the cotton reaches the three-leaf stage, some of the plants were watered with water as the control and the others were treated with 18% PEG6000 solutions until the phenotypes appeared.

Determination of Enzyme Activity in Gene-Silencing Cotton after Drought Stress Treatment
The MDA content was measured using the Malondialdehyde (MDA) assay kit (BC0025, Solarbio, Beijing, China). The total antioxidant capacity (T-AOC) assay kit (A015-3-1, Nanjing Jiancheng Bioengineering Institute, Nanjing, China) was used to determine the T-AOC content. The experimental procedure was carried out in accordance with the instructions. During the measurement process, a spectrophotometer was used to determine the content of about 0.1 g of leaf samples from the wild type control and the GhRF2-32-silenced plants after 72 h of drought treatment. The experiment was repeated three independent times for three biological repetitions.

Discussion
Abiotic stress is a serious problem for plant growth and development, which threatens the growth and development of plants from all aspects, so it is urgent to study the response genes under stress. As an important economic crop in the world, cotton is also widely affected by almost all abiotic stresses, including drought, salt and high/low temperatures. RF2 gene is an important gene belonging to the bZIP transcription factor superfamily. It has been identified as being associated with disease and abiotic stress in rice and wheat, but it has not been reported in cotton.
In this study, we identified a total of 113 RF2 genes in four cotton species, including G. hirsutum, G. barbadense, G. arboretum and G. raimondii containing 36, 39, 20 and 18 RF2 genes, respectively. Our results show that the number of RF2 in tetraploid cotton is almost twice as much as that in diploid cotton, which confirms that allotetraploid cotton is produced via the natural hybridization of diploid cotton seeds containing A genome and D genome [39,40]. In the past, most of the results of family identification focused on the identification of bZIP superfamily, such as A. thaliana [41] and tobacco [42]. The physical and chemical properties of cotton RF2 protein showed that its sequence length, relative molecular weight and isoelectric point distribution were very wide. This may be related to the large-scale duplication of the genome of tetraploid cotton and the existence of a large number of RF2 genes in tetraploid cotton. In this study, we predicted the subcellular location of cotton RF2 and found that almost all genes were located on the nucleus, which was consistent with the results of previous studies. RF2b was clearly located on the nucleus in rice [15].
We constructed a phylogenetic tree of 113 RF2 proteins and divided them into eight subfamilies. Each subfamily contains four RF2 genes of cotton species. At the same time, we found that the genes in the same subgroup had similar motif distribution patterns. The results of gene structure analysis also showed that RF2 was highly conserved. The RF2 of the same subfamily has similar exon/intron number and exon length.
In our study, we found that the distribution of RF2 genes on each chromosome was uneven, but we found that most of the number and location of genes on each chromosome of the two tetraploid cotton species were one to one, this phenomena was also observed for other transcription factors in other plant species [43][44][45][46][47][48]. This may be due to the fact that allotetraploids share the common ancestors A0 (the possible extinct common ancestor of the existing A genome) and D5 (G. raimondii) [49].
Cis element analysis showed that RF2 gene was associated with abiotic stress responses such as salt, drought and low temperature. Heat map analysis revealed that all RF2 genes were expressed differently in upland cotton under different treatment times. Under salt stress, the growth of the A. thaliana lines overexpressing the RF2 gene was significantly worse than that of the wild type, indicating that TaRF2 (TabZIP3) gene played a positive regulatory role in plant salt stress response [17]. There are also many bZIP genes associated with abiotic stress. OsbZIP71 specifically binds to the promoter of OsNHX1 gene to enhance the salt stress resistance of rice by regulating cell osmosis [50]. In A. thaliana, AtbZIP17 binds specifically to TGACG elements to regulate downstream salt stress response genes. The transcriptional activation state of TabZIP2 increased in ploidy under severe drought conditions [51]. ZmbZIP72 improved drought resistance in transgenic maize plants [52]. The overexpression of TabZIP6 reduced the freezing resistance of transgenic A. thaliana seedlings [53]. bZIP6 gene improved the freezing resistance of overexpressed plants in A. thaliana [54].
Drought is one of the important abiotic factors that restrict cotton yield and fiber quality. Traditional breeding methods are difficult to meet the needs of contemporary production, so it is important to identify and explore drought resistance-related genes in cotton. RF2-32 gene was isolated from upland cotton TM-1. During the VIGS process of GhRF2-32, the emergence of an albino phenotype in plants indicates successful gene silencing. In addition, the expression of GhRF2-32 was significantly reduced after silencing. After drought stress treatment, compared with the wild group plants and the control group plants, the silenced cotton plants showed a lighter wilt phenomenon. Gene silencing indicates that GhRF2-32 has a negative regulatory effect on cotton drought resistance. There are various antioxidant macromolecules, small molecules and enzymes in plants. Their role is to eliminate various reactive oxygen species produced by plants under stress, thereby preventing the production of oxidative stress induced by reactive oxygen species. The total level of these antioxidant macromolecules, small molecules and enzymes reflects the total antioxidant capacity within the system. MDA is the final product produced by the decomposition of oxygen free radicals after acting on unsaturated fatty acids in lipids, and this reflects the level of lipid oxidation in plants. Drought stress can cause the accumulation of a large amount of reactive oxygen species, cause cell membrane lipid peroxidation and cause plant cell damage. Therefore, we evaluated the resistance of plants to environmental stress by measuring their total antioxidant capacity and MDA content under drought stress. The analysis after gene silencing showed that the total antioxidant capacity of cotton plants was significantly higher than that of the control cotton plants, and the content of malondialdehyde was significantly lower than that of non-silenced cotton plants, indicating that the metabolism in cotton leaves was slow, the degree of severity increased, which was consistent with the phenotype after drought treatment. It can be inferred that the silencing of the GhRF2-32 gene improves the ability of cotton plants to eliminate the damage caused by reactive oxygen species, making them less susceptible to drought stress. Through our research on the function of GhRF2-32 gene, including qRT-PCR, gene silencing techniques in cotton and the detection of physiological indicators after silencing, it can be inferred that GhRF2-32 gene is involved in the regulation of cotton drought resistance.

Conclusions
Advanced genome sequencing and bioinformatics analysis have been bringing the related research into the next era of genome wide identification and gene function analysis [55]. In this study, a total of 113 RF2 genes were identified in four different cotton species. The RF2 family was analyzed at the evolutionary level by physical and chemical property analysis, phylogenetic tree construction, conservative motif analysis, gene structure analysis and collinearity analysis. The results showed that RF2 genes were conserved in the process of evolution and expanded by gene duplication in different ploidy cotton species. Twelve key candidate genes were screened by heat map and qRT-PCR analysis, and GhRF2-32 of the most significantly expressed core genes was identified for follow-up research. In the VIGS experiment, wild-type and control plants showed significant wilt compared to gene silenced cotton plants under the same drought treatment. The T-AOC of gene silenced cotton plants is higher than that of wild type and control group, and the MDA content is lower than that of wild type and control group, which indirectly indicates that GhRF2-32 gene is related to drought resistance. Although the further functional analysis of these genes using the more advanced technologies may be required, such as CRISPR/Cas genome editing [56,57], this study lays a foundation for the mining of stress-resistant genes in cotton.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants12142613/s1, Table S1: Use NCBI designed differentially expressed gene-specific primer sequences; Table S2: Physico-chemical and biochemical characteristics of RF2 genes in cotton; Table S3: RNA-Seq data of RF2 gene family under different abiotic stress. Funding: This work was supported by the National Natural Science Foundation of China (31471548 to R.P.), the Science and Technology Research Project of Anyang City (2022C01NY001 to R.P.), the Workstation of Central Plains Scholars (ZYGZZ2021050 to R.P.), the Postgraduate Education Reform Research and Practice Project of Henan Province (2021SJGLX067Y to R.P.), and the Postgraduate Education Reform and Quality Improvement Project of Henan Province (YJS2022JD47 to R.P.). This work was also partially supported by Cotton Incorporated and the National Science Foundation (1658709) to B.Z.