Identification and Comprehensive Structural and Functional Analyses of the EXO70 Gene Family in Cotton

The EXO70 gene is a vital component of the exocytosis complex and participates in biological processes ranging from plant cell division to polar growth. There are many EXO70 genes in plants and their functions are extensive, but little is known about the EXO70 gene family in cotton. Here, we analyzed four cotton sequence databases, identified 165 EXO70 genes, and divided them into eight subgroups (EXO70A–EXO70H) based on their phylogenetic relationships. EXO70A had the most exons (≥11), whereas the other seven each had only one or two exons. Hence, EXO70A may have many important functions. The 84 EXO70 genes in Asian and upland cotton were expressed in the roots, stems, leaves, flowers, fibers, and/or ovules. Full-length GhEXO70A1-A cDNA was homologously cloned from upland cotton (Gossypium hirsutum, G. hirsutum). Subcellular analysis revealed that GhEXO70A1-A protein was localized to the plasma membrane. A yeast two-hybrid assay revealed that GhEXO70A1-A interacted with GhEXO84A, GhEXO84B, and GhEXO84C. GhEXO70A1-A silencing significantly altered over 4000 genes and changed several signaling pathways related to metabolism. Thus, the EXO70 gene plays critical roles in the physiological functions of cotton.


Introduction
Vesicle transport is an extremely important cytological process in eukaryotes. It moves proteins, lipids, and other substances between the inner membrane system and the cells, and establishes cell polarity, secretion, growth, division, and wall formation [1]. Tethering is a key step in vesicle transport. Large multi-subunit tethering complexes were first discovered in yeast [2]. Exocysts tether different vesicles to the exocytosis site required for cellular secretion [3]. They are evolutionarily conserved octameric protein complexes composed of Sec3, Sec5, Sec6, Sec8, Sec10, Sec15, EXO70, and EXO84 [4,5]. EXO70 plays a key role in exocyst assembly [6]. It recruits exocysts on the target membrane and interacts with Rho protein to regulate SNARE complex assembly and activation there via SEC6. In this manner, EXO70 mediates polar exocytosis [7,8].
The exocyst subunits are encoded by a single gene in yeast and just a few genes in metazoans. However, 23 EXO70 subunits encoded by various loci have been identified in Arabidopsis [9,10]. The EXO70s in terrestrial plant genomes have even more copies. This phenomenon is unique to the EXO70 subunit of the exocyst [11]. In the fungal and animal genomes sequenced to date, only one EXO70 coding gene was found. Hence, multiple EXO70 gene copies are unique to higher terrestrial plants [12]. Certain EXO70 functions Genes 2021, 12, 1594 2 of 21 might have been alienated during evolution and participated in other biological processes besides membrane vesicle transport. Alternatively, various EXO70 functions are specialized and form different exocysts from other subunits that participate in specific membrane vesicle transport processes in the organization, carrier substrate, or transport link [12]. The expression profiles of the 23 members of the Arabidopsis EXO70 family have been analyzed. Expression of this gene family has the following characteristics: spatiotemporal expression specificity at the cell and tissue levels; no constitutive expression; and specific expression in dividing, growing, differentiating, and secreting cells [12]. Plant EXO70 gene family members participate at the transcriptional level in the biological processes of different cell types via cell-and tissue-specific expression patterns.
Cotton is a major global economic crop. It is a source of seed, fiber, oil, and medicine [30,31]. The development of novel high-quality cotton varieties is of great commercial importance. The high copy numbers and tissue-specific functions of the EXO70 gene in plants suggest that targeting EXO70 to construct high-quality cotton is feasible. To date, however, few studies have investigated the cotton EXO70 gene. The only report is that GhEXO70B1 may respond to stressors by mediating cell autophagy [32]. Here, we conducted evolutionary analyses, systematically named cotton EXO70 gene family members, and examined the functions of GhEXO70A1-A. It is believed that the discoveries of the present work will provide theoretical and empirical references for future research on cotton EXO70.

Identification of EXO70 Family Members in Cotton
Upland cotton (G. hirsutum, NAU), sea island cotton (Gossypium barbadense, G. barbadense, HAU), their common ancestor Asian cotton (Gossypium arboretum, G. arboretum, CRI), and Raymond cotton (Gossypium raimondii, G. raimondii, JGI) genomes and gene structure annotation file were downloaded from the CottonFGD genome database (https:// cottonfgd.org/, accessed on 15 December 2020). The tobacco genome and gene structure annotation file were downloaded from the Ensembl website (http://ensembl.org/index.html, accessed on 14 May 2021). The Arabidopsis EXO70 protein sequences were obtained from the TARE website (https://www.arabidopsis.org/, accessed on 12 August 2020). The CDS sequence from the downloaded genome and gene structure annotation file was extracted, and then translated into amino acid sequences. A local BLASTP search was completed to identify complete EXO70 genes, using A. thaliana EXO70 sequences as queries. A BLASTP search (e value: 1e−5) was used to obtain a dataset of EXO70 proteins. To identify all members of the EXO70 gene family, previous EXO70 gene sequences were analyzed using the Simple Modular Architecture Research Tool (NCBI CDD; https://www.ncbi.nlm.nih.gov/cdd/, accessed on 15 August 2020) and Pfam (http://pfam.janelia.org/, accessed on 15 August 2020). The EXO70 genes were annotated according to their corresponding orthologs in Arabidopsis and their chromosomal posi-tions in cotton. The genes in G. hirsutum and G. barbadense were named according to their homologs in each subgenome, where "A" and "D" represent the At and Dt subgenomes, respectively [33]. Assignment of each gene to the At subgenome or Dt subgenome was based on its DNA sequence homology with an A genome diploid species (G. arboretum) or a D genome diploid species (G. raimondii). The gene positions of GaEXO70, GrEXO70, GhEXO70, and GbEXO70 on the chromosome were determined from the genome data and were displayed using Tbtools software (https://www.tbtools.com/user-guide/installation, accessed on 20 August 2020). ExPASy software (https://www.expasy.org, accessed on 10 October 2020) was used to predict the theoretical molecular weights (MW) and isoelectric points (pI) [34]. The subcellular location of the EXO70 protein was established via the Cell-PLoc 2.0 (http://www.csbio.sjtu.edu.cn/bioinf/euk-multi-2/, accessed on 20 October 2020). The 2-kb region upstream of the EXO70 start codon was found on the CottonFGD website and submitted to the PlantCARE database (http://bioinformatics.psb. ugent.be/webtools/plantcare/html/, accessed on 10 November 2020) to localize cis-acting elements [35].
A total RNA extraction kit for plant polysaccharides and polyphenols (No. DP441; Tiangen Biotech, Beijing, China) was used to extract total RNA from roots, stems, leaves, flowers, cotyledons, and ovules at the 0 DPA, 10 DPA, and 20 DPA stages of upland cotton. A reverse transcription kit (No. E047; Novoprotein Scientific Inc., Summit, NJ, USA) was used according to the manufacturer's instructions to reverse transcribe the extracted RNA into cDNA. Quantitative GhEXO70A1-A primers (Table S1) were designed, and qRT-PCR was performed to analyze relative GhEXO70A1-A expression in different cotton tissues.

Cloning of GhEXO70A1-A Gene
AtEXO70A1 homologous genes were sought in CottonFGD. Vector NTI Advance (Thermo Fisher Scientific, Waltham, MA, USA) was used to compare the downloaded sequences. The group A genes with the highest homology were named GhEXO70A1-A. Primers were designed using Vector NTI Advance software according to the GhEXO70A1-A gene coding sequence (Table S1). The cDNA was amplified by PCR to obtain the GhEXO70A1-A gene, and the PCR product was purified and recovered. It was cloned into the B-zero vector, and it was sent out for sequencing (Sangon Biotech (Shanghai, China)).

Real-Time Quantitative PCR (RT-qPCR) Analysis
NovoStart SYBR qPCR SuperMix Plus (E096; Novoprotein) was used to perform qRT-PCR detection of the foregoing cDNA. GhUBQ7 was the internal reference gene. The quantitative primer statistics are shown in Table S1. qRT-PCR was performed in a LightCycler 480 II (Roche Diagnostics, Basel, Switzerland). The parameters of the PCR system and procedures were set according to the quantitative fluorescence kit instructions. There were three biological replicates per reaction. Relative gene expression was calculated by the 2 − Ct method [40]. GhEXO70A1-A gene expression patterns were analyzed in different cotton tissues.

Subcellular Localization Detection of GhEXO70A1-A
The GhEXO70A1-A coding region sequence was cloned into a 2300-YFP vector to obtain the p35S-GhEXO70A1-A-YFP plasmid (the primers are listed in Table S1), which was transformed into an Agrobacterium GV3101 competent strain (No. AC1001; Shanghai Weidi Biotechnology). A single clone was selected to identify the positive strain. Expanded culture was performed to inject tobacco leaves, which were then incubated in the dark at 28 • C for 16 h, illuminated for 48 h, and observed under a confocal laser microscope (No. FV1200; OLYMPUS, Tokyo, Japan). Green fluorescent protein signals in the tobacco leaf epidermal cells were photographed. The 2300-YFP empty vector was the control, and there were three independent biological replicates.

Virus-Induced GhEXO70A1-A Gene Silencing
An online website (http://vigs.solgenomics.net/, accessed on 20 March 2021) was used to determine the GhEXO70A1-A silencing sequence. Primers were designed for PCR amplification and ligated to the pCLCRVA vector to obtain the pCLCRV-GhEXO70A1-A recombinant plasmid. The pCLCRV-GhEXO70A1-A, pCLCRVB, and pCLCRVA plasmids were transformed into Agrobacterium tumefaciens GV3101, and the positive clones were selected. The pCLCRV-GhEXO70A1-A and pCLCRVA plasmids were mixed with equal amounts of pCLCRVB resuspension. After 3 h, a sterile syringe was used to inject the mixture into 14-d cotton seedling cotyledons. A CLCRV empty-load injection was the control (CK). After the positive seedling leaf yellowing phenotype appeared, total RNA was extracted and GhEXO70A1-A gene expression was determined by quantitative fluorescence PCR detection. The quantitative primers are listed in Table S1.

Transcriptome Analysis
Based on the results of the GhEXO70A1-A gene silencing expression analysis, total RNA was extracted from the fifth and sixth leaves of the virus-induced GhEXO70A1-A gene-silenced and control plants. The samples were subjected to gel electrophoresis, and transcriptome sequencing was performed at Beijing Nuohe Zhiyuan Technology Co. Ltd. (Beijing, China).

Identification and Analysis of the Phylogenetic Relationship of the EXO70 Gene Family in Cotton
The exocytosis complex subunits comprise mostly EXO70 gene family members. There are 23 and 47 EXO70 genes in the model dicotyledon Arabidopsis thaliana and the monocotyledon rice, respectively. Here, we identified 165 EXO70 genes among the four cotton subspecies included in the CottonFGD database, namely, G. hirsutum, G. barbadense, G. arboretum, and G. raimondii. There were 27, 26, 55, and 57 genes in G. arboretum, G. raimondii, G. barbadense, and G. hirsutum, respectively. We also identified 48 EXO70 genes in the tobacco database ( Figure 1A). A phylogenetic analysis of the evolutionary relationships of 23 Arabidopsis EXO70s, 41 rice EXO70s, 165 cotton EXO70s, and 48 tobacco EXO70s ( Figure 1C) showed that cotton EXO70 resembled Arabidopsis EXO70. Both of these plants are dicotyledons [9,41] and their EXO70s could be divided into eight categories. Relative to the monocotyledon rice, the dicotyledons lacked four EXO70 categories, such as EXO70I-EXO70L [42]. Hence, the EXO70 gene may be markedly differentiated between monocotyledons and dicotyledons. Monocotyledons possess more EXO70 genes than dicotyledons. Based on the phylogenetic tree, the grouping and naming of Arabidopsis, cotton EXO70s can be divided into eight subgroups (EXO70A-EXO70H) containing 12, 6, 29, 12, 27, 12, 22, and 45 genes, respectively ( Figure 1B).
According to the cotton EXO70 gene classification, we named the 57 EXO70 genes in upland cotton as GhEXO70A1-GhEXO70A2, GhEXO70B, GhEXO70C1-GhEXO70C5, GhEXO70D1-GhEXO70D2, GhEXO70E1-GhEXO70E6, GhEXO70F1-GhEXO70F2, GhEXO7 0G1-GhEXO70G4, and GhEXO70H1-GhEXO70H8. Groups A and D were represented by -A and -D, respectively. We predicted their genome locations, protein lengths, numbers of exons, isoelectric points, protein molecular weights, and subcellular locations. The numbers of exons widely varied among GhEXO70 genes. All four GhEXO70A genes had the most exons (≥11 each) ( Table 1). Further analysis of the exons of the EXO70 gene in Arabidopsis and rice showed that only Group A contained more exons in Arabidopsis and rice. The number of EXO70 in Arabidopsis A group was ≥9, and the number of EXO70 in rice A group was ≥12. Therefore, this phenomenon is not unique to cotton EXO70s, but is conserved in plants (Tables S2 and S3). The subcellular localization prediction results showed that GhEXO70 was mostly localized in the cell membrane, cytoplasm, or nucleus, which was consistent with the reported subcellular localization results of EXO70 from H. villosa [43] (Table 1).
Statistical analysis of the EXO70 gene distributions on the chromosomes revealed that there were relatively more EXO70 genes on chromosomes 5 and 9 in G. arboretum, G. barbadense, and G. hirsutum, but no EXO70 genes on chromosome 6 or 8. The EXO70 gene on chromosome 9 was distributed in G. raimondii, but that which was on chromosome 5 was not distributed. The EXO70 gene distributions on chromosomes 6 and 8 of G. raimondii (four and two, respectively) were the opposite of those for the other three cotton species (Figure 2; Table 2). The number of EXO70 genes in tetraploid cotton was nearly twice that in diploid cotton. The diploid cotton species (G. arboretum and G. raimondii) contained two EXO70As, one EXO70B, two EXO70Ds, and two EXO70Fs, whereas the tetraploid cotton species had twice these EXO70 gene copy numbers ( Table 3). The numbers of EXO70Cs, EXO70Gs, and EXO70Hs in tetraploid cotton were twice those in the autodiploid species and equal to the sum of the number in the allodiploid species (Table 3). In polyploid cotton, then,  The number of EXO70 genes in tetraploid cotton was nearly twice that in diploid cotton. The diploid cotton species (G. arboretum and G. raimondii) contained two EXO70As, one EXO70B, two EXO70Ds, and two EXO70Fs, whereas the tetraploid cotton species had twice these EXO70 gene copy numbers ( Table 3). The numbers of EXO70Cs, EXO70Gs, and EXO70Hs in tetraploid cotton were twice those in the autodiploid species and equal to the sum of the number in the allodiploid species (Table 3). In polyploid cotton, then, the number of EXO70 genes increases via genome polyploidization. Most GhEXO70 genes are highly parallel in the At group and Dt subgenome. The exception is that GhEXO70E2-D and GhEXO70E4-D have no homologs in the At subgenome, while GhEXO70H8-A has no homologs in the Dt subgenome, indicating that they may be lost during evolution.

Analysis of EXO70 Gene Structure in G. hirsutum
G. hirsutum is the major global cotton variety and was the focus of research attention here. Structural analysis of its 57 GhEXO70 genes showed that all of them had one or two exons except for GhEXO70A, which had 10 or 11 exons. All GhEXO70 genes with similar structures are grouped in the same clade. Moreover, the genes with closely related phylogeny in the same subgroup also had similar structures. Within the same subgroup, however, certain genes exhibited entirely different structures. GhEXO70E2-D contained two exons, while the other genes within the same subgroup had only one. Similarly, GhEXO70G4 contained two exons, whereas GhEXO70G1-GhEXO70G3 each contained a single exon. GhEXO70H2-A and GhEXO70H8-A each contained two exons while the other genes within the same subgroup had only one (Figure 3).
We used MEME online software to analyze the conserved motifs in the GhEXO70 protein and study its motif composition diversity and conservation. Figure 3 shows that 10 motifs (1-10) were identified, and each one was localized mainly to the C-terminal of the gene. Therefore, the C-terminal sequence of the GhEXO70 protein is highly conserved. The motif types revealed that the GhEXO70 gene members in subgroups A, B, C, and D were highly conserved and included all motifs. The GhEXO70s gene members in the other subgroups presented with obvious differences in motif type distribution, and some of them were lost. GhEXO70E2-D, GhEXO70E4-D, GhEXO70H1-D, and GhEXO70G2-D contained two, three, five, and six motifs, respectively. The functions of Motifs 1-10 have not been elucidated. Nevertheless, analysis of the conserved domains via the NCBI Conserved Domain Database (CCD) disclosed that they comprise the Exo70 domain ( Figure 3). and none of the other rice EXO70s has a transmembrane domain ( Figure S3). Among the prediction results of the transmembrane domain of cotton EXO70, only GhEXO70E2-D has a transmembrane domain, and the others have no transmembrane domain. Both Arabidopsis and cotton contain fewer EXO70s with transmembrane domains. As rice is a monocot, it may be evolving to have more EXO70, and there are more EXO70s with transmembrane domains.

Analysis of EXO70 Gene Expression Patterns in G. arboretum and G. hirsutum
Gene expression has spatiotemporal properties. The expression patterns of the various members of the EXO70 gene family may indicate the potential biological effects of these genes. We analyzed expression profile data in the CottonFGD and Cottongen (https://www.cottongen.org/, accessed on 22 March 2021) databases to clarify the spatiotemporal expression characteristics of the EXO70 gene. In G. hirsutum and G. arboretum, the EXO70 gene is commonly expressed in the roots, stems, leaves, flowers, fibers, and ovules and has spatiotemporal properties (Figures 4 and S4).  The PFam03081 domain at the C-terminus of the EXO70 protein is characteristic of the EXO70 superfamily [12], and all 165 predicted homologous clone EXO70 proteins possess it. However, the amino acid sequence lengths differed among EXO70 proteins and were in the range of 134-735 aa (average length = 618.736842105263 aa) ( Table 1). It was discovered that most GhEXO70 genes lacked a transmembrane (TM) structure. Only GhEXO70E2-D might possess a transmembrane region. Therefore, it may have evolved along with eukaryote evolution ( Figure S1). For the prediction of the transmembrane domain of Arabidopsis EXO70, the results showed that AtEXO70C1, AtEXO70C2, AtEXO70H5, AtEXO70H8, and AtEXO70A3 have transmembrane domains, but they are not obvious, and the other EXO70s have no transmembrane domains ( Figure S2). The prediction results of rice EXO70 show that OsEXO70A3, OsEXO70A4, OsEXO70H1a, OsEXO70H1b, OsEXO70H2, OsEXO70H3, OsEXO70H4, OsEXO70I3, OsEXO70I4, OsEXO70L1, OsEXO6K1, OsEXO70L, OsEXO70J1, OsEXO70J1, OsEXO70J2, OsEXO70J6, OsEXO70J8, OsEXO70K1, OsEXO7K2, and OsEXO70L1 have a transmembrane domain. In addition, OsEXO70A4 has a more obvious transmembrane domain at the C-terminus, and none of the other rice EXO70s has a transmembrane domain ( Figure S3). Among the prediction results of the transmembrane domain of cotton EXO70, only GhEXO70E2-D has a transmembrane domain, and the others have no transmembrane domain. Both Arabidopsis and cotton contain fewer EXO70s with transmembrane domains. As rice is a monocot, it may be evolving to have more EXO70, and there are more EXO70s with transmembrane domains.

Analysis of EXO70 Gene Expression Patterns in G. arboretum and G. hirsutum
Gene expression has spatiotemporal properties. The expression patterns of the various members of the EXO70 gene family may indicate the potential biological effects of these genes. We analyzed expression profile data in the CottonFGD and Cottongen (https://www.cottongen.org/, accessed on 22 March 2021) databases to clarify the spatiotemporal expression characteristics of the EXO70 gene. In G. hirsutum and G. arboretum, the EXO70 gene is commonly expressed in the roots, stems, leaves, flowers, fibers, and ovules and has spatiotemporal properties ( Figure 4 and Figure S4). GhEXO70A1-A, GhEXO70A1-D, GhEXO70B-A, GhEXO70B-D, GhEXO70D1-A, GhEXO70E1-A, GhEXO70E6-A, GhEXO70F2-D, and other genes in G. hirsutum are generally expressed at high levels and in various tissues. The GhEXO70H3-A gene is expressed mainly in the stamens, whereas the GhEXO70H5-A and GhEXO70H5-D genes are expressed mainly during the early stages of ovule development. In G. arboretum, the GaEXO70A1, GaEXO70E4, GaEXO70B, GaEXO70F1, GaEXO70F2, GaEXO70D1, GaEXO70E1 genes are generally highly expressed in different tissues, while GaEXO70A2 is expressed mainly in the 15D fibers ( Figure S4).
Ubiquitous EXO70 expression suggests that this gene is implicated in cotton growth and development. The GaEXO70A2 gene is expressed mainly in the fibers and might participate in cotton fiber development. The GhEXO70H3-A gene is expressed mainly in the stamens and may be associated with cotton fertility. The GhEXO70H5-A and GhEXO70H5-D genes are expressed mainly in the early stages of ovule development and could be involved in cotton seed formation.

EXO70 Gene Transcription Regulation Analysis
Spatiotemporal gene expression is regulated mainly by transcription factors (TFs) and epigenetics [44]. The observed differences in spatiotemporal expression of the various EXO70 genes may be related to their promoter specificity. We intercepted the 2-kb sequence upstream of the cotton EXO70 gene start codon and used the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/, accessed on 20 April 2021) to analyze the cis-elements in the promoter region. A total of 1081 cis-elements were predicted in the 57 GhEXO70 gene promoter regions. Of these, 10 and 11 categories were related to phytohormones and environmental stressors, respectively. The functions of the cis-elements in phytohormone and environmental stress response are highlighted in Figure 5. Among the predicted phytohormone response elements, the ERE, ABRE, and CGTCA motifs were the most abundant. Hence, the GhEXO70 gene might respond to ethylene, abscisic acid, and methyl jasmonate (MeJA) ( Figure 5A). Ten environmental stressrelated elements were identified and mainly involved drought stress (MYC), stress response (STRE), and anaerobic induction (ARE) ( Figure 5B). Therefore, the EXO70 gene may participate in the response to adversity. To further verify whether the above cis-acting elements are unique to cotton EXO70s, we also analyzed the EXO70s gene promoters in Arabidopsis and rice. The results indicate that the promoters of EXO70 genes in Arabidopsis and rice also contain cis-acting elements that respond to environmental stress and plant hormones ( Figures S5 and S6). It shows that this phenomenon is not unique to cotton EXO70, but is conserved in plants. 21, 12, x FOR PEER REVIEW 13 of 21 highly expressed in different tissues, while GaEXO70A2 is expressed mainly in the 15D fibers ( Figure S4).  sponse (STRE), and anaerobic induction (ARE) ( Figure 5B). Therefore, the EXO70 gene may participate in the response to adversity. To further verify whether the above cis-acting elements are unique to cotton EXO70s, we also analyzed the EXO70s gene promoters in Arabidopsis and rice. The results indicate that the promoters of EXO70 genes in Arabidopsis and rice also contain cis-acting elements that respond to environmental stress and plant hormones (Figures S5 and S6). It shows that this phenomenon is not unique to cotton EXO70, but is conserved in plants.

Expression Analysis and Subcellular Location of GhEXO70A1-A
The EXO70A1 gene is the most widely studied of all plant EXO70 genes. In Arabidopsis, AtEXO70A1 differentiates tubular molecules and regulates seed coat, root hair, stigma papillae development, and Kjeldahl band formation [45][46][47]. OsEXO70A1 plays important roles in vascular bundle differentiation and mineral nutrient assimilation [28]. In this study, we used GhEXO70A1-A in an experimental study on cotton EXO70 genes. We tested the GhEXO70A1-A gene expression patterns. GhEXO70A1-A was predominantly expressed in the stems, leaves, and flowers but its expression levels were low in the roots, ovules, and cotyledons ( Figure 6A). Subcellular GhEXO70A1-A protein localization predicted its roles in biological processes. Transient 35S-GhEXO70A1-A-GFP expression in tobacco produced a fluorescent signal. GhEXO70A1-A induced signals on the plasma membrane ( Figure 6B). Thus, GhEXO70A1-A was localized to the endomembrane system. This discovery was consistent with the roles of EXO70s in vesicle transport. sis, AtEXO70A1 differentiates tubular molecules and regulates seed coat, root hair, stigma papillae development, and Kjeldahl band formation [45][46][47]. OsEXO70A1 plays important roles in vascular bundle differentiation and mineral nutrient assimilation [28]. In this study, we used GhEXO70A1-A in an experimental study on cotton EXO70 genes. We tested the GhEXO70A1-A gene expression patterns. GhEXO70A1-A was predominantly expressed in the stems, leaves, and flowers but its expression levels were low in the roots, ovules, and cotyledons ( Figure 6A).

GhEXO70A1-A Protein Interaction Analysis
We used a yeast two-hybrid (Y2H) assay to explore the interactions among GhEXO70A1-A and the other subunits of the exocytosis complex. Plasmids containing GhEXO70A1-A and the other subunits of the exocytosis complex were co-transformed into Y2H Gold cells, which can grow on SD/-Leu-Trp. However, the cells were inoculated onto SD/-Ade/-His/-Leu/-Trp medium and only GhEXO70A1-A and GhEXO84A, Gh EXO84B, GhEXO84C co-transformed cells could grow on it and express X-α-Gal activity. GhEXO70A1-A interacted with EXO84A, EXO84B, and EXO84C (Figure 7), which means that it may function as a subunit of the exocytosis complex.
GhEXO70A1-A and the other subunits of the exocytosis complex. Plasmids containing GhEXO70A1-A and the other subunits of the exocytosis complex were co-transformed into Y2H Gold cells, which can grow on SD/-Leu-Trp. However, the cells were inoculated onto SD/-Ade/-His/-Leu/-Trp medium and only GhEXO70A1-A and GhEXO84A,Gh EXO84B, GhEXO84C co-transformed cells could grow on it and express X-α-Gal activity. GhEXO70A1-A interacted with EXO84A, EXO84B, and EXO84C (Figure 7), which means that it may function as a subunit of the exocytosis complex.

VIGS Silencing of GhEXO70A1-A Causes Changes in Signaling Pathways and Gene Expression
Gene silencing is an effective method of studying gene function. To explore the functions of GhEXO70A1-A in cotton, we constructed GhEXO70A1-A-gene-silenced cotton plants by virus-induced gene silencing (VIGS). qPCR demonstrated that the GhEXO70A1-A gene was successfully knocked down ( Figure 8A). We then used next-generation sequencing (NGS) technology to detect any changes in the transcriptome of GhEXO70A1-Asilenced leaves. We sorted out the expression of other EXO70 genes in the NGS data, and the results are shown in Figure S7; GhEXO70B-D, GhEXO70B-A, GhEXO70E6-D, GhEXO70F1-D, GhEXO70E3-D, GhEXO70H6-D, GhEXO70E1-D, GhEXO70H6-A, GhEXO70H3-D, GhEXO70D1-D, GhEXO70E3-A, GhEXO70C1-A, GhEXO70E5-A, and GhEXO70C5-A have significant changes, and there are no significant changes in other EXO70 genes. However, except for GhEXO70C1-A, GhEXO70H6-D, and GhEXO70H6-A, which decreased to 32.1%, 46.9%, and 56.6% of the control, all other genes fell to more than 60% of the control, and the fold increase was also less than 1. Although the three genes GhEXO70C1-A, GhEXO70H6-D, and GhEXO70H6-A declined slightly, their

VIGS Silencing of GhEXO70A1-A Causes Changes in Signaling Pathways and Gene Expression
Gene silencing is an effective method of studying gene function. To explore the functions of GhEXO70A1-A in cotton, we constructed GhEXO70A1-A-gene-silenced cotton plants by virus-induced gene silencing (VIGS). qPCR demonstrated that the GhEXO70A1-A gene was successfully knocked down ( Figure 8A). We then used next-generation sequencing (NGS) technology to detect any changes in the transcriptome of GhEXO70A1-A-silenced leaves. We sorted out the expression of other EXO70 genes in the NGS data, and the results are shown in Figure S7; GhEXO70B-D, GhEXO70B-A, GhEXO70E6-D, GhEXO70F1-D, GhEXO70E3-D, GhEXO70H6-D, GhEXO70E1-D, GhEXO70H6-A, GhEXO70H3-D, GhEXO7 0D1-D, GhEXO70E3-A, GhEXO70C1-A, GhEXO70E5-A, and GhEXO70C5-A have significant changes, and there are no significant changes in other EXO70 genes. However, except for GhEXO70C1-A, GhEXO70H6-D, and GhEXO70H6-A, which decreased to 32.1%, 46.9%, and 56.6% of the control, all other genes fell to more than 60% of the control, and the fold increase was also less than 1. Although the three genes GhEXO70C1-A, GhEXO70H6-D, and GhEXO70H6-A declined slightly, their expression abundance was also very low. The above results show that the knockdown of GhEXO70A1-A by VIGS does affect the expression of other EXO70 genes, but the effect is not significant after analysis. The changes in differential genes should be mainly caused by the changes in GhEXO70A1-A.
Correlation analyses among samples disclosed significant differences between the GhEXO70A1-A-silenced (EXO70A1) and the control (VIGS-CK) groups ( Figure 8B). Thus, GhEXO70A1-A silencing in cotton altered the gene expression profiles. Differentially expressed genes (DEG) were those that met the criteria of |log2(Fold Change)| ≥ 1 and p ≤ 0.05. A total of 3264 upregulated and 1103 downregulated genes were screened, as shown in a volcano graph ( Figure 8C) and a heat map ( Figure 8D). Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment of the DEGs ( Figure 8E) displayed 13 pathways with p < 0.01. These included photosynthesis antenna protein, phenylpropane biosynthesis, flavonoid biosynthesis, starch and sucrose metabolism, circadian rhythmplant, keratin, cork and wax biosynthesis, steroid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, glutathione metabolism, cyano-amino acid metabolism, photosynthesis, and glucosinolate biosynthesis (Table 4). GSEA results showed that GhEXO70A1-A was significantly related to photosynthesis antenna protein, photosynthesis, and circadian rhythm-plants ( Figure 8F-H). Of the 13 significantly different pathways, all except for circadian rhythm-plants were related to metabolism. Therefore, cotton leaf GhEXO70A1-A may regulate biochemical anabolism and catabolism.

Evolutionary Relationships of the Cotton EXO70 Gene Family
The evolutionary relationships of the EXO70 gene family in Asian, Raymond, upland, and sea island cotton were inferred according to their total numbers, classifications, chromosome distributions, and structures. Allotetraploid cotton was crossed from its ancestors G. raimondii and G. arboretum about 5 million years ago [48]. There are about twice as many EXO70 genes in tetraploid cotton (57 and 55, respectively) as there are in diploid cotton.
Hence, the former underwent a round of polyploidization event. Polyploidization leads to rapid, extensive genetic and epigenetic changes in the genome. It is associated with many molecular and physiological adjustments and significant gene losses in heterotetraploid cotton after domestication. In tetraploid cotton, the EXO70 gene family underwent neither large-scale amplification nor reduction according to chromosome location analyses.
Phylogenetic analyses revealed that all eight EXO70 subgroups (EXO70A-EXO70H) are represented in all four cotton species. Genes within the same subgroup had similar structures. Subgroup EXO70A consisted of numerous exons, whereas there were relatively fewer in the other seven subgroups. During long-term natural selection, many EXO70 genes differentiated and evolved to contend with various stressors. A study on Arabidopsis EXO70 demonstrated that recurring gene loss was nonrandom. Genes involved in DNA repair were comparatively more prone to loss, while those implicated in signal transduction and transcription were preferentially retained [49]. Hence, the DNA repair function appears to be absent in the cotton EXO70I subgroup. Other plant species should be examined to determine whether the EXO70I branch is unique to monocots.

Biological Processes Implicating GhEXO70A1-A
EXO70 is one of eight exocyst subunits. It participates in the movement of membranerelated substances and is vital for vesicle transport, cell secretion, growth, division, and other processes [1]. EXO70 genes are vital to plant growth, development, and metabolism. The study of EXO70 is of great significance to cotton development and metabolism, as this crop is commercially important. Unlike yeast and animals, plants harbor dozens of EXO70 genes. The latest research shows that AtEXO70A1 recruits the entire complex to the cytoplasmic membrane by binding to negatively charged phospholipids, providing an important research basis for the study of plant cell polarity and morphogenesis [50]. A functional defect in the EXO70C2 gene of Arabidopsis affected pollen tube growth and led to male-specific transmission defects [22]. EXO70H4 and PMR4-dependent corpus callosum deposition in trichomes is necessary for cell wall silicification [23]. An EXO70B1 knockout mutant in Arabidopsis exhibited impaired light-induced stomatal opening [24]. An OsEXO70L2 mutation caused defects in rice root development [29]. Evidently, the EXO70 gene has been continuously amplified over the course of plant evolution and displays different physiological functions in various plant tissues and organs.
Transcriptome sequencing indicated that after GhEXO70A1-A knockdown in cotton leaf, the expression levels of over 4000 genes significantly changed. Therefore, GhEXO70A1-A-gene silencing triggered substantial alterations in the cotton leaf gene expression profile. A DEG function enrichment analysis demonstrated that GhEXO70A1-A function in cotton leaf is strongly correlated with metabolism and especially photosynthesis. Plant metabolism is regulated mainly by various phytohormones, such as auxins, cytokinins, jasmonic acid, abscisic acid, and so on [51,52]. GhEXO70A1-A silencing in cotton induces extensive modifications of the metabolic pathways. Thus, GhEXO70A1-A loss could affect phytohormone biosynthesis and/or distribution. Studies on Arabidopsis reported that EXO70 is vital for recycling plasma membrane proteins, including the auxin efflux carrier PIN [21,27,53]. The differential functions of various EXO70 genes in various plant tissues suggest that selective control of specific plant organs is feasible through targeted EXO70 gene regulation. However, further experimental research in this area must be conducted.

Conclusions
In the present study, the EXO70 genes of Asian cotton, Raymond cotton, G. hirsutum, and sea island cotton were analyzed by bioinformatics. A total of 165 EXO70 genes in eight categories were identified. An evolutionary analysis showed that the EXO70 gene multiplication patterns were nearly identical in upland cotton and sea island cotton, and the same EXO70 gene types were highly conserved in both varieties. However, EXO70 gene multiplication in these varieties underwent divergent evolution. Representative research on GhEXO70A1-A revealed that the EXO70 genes in cotton may be widely involved in metabolism and could also affect phytohormone biosynthesis and/or distribution.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/genes12101594/s1, Table S1. Primers used in this study. Table S2. Predicted EXO70 genes from Arabidopsis thaliana and the corresponding proteins. Table S3. Predicted EXO70 proteins and the corresponding gene from Oryza sativa. Figure S1. Analysis of transmembrane domains of EXO70 gene family in upland cotton. Figure S2. Analysis of transmembrane domains of EXO70 gene family in Arabidopsis thaliana. Figure S3. Analysis of transmembrane domains of EXO70 gene family in rice. Figure S4. EXO70 family gene expression patterns in different tissues and organs of asian cotton. Figure S5. Cis-acting elements in Arabidopsis EXO70s promoter. Figure S6. Cis-acting elements in rice EXO70s promoter. Figure S7. Expression analysis of other changed EXO70 genes in NGS data.