Evolutionary and Characteristic Analysis of RING-DUF1117 E3 Ubiquitin Ligase Genes in Gossypium Discerning the Role of GhRDUF4D in Verticillium dahliae Resistance

Verticillium wilt, primarily induced by the soil-borne fungus Verticillium dahliae, is a serious threat to cotton fiber production. There are a large number of really interesting new gene (RING) domain-containing E3 ubiquitin ligases in Arabidopsis, of which three (At2g39720 (AtRHC2A), At3g46620 (AtRDUF1), and At5g59550 (AtRDUF2)) have a domain of unknown function (DUF) 1117 domain in their C-terminal regions. This study aimed to detect and characterize the RDUF members in cotton, to gain an insight into their roles in cotton’s adaptation to environmental stressors. In this study, a total of 6, 7, 14, and 14 RDUF (RING-DUF1117) genes were detected in Gossypium arboretum, G. raimondii, G. hirsutum, and G. barbadense, respectively. These RDUF genes were classified into three groups. The genes in each group were highly conserved based on gene structure and domain analysis. Gene duplication analysis revealed that segmental duplication occurred during cotton evolution. Expression analysis revealed that the GhRDUF genes were widely expressed during cotton growth and under abiotic stresses. Many cis-elements related to hormone response and environment stressors were identified in GhRDUF promoters. The predicted target miRNAs and transcription factors implied that GhRDUFs might be regulated by gra-miR482c, as well as by transcription factors, including MYB, C2H2, and Dof. The GhRDUF genes responded to cold, drought, and salt stress and were sensitive to jasmonic acid, salicylic acid, and ethylene signals. Meanwhile, GhRDUF4D expression levels were enhanced after V. dahliae infection. Subsequently, GhRDUF4D was verified by overexpression in Arabidopsis and virus-induced gene silencing treatment in upland cotton. We observed that V. dahliae resistance was significantly enhanced in transgenic Arabidopsis, and weakened in GhRDUF4D silenced plants. This study conducted a comprehensive analysis of the RDUF genes in Gossypium, hereby providing basic information for further functional studies.


Introduction
Cotton (Gossypium) is well-known as the most important natural fiber. The cultivated tetraploids upland cotton (Gossypium hirsutum, AD 1 ) and sea island cotton (G. barbadense, AD 2 ) have resulted from hybridization and genome doubling from the diploid ancestors of G. herbaceum (A 1 ) or G. arboreum (A 2 ) and G. raimondii (D 5 ) [1][2][3]. More than 90% of cultivated cotton crops comprise upland cotton, owing to its high fiber production volumes and high tolerance to various environmental conditions. Verticillium wilt is primarily caused by the soil-borne fungus V. dahliae, and is a serious disease that poses a significant threat to the fiber industry. Because of cotton's limited germplasm diversity, V. dahliaeresistant cotton varieties are difficult to obtain via traditional breeding practices [4]. Efforts CRI_v1.0) [37] and G. raimondii (D 5 , JGI_v2.1) [38] were downloaded from the CottonGen website [39]. The Basic Local Alignment Search Tool (BLAST) program BLASTp was used to detect candidate RDUF genes using the Arabidopsis RDUF amino acid (aa) sequences of AtRDUF1 (AT3G46620), AtRDUF2 (AT5G59550), and AtRHC2A (AT2G39720) as queries. Then, the Pfam RDUF1117 domain (PF06547) was used for the examination of RDUF members via a HMMER search. The RDUF genes were named based on their chromosomal location. Genes in tetraploid cotton were named as per the homologous relationship of each subgenome, with "A" and "D" indicating the homologous genes in At and Dt subgenomes, respectively. In addition, the theoretical molecular weight (MW) and isoelectric point (pI) were estimated by ExPASy [40], and the subcellular localization was evaluated by TargetP-2.0 Server (available online: http://www.cbs.dtu.dk/services/TargetP/ accessed on 30 December 2020). To detect cis-elements in the RDUF promoter regions, the 2 kb upstream regions of the start codon of the RDUF genes were submitted to the PlantCARE database [41].

Phylogenetic, Gene Structure and Conserved Domain, and Motif Analysis
A maximum likelihood (ML) phylogenetic tree was built using the bootstrap test method with 1000 replicates [42], and then the tree was visualized via the Interactive Tree of Life (iTOL) online tool [43]. Subsequently, gene structure and the conserved domain of the RDUF proteins were detected and displayed using the Gene Structure Display Server (GSDS, v2.0) [44], SMART [45], and TBtools (v1.0983) [46] softwares. Moreover, the conserved motifs in the RDUF proteins were identified using the MEME program [47], and the functional annotations were confirmed using the InterProScan website.

Chromosomal Location, Gene Synteny, and RNA-Seq Dataset Analysis
The chromosomal localization of the RDUF genes were determined and displayed using the Mapchart 2.2 software [48], based on the downloaded genomic datasets. MC-ScanX [49] software was used for gene synteny analysis, whereas Circos (v0.69) [50] software was used for graphical depiction. The expression profiles of GhRDUFs were analyzed from publicly released RNA-Seq datasets in silico, with the BioProject accession number PRJNA490626 [2]. Transcript levels were calculated with HISAT and StringTie softwares [51], in the form of fragments per kilobase million (FPKM) values.

Transcription Factors (TFs) and miRNAs in Targeting GhRDUF Genes
TFs involved in the regulation of the GhRDUF genes were predicted using PlantRegMap [52]. The cDNA of the GhRDUF homologs were submitted to the psRNATarget website to search for potential miRNAs [53]. Then, Cytoscape (v3.7.2) was used to display the regulatory relationships between the predicted TFs, miRNAs, and the targeted GhRDUF genes [54].

Cotton Seedlings under Hormone Treatments
G. hirsutum cv. Zhongzhimian No. 2 were stored in our laboratory and used for gene expression analysis in response to plant hormone treatments. The cotton seedlings were cultivated in Hoagland liquid medium. The seedlings were sprayed uniformly at the three-leaf stage with 2 mM methyl salicylic acid (MeSA), 100 µM methyl jasmonate (MeJA), or ethylene (ET) released from 5 mM ethephon. The leaves from three cotton seedlings were collected at each time points of 0, 3, 6, 12, 24, and 48 h after spraying with MeSA, MeJA, and ET, respectively, and then immediately frozen in liquid nitrogen for total RNA extraction.

GFP Vector Construction and Fluorescent Visualization
The full-length cDNA of GhRDUF4D without the termination codon was cloned into the pCAMBIA2300-GFP vector, via the ClonExpress ® MultiS One Step Cloning Kit (Vazyme, Nanjing, China). The precursor sequence (gra-miR482c) and its mutant (gra-miR482c-mut) Biomolecules 2021, 11, 1145 4 of 24 was biosynthesized by BGI Biotechnology Co., Ltd. (Wuhan, China), and then inserted into the pBI121 vector. After transformation into the A. tumefaciens strain GV3101, the transient expression assay was performed on tobacco leaves using an MMA solution (10 mM MgCl 2 , 10 mM N-morpholino ethanesulfonic acid and 200 mM acetosyringone), to validate the interactional relationships between gra-miR482c and GhRDUF4D using an in vivo plant imaging system (NightSHADE LB 985, Berthold, Germany).

Overexpression Vector Construction and Screening of Transgenic Arabidopsis
The full cDNA length of GhRDUF4D was cloned to confirm its sequence. The ClonExpress ® MultiS One Step Cloning Kit (Vazyme, Nanjing, China) was used to construct GhRDUF4D to the pCAMBIA2300 vector under the control of a CaMV35S promoter, the recombinant vector was named as 35S::GhRDUF4D in the present study. The 35S::GhRDUF4D construction was introduced into the A. tumefaciens strain GV3101, and then transformed into A. thaliana using the floral-dip method. Positive transgenic Arabidopsis were confirmed using the reverse transcription polymerase chain reaction (PCR) methods. The specific primers used in this study are provided in Supplementary Table S1. The T 0 -T 3 seeds were screened in MS medium containing antibiotics. The T 3 transgenic Arabidopsis lines with the correct segregation ratio were used for the V. dahliae resistance analysis.

VIGS Vector Construction and Assays Performing
Based on the expression module, GhRDUF4D was selected for functional analysis using the virus-induced gene silencing (VIGS) approach. In brief, the 300 bp fragment of GhRDUF4D was cloned from the root tissue of G. hirsutum cv. Zhongzhimian No.

Then, ClonExpress ® MultiS One
Step Cloning Kits (Vazyme, Nanjing, China) were used to link the tobacco rattle virus (TRV) vector (pYL156) with the Xba I and BamH I restriction sites, the recombined vector was named as TRV::GhRDUF4D in this study. Next, the recombined vector TRV::GhRDUF4D, the pYL156 empty vector TRV::00, the marker vector TRV::CLA1 (Cloroplastos alterados 1), and the auxiliary vector pLY192 were introduced into the A. tumefaciens strain GV3101. Subsequently, a single colony of the strain was resuspended in an MMA solution to an OD600 of 1.0. Finally, the suspension of TRV::00 and TRV::GhRDUF4D were infected onto the cotyledons of the cotton seedlings. When the photobleaching phenotype was presented in the marker construct of TRV::CLA1, the true leaves were collected from each treatment to detect GhRDUF4D expression using RT-qPCR. Subsequently, the GhRDUF4D silenced plants and the control were subjected to V. dahliae inoculation.

V. dahliae Infection and Disease Evaluation
The highly aggressive, defoliating V. dahliae strain, Vd991, was used for inoculation. In brief, the spore suspensions of Vd991 were prepared at concentrations of 1 × 10 7 spore mL −1 with sterile water. The cotton and Arabidopsis plant roots were then incubated with the conidial suspension for 10 min. The roots of Arabidopsis seedlings planted on MS medium were treated with 2 µL of 5 × 10 3 spore mL −1 conidial suspension for infection. The disease index test was calculated as previously described [5]. Cotton stems were cut from each line at the same position, using a microscope to determine vascular wilt symptoms. The fungal biomass was detected according to the method described by Liu et al. [55]. In brief, total DNA was extracted from cotton seedling stems at 3 weeks post-inoculation, and then the specific primers were employed to detect fungal DNA using quantitative PCR (Table S1). The V. dahliae recovery assay was performed as per Zhang et al. [11]. The significant differences for disease index and fungal biomass analysis were determined using the Student's t test.

Reactive Oxygen Species (ROS) Analysis and Callose
ROS was detected using 3,3 -diaminobenzidine (DAB) staining. Leaves from TRV::0 and TRV::GhRDUF4D cotton plants from 12 h after inoculation were suspended in DAB Biomolecules 2021, 11, 1145 5 of 24 staining solution (1 mg/mL, pH = 7.5) in darkness for 8 h. Then, leaves were decolorized in 95% ethanol and absolute ethanol in boiling water, until the green color faded. Next, 70% glycerin was added to soak the leaves; they were then observed under a stereomicroscope. At least three plants were sampled for each treatment.
Callose depositions were observed using aniline blue staining. Leaf samples at 7 days post-inoculation (dpi) were first destained in a fixative solution (3:1 ethanol/acetic acid) for 3 h, and then soaked in 70% and 50% ethanol for 2 h. Afterwards, the leaves were transferred into water for 12 h and destained in NaOH solution (10%, w/v) for 2 h. Finally, the leaves were stained with 0.01% (w/v) aniline blue for 3 h, and visualized using fluorescence microscopy.

RT-qPCR Analysis
Total RNA was isolated using RNA extraction kits (Vazyme, Nanjing, China). cDNA was reverse-transcribed using the PrimeScript™ RT reagent kit (TaKaRa, Dalian, China), as per the manufacturer's instructions. Specific primers were designed via the Primer6 software and are listed in Table S1. RT-qPCR analysis was performed as per Zhao et al. [56]. The housekeeping genes, GhUBQ7 in cotton and AtSAND in Arabidopsis, were used as internal references. The relative expression levels were calculated using the 2 −∆∆CT method. The assays were performed in triplicate. The significant differences were determined using the Student's t test.

Identification and Phylogenetic Analysis of RDUF Family Genes
The Pfam family of DUF1117 (PF06547) was used as a query to search for RDUF. At the same time, the AtRDUF amino acid (aa) sequences of AtRDUF1 (AT3G46620), AtRDUF2 (AT5G59550), and AtRHC2A (AT2G39720) in Arabidopsis were used for the BLASTp search in the Gossypium protein database. We detected 6, 7, 14, and 14 RDUF genes in the A 2 , D 5 , AD 1 , and AD 2 genomes, respectively. We observed that the RDUF gene number in the tetraploid G. hirsutum and G. barbadense was probably twice than that observed in the diploid cottons (G. arboreum and G. raimondii), suggesting hybridization and whole genome duplication (WGD) events. Meanwhile, 3 TcRDUF genes in cacao (Theobroma cacao), 10 GmRDUF genes in soybean (Glycine max), 12 BnRDUF genes in oilseed (Brassica napus), and 7 HaRDUF genes in sunflower (Helianthus annuus) were also detected. The phylogenetic tree showed that the 76 RDUF proteins were categorized into 3 groups: Group I, Group II, and Group III ( Figure 1). In addition, we also found that the RDUF family proteins were clustered closely in the different species, indicating that RDUF proteins evolved independently after separation from other species. Another phylogenetic tree containing monocotyledonous crop RDUFs was built, including three OsRDUF proteins from Oryza sativa, five ZmRDUF proteins from Zea mays, and eight TaRDUF proteins from Triticum aestivum (Table S2). The RDUFs in monocotyledonous crops were clustered far from those found in dicotyledonous crops ( Figure S1, in green clade), suggesting that the RDUFs in monocotyledonous and dicotyledonous crops likely evolved independently. The RDUFprotein phylogenetic trees from the four cotton species were also analyzed. We found that the RDUF members in all four species were categorized into three groups ( Figure S2).
Data from the 41 RDUFs in Gossypium were also investigated. The coding sequence (CDS) length of the 41 RDUF genes ranged from 953 bp (GrRDUF4) to 1312 bp (GbRDUF7A). Most of the RDUF genes did not contain an intron, but GaRDUF3, GrRDUF4, GbRDUF2D, and GbRDUF7A had two exons each. The amino acid sequences of the RDUF proteins were also characterized. The MW ranged from 34.13 kDa (GbRDUF3A) to 46.92 kDa (GbRDUF7A), and the pI ranged from 5.17 (GbRDUF4A) to 9.27 (GaRDUF4). Subcellular location analysis showed that 23 of the 41 RDUF proteins were predicted to be located in the nucleus, and 13 RDUF proteins were predicted to be located in the chloroplasts (Table 1). Data from the 41 RDUFs in Gossypium were also investigated. The coding sequence (CDS) length of the 41 RDUF genes ranged from 953 bp (GrRDUF4) to 1312 bp (GbR-DUF7A). Most of the RDUF genes did not contain an intron, but GaRDUF3, GrRDUF4, GbRDUF2D, and GbRDUF7A had two exons each. The amino acid sequences of the RDUF proteins were also characterized. The MW ranged from 34.13 kDa (GbRDUF3A) to 46.92 kDa (GbRDUF7A), and the pI ranged from 5.17 (GbRDUF4A) to 9.27 (GaRDUF4). Subcellular location analysis showed that 23 of the 41 RDUF proteins were predicted to be located in the nucleus, and 13 RDUF proteins were predicted to be located in the chloroplasts (Table 1).

Chromosomal Localization and Gene Synteny Analysis of RDUF Genes
To reveal the homolog relationships between the RDUF genes in Gossypium, gene localization on chromosomes as well as gene synteny analysis were performed. The results showed that the RDUF genes in diploid cotton species were distributed on six chromosomes, with two GrRDUF genes localized on chr12 in G. raimondii. In tetraploid cotton species, we observed that seven RDUF genes of the At subgenome were distributed on seven chromosomes, with each RDUF gene located on one chromosome. The seven RDUF genes of the Dt subgenome were located on six chromosomes, with RDUF2D and RDUF3D located on the D04 chromosome ( Figure S3). Gene localization analysis revealed that most of the RDUF loci were highly symmetric between the At and Dt subgenomes. RDUF gene localization in the At subgenomic chromosomes consistently agreed with their homologs in the Dt subgenomic chromosomes. Nevertheless, RDUF2A and RDUF3A were located on the A03 and A04 chromosomes, respectively, but their homologs, RDUF2D and RDUF3D, were all located on the D04 chromosome. Gene synteny analysis was performed to understand the relationships of the RDUF genes between G. hirsutum and G. barbadense. We observed that homologous genes in or between each cotton genome were highly syntenic. The RDUF genes in G. hirsutum were consistent with their orthologs in G. barbadense ( Figure 2). We also investigated whether GhRDUF1A/D and GhRDUF4A/D; GhRDFU2A/3D and GhRDUF5A/D; and GhRDUF3A/2D, GhRDUF6A/D and GhRDUF7A/D might be duplicated genes ( Figure S4). The RDUF genes in the diploid cotton species (G. arboreum and G. raimondii) were twice more than those in cacao, reflecting that at least one WGD event occurred after cotton diverged from cacao. To explore the different selective constrains on the RDUF genes, the Ka/Ks ratios for the duplicated genes were calculated (Table S3). The majority of the Ka/Ks values of the RDUF gene pairs were <1, suggesting that the RDUF genes primarily experienced purified selection. Furthermore, 3 RDUF gene pairs exhibited neutral selection, whereas 13 RDUF duplicate gene pairs might have experienced positive selection in the 4 cotton species.

Conserved Structure and Domains in GhRDUF Proteins
A total of 41 RDUF proteins were used to understand phylogenetic relationships and gene structures. They were classified into three groups, in agreement with those in each cotton species ( Figure 3A). The RDUF proteins in the A genome and At subgenome, as well as their orthologs in the D genome and the Dt subgenome tended to consistently form one clade, indicating that the tetraploid cottons might have originated from a hybridization and WGD event of the two diploid cottons. Gene structure analysis displayed that most of the RDUF genes contained only one exon ( Figure 3B). Consistently, RDUF members with similar structures were grouped in the same clade.

Conserved Structure and Domains in GhRDUF Proteins
A total of 41 RDUF proteins were used to understand phylogenetic relationships and gene structures. They were classified into three groups, in agreement with those in each cotton species ( Figure 3A). The RDUF proteins in the A genome and At subgenome, as well as their orthologs in the D genome and the Dt subgenome tended to consistently form one clade, indicating that the tetraploid cottons might have originated from a hybridization and WGD event of the two diploid cottons. Gene structure analysis displayed that most of the RDUF genes contained only one exon ( Figure 3B). Consistently, RDUF members with similar structures were grouped in the same clade. The conserved domains were analyzed to understand the function of RDUF proteins in cotton. All the RDUF proteins in the four cotton species contained the zinc_ribbon_9 domain at the C-terminal, the RDUF117 domain at the N-terminal, the RING domain, and a low_complexity_region domain. Characteristically, the GrRDUF5, GbRDUF3A, GhRDUF3A, GbRDUF2D, and GhRDUF2D proteins contained an abbreviated RDUF117 domain. GbRDUF2D and GhRDUF2D contained the low_complexity_region domain at the N-terminal. GbRDUF2A, GhRDUF2A, GaRDUF2, GhRDUF1D, and GbRDUF1D possessed the low_complexity_region domain between the RING and the RDUF117 domain ( Figure 3C). Subsequently, the conserved motifs were detected using the MEME website. All 41 RDUF proteins contained motif_3, motif_5, motif_6, and motif_7. Most of the RDUF proteins contained motif 2, except for GaRDUF3. GrRDUF4 did not contain motif 1, whereas the RDUF3A/2D clade was without motif_4. GrRDUF4 and GrRDUF7 did not have motif_8, and the RDUF1A/D clade was without motif_9. However, only the The conserved domains were analyzed to understand the function of RDUF proteins in cotton. All the RDUF proteins in the four cotton species contained the zinc_ribbon_9 domain at the C-terminal, the RDUF117 domain at the N-terminal, the RING domain, and a low_complexity_region domain. Characteristically, the GrRDUF5, GbRDUF3A, GhRDUF3A, GbRDUF2D, and GhRDUF2D proteins contained an abbreviated RDUF117 domain. GbRDUF2D and GhRDUF2D contained the low_complexity_region domain at the N-terminal. GbRDUF2A, GhRDUF2A, GaRDUF2, GhRDUF1D, and GbRDUF1D possessed the low_complexity_region domain between the RING and the RDUF117 domain ( Figure 3C). Subsequently, the conserved motifs were detected using the MEME website. All 41 RDUF proteins contained motif_3, motif_5, motif_6, and motif_7. Most of the RDUF proteins contained motif 2, except for GaRDUF3. GrRDUF4 did not contain motif 1, whereas the RDUF3A/2D clade was without motif_4. GrRDUF4 and GrRDUF7 did not have motif_8, and the RDUF1A/D clade was without motif_9. However, only the RDUF6A/D clade had motif_10 ( Figure S5). The motif functions were annotated in the Pfam software, with motif_1 belonging to the RING finger domain, motif_3 and motif_6 belonging to the zinc-ribbon domain, and motif_4 belonging to the RDUF1117 domain (Table S4). In addition, the amino acid sequences of the 14 GhRDUF proteins and their homolog alignments were recorded using the Clustalx software. These results showed that the sequences of the conserved domains were highly similar. The conserved metal ligand positions and zinc (Zn 2+ ) coordinating amino acids were also investigated. We observed that all 14 GhRDUF proteins contained a His at metal ligand positions four and five, suggesting that the GhRDUF proteins belong to the RING-H2 type ( Figure S6).

Cis-Elements in the GhRDUFs Promoter and the Targeting TFs
Globally, over 90% of cotton crops are upland cotton, owing to its higher fiber production and environmental adaptation. The promoters of the 14 GhRDUF genes were submitted in the PlantCARE database to investigate the cis-elements. In total, 75 cis-elements were predicted to exist in the GhRDUF promoter regions, with more cis-elements implicated in the light response, including Box 4, G-Box, and MRE elements (Table S5). The participation of cis-elements in the environmental and hormone stress responses are highlighted in Figure 4. Among the 11 hormone response elements predicted in the GhRDUF promoters, ERE (cis-acting ethylene responsive element), MYC (cis-acting element involved in MeJA stress), and ABRE (cis-acting element involved in abscisic acid response) elements were the most abundant, indicating that the GhRDUF genes may primarily respond to ET, MeJA, and ABA stress. Nine environmental stress-related elements were identified. The majority of them were involved in anaerobic induction (ARE), plant defense signaling (W-box), and stress responses (STRE). In particular, more W-box cis-elements were found in the GhRDUF4D promoter region, implying that GhRDUF4D participates in cotton's response to disease.
Biomolecules 2021, 11, x FOR PEER REVIEW 12 of 24 RDUF6A/D clade had motif_10 ( Figure S5). The motif functions were annotated in the Pfam software, with motif_1 belonging to the RING finger domain, motif_3 and motif_6 belonging to the zinc-ribbon domain, and motif_4 belonging to the RDUF1117 domain (Table S4). In addition, the amino acid sequences of the 14 GhRDUF proteins and their homolog alignments were recorded using the Clustalx software. These results showed that the sequences of the conserved domains were highly similar. The conserved metal ligand positions and zinc (Zn 2+ ) coordinating amino acids were also investigated. We observed that all 14 GhRDUF proteins contained a His at metal ligand positions four and five, suggesting that the GhRDUF proteins belong to the RING-H2 type ( Figure S6).

Cis-Elements in the GhRDUFs Promoter and the Targeting TFs
Globally, over 90% of cotton crops are upland cotton, owing to its higher fiber production and environmental adaptation. The promoters of the 14 GhRDUF genes were submitted in the PlantCARE database to investigate the cis-elements. In total, 75 cis-elements were predicted to exist in the GhRDUF promoter regions, with more cis-elements implicated in the light response, including Box 4, G-Box, and MRE elements (Table S5). The participation of cis-elements in the environmental and hormone stress responses are highlighted in Figure 4. Among the 11 hormone response elements predicted in the GhRDUF promoters, ERE (cis-acting ethylene responsive element), MYC (cis-acting element involved in MeJA stress), and ABRE (cis-acting element involved in abscisic acid response) elements were the most abundant, indicating that the GhRDUF genes may primarily respond to ET, MeJA, and ABA stress. Nine environmental stress-related elements were identified. The majority of them were involved in anaerobic induction (ARE), plant defense signaling (W-box), and stress responses (STRE). In particular, more W-box cis-elements were found in the GhRDUF4D promoter region, implying that GhRDUF4D participates in cotton's response to disease.  Cis-elements bind TFs to regulate the precise initiation of gene transcription. Thus, the targeting TFs of GhRDUFs were predicted using the PlantRegMap server, with a total of 188 relationships identified, including 35 TF families and 14 GhRDUF genes ( Figure S7). It appears that GhRDUF1 homologous genes are regulated by additional TFs, such as those of ERF, Dof, and MYB. Meanwhile, many GhRDUF genes were regulated by the stress TFs of MYB, C2H2, and Dof, indicating that GhRDUF might regulate STREs.

GhRDUF Targeting miRNAs Predict the Regulatory Network of GhRDUF4D and gra-miR482c
miRNAs have been widely investigated for their role in the regulation of plant development, as well as in the abiotic and biotic stress responses. In this study, 17 putative miRNAs targeting 14 GhRDUF genes were predicted by the psRNATarget website, including 27 interaction relationships. We observed that GhRDUF2D was the most targeted, and that it interacted with four miRNAs. Most of the homologous GhRDUF genes were regulated by the same miRNAs, indicating that they have similar functions. We highlighted that GhRDUF4D is possibly regulated simultaneously by gra-miR482c and gra-miR482d. Furthermore, gra-miR482c regulated GhRDUF4D expression with two binding sites, implying a strong regulatory relationship ( Figure 5A). The regulatory network between GhRDUF4D and gra-miR482c was thus examined using a transient transformation assay. The five point mutations (gra-miR482c-mut) of gra-miR482 were used as the control ( Figure 5B). Different vector groups were infiltrated in tobacco leaves and photographed under an ultraviolet light 3 days after infiltration ( Figure 5C). We observed that the average fluorescent intensity in the GhRDUF4D-GFP and gra-miR482c precursor group was significantly impaired, compared with the control (Figure 5D,E). Moreover, the coexpression of the GhRDUF4D-GFP and gra-miR482c precursors generated an apparent reduction in the expression level of GFP ( Figure 5F). These results indicated that gra-miR482c mediates the degradation of GhRDUF4D.

Expression Profiling of GhRDUF Genes in Upland Cotton
Gene expression models are important for gene function analysis. The public expression datasets in upland cotton were used for gene expression analysis, including the different tissues, ovule, and fiber development stages [2]. Most of the homologous genes showed similar expression patterns. GhRDUF3A and GhRDUF7D genes were barely expressed in upland cotton, whereas GhRDUF1A/D, GhRDUF4A/D, and GhRDUF6A/D homologous genes were abundantly expressed in the root, stem, sepal, bract, early developing ovules (−3, 0, and 1 DPA), and the 10 and 20 DPA ovules. These results indicate that these genes are involved in cotton reproduction and vegetable development. Meanwhile, GhRDUF1A/D was preferentially expressed in 20 DPA fibers, suggesting its role in the secondary wall thickening of fiber ( Figure 6A).
The expression profiles of the GhRDUF genes under abiotic stress were also investigated, including in cold, heat, drought, and salt conditions. Extensive expression changes were detected under these stresses. The expression of the GhRDUF2A/3D genes was induced after 1 h of heat stress, and at 12 h under drought and salt stress. There was an impaired expression after 6, 12, and 24 h of cold stress. GhRDUF5A/D genes were upregulated at 1 h of heat and cold stress, at 3 and 12 h of cold stress, and at 24 h of cold and heat stress. GhRDUF1A/D genes were significantly elevated after 3 h of cold stress, and at 6 h under heat stress. In particular, GhRDUF4A/D genes were downregulated at 6 h under heat stress. GhRDUF6A/D genes responded to cold stress at 3, 6, 12, and 24 h. GhRDUF7A/D genes were abundantly expressed at 1, 3, 6, and 12 h when exposed to cold conditions ( Figure 6B). The expression patterns of GhRDUF genes under abiotic stress were consistent with the many environment response elements that were predicted in the promoter region ( Figure 4C).

Expression Profiling of GhRDUF Genes in Upland Cotton
Gene expression models are important for gene function analysis. The public expression datasets in upland cotton were used for gene expression analysis, including the different tissues, ovule, and fiber development stages [2]. Most of the homologous genes showed similar expression patterns. GhRDUF3A and GhRDUF7D genes were barely expressed in upland cotton, whereas GhRDUF1A/D, GhRDUF4A/D, and GhRDUF6A/D homologous genes were abundantly expressed in the root, stem, sepal, bract, early developing ovules (−3, 0, and 1 DPA), and the 10 and 20 DPA ovules. These results indicate that these genes are involved in cotton reproduction and vegetable development. Meanwhile, GhRDUF1A/D was preferentially expressed in 20 DPA fibers, suggesting its role in the secondary wall thickening of fiber ( Figure 6A). Plant hormone signals play a pivotal role in plant autoimmunity, particularly MeJA, MeSA, and ethylene. The response of the GhRDUF genes to these three hormones was investigated in this study. We observed that the highest expression of the GhRDUF genes was induced under the ethylene treatment (Figure 7). The transcript levels of GhRDUF1A/D, GhRDUF4A/D, and GhRDUF6A/D were upregulated 48 h after ethylene spraying. Many GhRDUF genes had impaired expression when exposed to the MeJA and MeSA conditions.  The expression profiles of the GhRDUF genes under abiotic stress were also investigated, including in cold, heat, drought, and salt conditions. Extensive expression changes were detected under these stresses. The expression of the GhRDUF2A/3D genes was induced after 1 h of heat stress, and at 12 h under drought and salt stress. There was an impaired expression after 6, 12, and 24 h of cold stress. GhRDUF5A/D genes were upregulated at 1 h of heat and cold stress, at 3 and 12 h of cold stress, and at 24 h of cold and heat stress. GhRDUF1A/D genes were significantly elevated after 3 h of cold stress, and at 6 h under heat stress. In particular, GhRDUF4A/D genes were downregulated at 6 h under heat stress. GhRDUF6A/D genes responded to cold stress at 3, 6, 12, and 24 h. GhRDUF7A/D genes were abundantly expressed at 1, 3, 6, and 12 h when exposed to cold conditions ( Figure 6B). The expression patterns of GhRDUF genes under abiotic stress were consistent with the many environment response elements that were predicted in the promoter region ( Figure 4C).
Plant hormone signals play a pivotal role in plant autoimmunity, particularly MeJA, MeSA, and ethylene. The response of the GhRDUF genes to these three hormones was investigated in this study. We observed that the highest expression of the GhRDUF genes was induced under the ethylene treatment (Figure 7). The transcript levels of

Overexpression of GhRDUF4D Enhanced the V. dahliae Resistance in Arabidopsis
The expression level of the homologous genes GhRDUF1A/D, GhRDUF4A/D, and GhRDUF7A/D were upregulated after V. dahliae infection based on RNA-Seq data, which GhRDUF4D was significantly increased at 12 and 24 h after V. dahliae infection ( Figure S8A) [11]. RT-qPCR revealed that GhRDUF4D expression was increased at 24 and 48 h after V. dahliae infection ( Figure S8B). Thus, GhRDUF4D was selected to validate the function by overexpression in Arabidopsis. We observed that GhRDUF4D transgenic plants were more resistant to V. dahliae infection than the control ( Figure 8A,C,D). The statistical disease index result showed that V. dahliae resistance in GhRDUF4D transgenic Arabidopsis was improved significantly, not only in MS medium but also in the soil ( Figure 8B,E). The relative fungal biomass was reduced in GhRDUF4D overexpressed Arabidopsis ( Figure 8F). These results indicate that V. dahliae resistance was improved by the overexpression of GhRDUF4D in Arabidopsis.
Other GhRDUF genes showed intricate expression patterns. For example, GhRDUF2A/D was elevated under MeJA and MeSA treatments at many time points, GhRDUF5A/D was abundantly expressed under MeJA treatment at 3 h, and GhRDUF7A/D was expressed under MeSA treatment at 3 and 24 h.

Overexpression of GhRDUF4D Enhanced the V. dahliae Resistance in Arabidopsis
The expression level of the homologous genes GhRDUF1A/D, GhRDUF4A/D, and GhRDUF7A/D were upregulated after V. dahliae infection based on RNA-Seq data, which GhRDUF4D was significantly increased at 12 and 24 h after V. dahliae infection ( Figure  S8A) [11]. RT-qPCR revealed that GhRDUF4D expression was increased at 24 and 48 h after V. dahliae infection ( Figure S8B). Thus, GhRDUF4D was selected to validate the function by overexpression in Arabidopsis. We observed that GhRDUF4D transgenic plants were more resistant to V. dahliae infection than the control ( Figure 8A,C,D). The statistical disease index result showed that V. dahliae resistance in GhRDUF4D transgenic Arabidopsis was improved significantly, not only in MS medium but also in the soil ( Figure 8B,E). The

Knockdown of GhRDUF4D Compromise V. dahliae Resistance in Upland Cotton
To investigate the biofunction of GhRDUF4D in the relationship between upland cotton and V. dahliae, the VIGS approach based on TRV was used to knockdown GhRDUF4D transcription levels in upland cotton. When the true leaves of the cotton seedlings were inoculated with TRV::CLA1, Agrobacteria revealed a photobleaching phenotype (Figure 9A), and the expression of GhRDUF4D was downregulated in TRV::GhRDUF4D plants ( Figure 9B). The seedlings were inoculated with V991 upon reaching the three-leaf stage. The results showed that the upland cotton seedlings with GhRDUF4D downregulation were more susceptible to V. dahliae than the control (TRV::00). Indeed, GhRDUF4D-silenced plants displayed more wilting and yellow leaves, and a deeper vascular discoloration in stem tissues ( Figure 9C,D). The quantification of the disease index in GhRDUF4D downregulated seedlings was higher than that of the control ( Figure 9E). Fungal DNA abundance detection and fungal recovery assays displayed that greater amounts of V. dahliae had colonized the stem tissue of GhRDUF4D-silenced plants than the control ( Figure 9F,G). relative fungal biomass was reduced in GhRDUF4D overexpressed Arabidopsis ( Figure 8F). These results indicate that V. dahliae resistance was improved by the overexpression of GhRDUF4D in Arabidopsis.

Knockdown of GhRDUF4D Compromise V. dahliae Resistance in Upland Cotton
To investigate the biofunction of GhRDUF4D in the relationship between upland cotton and V. dahliae, the VIGS approach based on TRV was used to knockdown GhRDUF4D transcription levels in upland cotton. When the true leaves of the cotton seedlings were inoculated with TRV::CLA1, Agrobacteria revealed a photobleaching phenotype ( Figure  9A), and the expression of GhRDUF4D was downregulated in TRV::GhRDUF4D plants ( Figure 9B). The seedlings were inoculated with V991 upon reaching the three-leaf stage. The results showed that the upland cotton seedlings with GhRDUF4D downregulation were more susceptible to V. dahliae than the control (TRV::00). Indeed, GhRDUF4D-silenced plants displayed more wilting and yellow leaves, and a deeper vascular discoloration in stem tissues ( Figure 9C,D). The quantification of the disease index in GhRDUF4D downregulated seedlings was higher than that of the control ( Figure 9E). Fungal DNA abundance detection and fungal recovery assays displayed that greater amounts of V. dahliae had colonized the stem tissue of GhRDUF4D-silenced plants than the control ( Figure 9F,G).
ROS is regarded as an indicator of disease resistance. Thus, the ROS levels of TRV::00 and TRV::GhRDUF4D cotton leaves were stained with DAB at 12 h after V. dahliae inoculation. The results showed that ROS accumulation in TRV::GhRDUF4D plants was significantly less than that of the TRV::00 plants ( Figure 9H). Callose deposition was evaluated by aniline blue staining. Compared with the control, the density of callose depositions ROS is regarded as an indicator of disease resistance. Thus, the ROS levels of TRV::00 and TRV::GhRDUF4D cotton leaves were stained with DAB at 12 h after V. dahliae inoculation. The results showed that ROS accumulation in TRV::GhRDUF4D plants was significantly less than that of the TRV::00 plants ( Figure 9H). Callose deposition was evaluated by aniline blue staining. Compared with the control, the density of callose depositions was reduced, with more dead cells observed in GhRDUF4D-silenced cotton leaves ( Figure 9I). These results indicate that the reduced expression of GhRDUF4D compromises V. dahliae resistance in upland cotton, reflecting the positive role of GhDUF4D in plant responses to pathogenic fungal infection.
was reduced, with more dead cells observed in GhRDUF4D-silenced cotton leaves ( Figure  9I). These results indicate that the reduced expression of GhRDUF4D compromises V. dahliae resistance in upland cotton, reflecting the positive role of GhDUF4D in plant responses to pathogenic fungal infection.

Evolution and Functional Diversification of GhRDUF Genes
The RING-containing proteins ubiquitously function in light signaling, phosphate homeostasis [57], anther development [58], plant defense, and pathogen response [25,[29][30][31][32][33]. We analyzed the evolution and expression of GhRDUF members in cotton plants to investigate the potential role of the RDUF genes. In addition, phylogenetic trees were constructed to determine the evolutionary relationship of RDUF proteins in plants. Gene duplication plays a major role in plant genome evolution. Based on our findings, the GhRDUF genes were duplicated frequently during cotton evolution. Only one RDUF gene was detected in each group of Arabidopsis thaliana and Theobroma cacao, however, two or three members were identified in diploid cottons (G. arboretum and G. raimondii). Of these, both the RDUF2A/3D and RDUF5A/D clades, and the RDUF1A/D and RDUF4A/D clades duplicated only once, whereas the RDUF3A/2D, RDUF6A/D, and RDUF7A/D clades duplicated twice during cotton evolution. Based on their chromosomal localization, the RDUF genes may have undergone segmental duplication in cotton evolution (Figure 3 and Figure S2). These results indicate that the GhRDUF genes duplicated after the cotton genus diverged from cacao.
The duplicated gene pairs possibly underwent three alternative fates during evolution, including nonfunctionalization, neofunctionalization, and subfunctionalization [56]. Thus, we investigated the expression patterns of the GhRDUF duplicated genes. The diverse expression patterns indicated that the GhRDUF genes have experienced functional diversification. We observed that the GhRDUF2A/3D genes were abundant during cotton growth and development, but its duplicated gene GhRDUF5A/D was barely expressed in the developing ovules or fibers ( Figure 6). The GhRDUF2A/3D and GhRDUF7A/D genes had reduced expression levels in different cotton tissues and the developing ovules and fibers, whereas its duplicated gene GhRDUF6A/D had increased expression in many tissues and developing ovules. These results imply that the duplicated GhRDUF genes may have experienced neofunctionalization and subfunctionalization during the evolution of upland cotton.

RDUF Function in Plants Adapting to Abiotic Stress
RING-containing proteins are key players in plant adaptation to abiotic stress. Many U-box E3 ubiquitin ligases participate in stress management, including TaPUB1 in salt stress [59], CaPUB1 in cold and drought stress [60], and AtPUB48 in thermotolerance response [61]. The RING-type E3 ligase, SDIR1, modulates the salt and drought stress responses via ABA signaling in Arabidopsis [29,30]. The RING-H2 type E3 ligase, OsSIRH2-14, degrades the salt-related protein OsHKT2;1 via the ubiquitin/26S proteasome and enhances salinity tolerance in rice [62]. The function of RING-DUF1117-containing proteins has rarely been reported. Expression analysis revealed that AtRDUF1 is more abundant during cold, heat, and salt stress [63]. AtRDUF1 positively regulates salt stress, and the suppression of AtRDUF1 and AtRDUF2 reduces the plant's tolerance to ABA-mediated drought stress in Arabidopsis [34,35]. These results indicate that RDUF is involved in the abiotic stress response. Thus, we analyzed the potential role of GhRDUFs in upland cotton's adaptation to abiotic stress. We detected many cis-elements involved in environmental adaptation in GhRDUF promoter regions: more than four ARE cis-elements were detected in GhRDUF3A/D, GhRDUF5A/D, and GhRDUF6A promoter regions; three STRE elements were found in GhRDUF5D and GhRDUF7A/D promoters; and three WUN-motif elements were found in GhRDUF2D and GhRDUF7A/D promoter regions (Figure 4). The TFgene regulatory network showed that many stress-related TFs are predicted to regulate GhRDUFs expression, including MYB, C2H2, and Dof ( Figure S7). It has been previously reported that ThMYB8 enhanced salt stress tolerance in Tamarix hispida and Arabidopsis [64], a C2H2-type ZFP gene, PeSTZ1, enhances freezing tolerance in Poplar [65], and MdDof54 promotes drought resistance in apples [66]. Thus, we investigated the expression patterns of the GhRDUF genes in response to cold, heat, drought, and salt stress. The GhRDUF7A/D homologous genes responded to cold and heat stress at 3 and 6 h, respectively, and showed increased expression under drought and salt stresses at 24 h, with three STRE elements detected in its promoter regions. In addition, the expression level of GhRDUF3D was upregulated at 1, 6, and 12 h when exposed to drought conditions, and two MBS elements were present in its promoter. These results suggest that the GhRDUF genes participate in cotton's adaptation to abiotic stress.

GhRDUF4D Participates in Resistance to V. dahliae
Hormones play a vital role in plant resistance to microorganisms, including pathogenic fungi. JA, SA, and ET are the primary hormones involved in plant inner immunity [6][7][8].
Many cis-elements involved in hormones were detected in the GhRDUF promoter regions (Figure 4). In this study, the expression levels of many GhRDUF genes were increased during hormonal responses (Figure 7). The expression of the GhRDUF2A/D homologous genes was elevated at 3 and 6 h under MeJA conditions, and induced at 12 and 24 h under MeSA conditions. In addition, we observed a decrease in the expression of GhRDUF4A/D, with respect to the MeJA and MeSA responses at many time points (Figure 7). RING-containing proteins have been widely reported to play a substantial role in plant disease resistance [67]. ATL9, a RING zinc finger protein which expression levels are induced by chitin, as well as the RING domain, are important for the resistant phenotype of ATL9 [68,69]. XB3 contains a RING finger motif and displays E3 ubiquitin ligase activity. Reduced XB3 expression is involved in resistance to Xanthomonas oryzae pv oryzae [32]. The overexpression of OsPUB41, a rice E3 ubiquitin ligase, enhanced tolerance to Xoo and Rhizoctonia solani infections in rice and Arabidopsis [70]. The mRNA level of RING-H2 type BRH1 had induced rapidly by the pathogen elicitor chitin [71]. The transcript level of another RING-H2 protein, ATL2, was induced after incubation with pathogen elicitors [72]. OsRFPH2-10, a RING-H2 finger E3 ubiquitin ligase, plays a role in the antiviral defense of rice against the rice dwarf virus infection [33]. In cotton, resistance to V. dahliae was shown to improve by the inhibiting of the E3 ubiquitin ligase activity of GhPUB17 [73]; however, the role of RING-RDUF1117 in biotic stress remains unclear. In the present study, we observed that the gene expression levels of GhRDUF1A/D, GhRDUF4A/D, and GhRDUF7A/D were elevated upon V. dahliae infection ( Figure S7). Consistently, more W-box elements were detected in the GhRDUF1A/D and GhRDUF4A/D promoter regions, particularly in the GhRDUF4D promoter ( Figure 4). This result suggests that GhRDUF4D is involved in cotton's defense to pathogen infection. To verify this speculation, overexpression and the VIGS approach were used to investigate the role of GhRDUF4D under V. dahliae infection. The results revealed that Arabidopsis plants that overexpress GhRDUF4D were more resistant to V. dahliae, and that GhRDUF4D downregulation in cotton plants made them more sensitive to V. dahliae infection, compared with the control. miRNA regulates target genes in plant development and aids in the plant's adaptation to the environment. In recent years, miRNAs have been characterized as key players in plants' responses to pathogens. Because NBS-LRR proteins have been associated with effector-triggered immunity, there is a good correlation between species with large numbers of NBS-LRR genes and those with miR482 superfamily members [74,75]. This result suggests that miR482 plays an important role in plant immunity. Studies have reported that miR482 suppressed Phytophthora and V. dahliae resistance in tomatoes and potatoes [76][77][78]. In cotton, it has been reported that miR482 regulates NBS-LRR defense genes during fungal pathogen infection [79]. In the present study, GhRDUF4D was predicted to be regulated by two miR482 (gra-miR482b and gra-miR482c). Furthermore, gra-miR482c regulated GhRDUF4D at two binding sites, suggesting that GhRDUF4D is regulated by miR482. Hence, a transient transformation assay was performed to examine the regulatory network of GhRDUF4D and gra-miR482c in tobacco leaves. The fluorescent intensity and GFP expression level in GhRDUF4D-GFP and gra-miR482c lanes displayed an apparent reduction compared with the GhRDUF4D-GFP and gra-miR482c-mut lanes ( Figure 5). These results indicate that GhRDUF4D might be regulated by gra-miR482c in vivo.

Conclusions
In conclusion, we detected and analyzed RDUF family members in cotton, and elucidated their roles in the biotic and abiotic stress responses. We found that the GhRDUF genes were mostly involved in the stress response, whereas the GhRDUF4D genes were involved in resistance to V. dahliae infections. The findings of this study contribute to the understanding of the RDUF genes in the development of resistant varieties. Moreover, GhRDUF4D might be a vital candidate gene applied for genetic transformation in cotton, or might be a marker gene used for V. dahliae resistance variety selection for cotton breeders.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/biom11081145/s1, Figure S1: Phylogenetic tree of RDUF protein family. At, Arabidopsis thaliana; Bn, Brassica napus; Ga, G. arboretum; Gr, G. rai-mondii; Gh, G. hirsutum; Gb, G. barbadense; Gm, Glycine max; Ha, Helianthus annuus; Os, Oryza sativa; Ta, Triticum aes-tivum; Tc, Theobroma cacao; Zm, Zea mays, Figure S2: Phylogenetic tree of RDUF family proteins in G. arboretum, G. raimondii, G. hirsutum, and G. barbadense, respectively, Figure S3: Gene location of RDUFs in G. arboretum, G. raimondii, G. hirsutum, and G. barbadense, Figure S4: The synteny relationship of RDUF genes in G. hirsutum, Figure S5: The conserved motifs identified in RDUF proteins; (A) phylogenetic tree of RDUF proteins; (B) the top ten conserved motifs in RDUF proteins; (C) logos of the ten conserved motifs in RDUF proteins, Figure S6: Amino acid sequences alignment of GhRDUF proteins. The conserved domains were circled by a rectangle with a dotted line; the eight metal ligands were highlighted by a red triangle, Figure S7: The predicted target TFs of GhRDUF genes. The predicted regulation TFs were identified with a yellow background in a circle, the target GhRDUF genes were marked with a blue background in a square. The interaction levels were displayed with different degrees, Figure S8 Table S1: Primers used in the present study, Table S2: Identification and nomenclature of RDUF members in other species, Table S3: Ka/Ks ratio of ortholog or paralog gene pairs, Table S4: Annotation of the conserved motifs in RDUF proteins, Table S5: cis-elements predicted in 14 GhRDUF promoter regions.