Genome-Wide Identification and Evolutionary Analysis of Gossypium YTH Domain-Containing RNA-Binding Protein Family and the Role of GhYTH8 in Response to Drought Stress

YTH domain-containing proteins are one kind of RNA-binding protein involved in post-transcriptional regulation and play multiple roles in regulating the growth, development, and abiotic stress responses of plants. However, the YTH domain-containing RNA-binding protein family has not been previously studied in cotton. In this study, a total of 10, 11, 22, and 21 YTH genes were identified in Gossypium arboreum, Gossypium raimondii, Gossypium barbadense, and Gossypium hirsutum, respectively. These Gossypium YTH genes were categorized into three subgroups by phylogenetic analysis. The chromosomal distribution, synteny analysis, structures of Gossypium YTH genes, and the motifs of YTH proteins were analyzed. Furthermore, the cis-element of GhYTH genes promoter, miRNA targets of GhYTH genes, and subcellular localization of GhYTH8 and GhYTH16 were characterized. Expression patterns of GhYTH genes in different tissues, organs, and in response to different stresses were also analyzed. Moreover, functional verifications revealed that silencing GhYTH8 attenuated the drought tolerance in the upland cotton TM-1 line. These findings provide useful clues for the functional and evolutionary analysis of YTH genes in cotton.


Introduction
The post-transcriptional regulation of gene expression is of great importance for various biological processes, such as plant growth and development, and environmental stress responses. Many processes are involved in post-transcriptional regulation, including pre-mRNA processing, polyadenylation, mRNA stability, and RNA transport [1,2]. These processes could be directly mediated by RNA-binding proteins (RBPs) through interacting with target RNA molecules or indirectly mediated by RBPs by regulating the function of other regulatory factors [3]. The study of RBPs in plants was slower than other organisms due to the lack of appropriate plant-derived in vitro systems for the research of posttranscriptional gene regulation. More than 200 RBPs were identified in Arabidopsis thanliana, and about 250 RBPs were identified in rice. Most of these identified RBPs are plant-specific and can act in plant-specific functions, such as the involvement in hormone responses and various environmental stresses [4,5].
YTH domain-containing proteins are one class of RBPs. In 1998, Imai Y et al. reported one member of the YTH-containing proteins in rats named YT521 [6]. YT521 was an RNA

Phylogenetic Analysis of YTH Genes
To study the evolutionary relationships among 64 cotton YTH genes, protein sequences from the 4 cotton species and Arabidopsis were used to build a phylogenetic tree using the ML method. The cotton YTH genes were categorized into three subgroups from Group I to Group III. Genes that were divided into the same group were more closely related. Group I to Group III contains 30, 18, and 16 cotton YTH genes, respectively. There were 5, 3, and 2 GaYTHs contained in Group I to Group III in G. arboretum, respectively, and 5, 3, and 3 GrYTHs in G. raimondii. While 10, 6, and 6 GbYTHs were classified into Group I to Group III, respectively, in G. barbadense, and 10, 6, and 5 GhYTHs, in G. hirsutum (Figure 1). The phylogenetic tree displayed that G. hirsutum and G. barbadense tend to be assigned to the same branch, suggesting a closer relationship between these two species. The phylogenetic tree resolved YTH genes into three groups. The phylogenetic ML tree was constructed using the amino acid sequence with a 1000 bootstrap value with MEGA-X software (version 7.0). YTH genes were categorized into threegroups from Group I to Group III. Ga, Gossypium arboretum; Gr, Gossypium raimondii; Gh, Gossypium hirsutum; Gb, Gossypium barbadense; At, Arabidopsis thaliana.

Motifs and Gene Structure Analysis
To explore the evolutionary patterns and classification of YTHs in cotton, the unrooted phylogenetic tree was constructed. The resultant phylogenetic tree of four cotton species divided the YTH proteins into three main groups (Groups I-III) ( Figure  2A). Based on the phylogenetic result, we analyzed the motifs and the sequence structure Phylogenetic analysis of YTH gene family in Gossypium and Arabidopsis. The phylogenetic tree resolved YTH genes into three groups. The phylogenetic ML tree was constructed using the amino acid sequence with a 1000 bootstrap value with MEGA-X software (version 7.0). YTH genes were categorized into threegroups from Group I to Group III. Ga, Gossypium arboretum; Gr, Gossypium raimondii; Gh, Gossypium hirsutum; Gb, Gossypium barbadense; At, Arabidopsis thaliana.

Motifs and Gene Structure Analysis
To explore the evolutionary patterns and classification of YTHs in cotton, the unrooted phylogenetic tree was constructed. The resultant phylogenetic tree of four cotton species divided the YTH proteins into three main groups (Groups I-III) (Figure 2A). Based on the phylogenetic result, we analyzed the motifs and the sequence structure of the identified YTHs. To explore the conserved motif organization in YTH proteins, the MEME program was used for the conserved motif analysis. A total of 10 motifs were identified in the YTH proteins. Motifs 1-4 and 6-8 were found in all YTH proteins. Motifs 1-9 were identified in Group I, and especially Motif 5 was only found in Group I ( Figure 2B). Furthermore, in most of the YTH proteins, the YTH domain is the only recognizable module that is in line with those of other species [13]. In addition, dnaA superfamily domain and PHA03377 superfamily domain were also found in several YTH proteins. The dnaA superfamily domain was only found in Group I, while the PHA03377 superfamily domain was only found in Group III (Supplementary Figure S1).
To reveal the consistency of the exon-intron pattern in the phylogenetic groups, we performed the gene structure analysis on the cotton YTH genes. Our data revealed that the number of exons varied from 6~9 in G. arboreum, 6~10 in G. raimondii, 6~11 in G. barbadense, and 6~12 in G. hirsutum in YTH gene family. The exons and introns arrangement uncovered the evolutionary relationships between different gene family members. The YTH genes that have a similar exon-intron pattern tend to be clustered into a similar branch. Particularly, in Group III, most of the exon numbers of the YTH genes were 8, except for GhYTH18 ( Figure 2C, Table S1). Our results showed that the closely associated genes are more structurally similar and differ in the length of intron and exons.

Chromosome Distribution and Synteny Analysis of YTH Genes
The chromosomal distribution of YTH genes was identified according to the genome annotation file. YTH genes were distributed on chromosomes 1-5, 9-11, and 13. GaYTH genes were localized across chromosomes 1-3, 8, 10, 11, and 13. GrYTH genes were localized across chromosomes 2, 4, 5, 7, 9, 11, and 13. GhYTHs were distributed on At chromosomes 1-3, 5, 8, 11, 13, and Dt chromosomes 1, 2, 5, 8, 10, 11, and 13 with 10 and 11 genes, respectively. GbYTHs were located on At chromosomes 1-3, 5, 8, 10, 11, 13, and Dt chromosomes 1, 2, 5, 8, 10, 13 with 11 and 9 genes, respectively. Another two GbYTHs (GbYPH20 and GbYPH21) were identified on unanchored scaffolds. Results showed that the chromosome distribution pattern between GbYTHs and GhYTHs was similar ( Figure 3). To further explore the phylogenetic mechanisms of the cotton YTH gene family, we performed the synteny analysis of YTH genes between G. hirsutum and the three other representative cotton species. A total of 32 GhYTH genes showed syntenic relationship with those in G. barbadense, followed by G. raimondii (19) and G. arboretum (18). Some To further explore the phylogenetic mechanisms of the cotton YTH gene family, we performed the synteny analysis of YTH genes between G. hirsutum and the three other representative cotton species. A total of 32 GhYTH genes showed syntenic relationship with those in G. barbadense, followed by G. raimondii (19) and G. arboretum (18). Some GhYTH genes were found to have at least three syntenic gene pairs, suggesting that these genes may have played an essential role during evolution ( Figure 4).
To better reveal G. hirsutum replication events, we identified all 32 homologous genes in G. hirsutum and compared them with those in G. barbadense. We calculated the levels of Ka and differences in Ks (Table S2). Ka/Ks = 1.0 suggested pseudogenes caused by neutral selection, Ka/Ks < 1 implied purifying selection, and Ka/Ks > 1 denoted positive selection and accelerated evolution. The Ka/Ks ratios of most duplicate gene pairs were less than one (Table S2), indicating that YTH gene purification and positive selection existed during cotton evolution.

Cis-Regulatory Element Analysis of YTH Genes
Cis-regulatory elements play important roles in gene expression regulation [22,23]. Understanding the cis-elements involved in the regulation of YTH genes will facilitate the study of the regulation mechanism and the putative functions of these YTH genes. The 2000 bp region upstream of the start codon of each YTH gene from 4 cotton species was extracted as the promoter region. The cis-acting elements, including salt responsive element, ethylene response element, and abscisic acid responsive element, were identified in the promoters of YTH genes of four cotton species. The number of these cis-acting elements in different species was basically consistent ( Figure 5A). Inducible promoters can respond immediately to stimulus signals and be involved in the related gene expression regulation under particular environmental and stimulus conditions. Furthermore, we identified a total of 15 inducible promoters of YTH genes in cotton ( Figure 5B).

Cis-Regulatory Element Analysis of YTH Genes
Cis-regulatory elements play important roles in gene expression regulation [22,23]. Understanding the cis-elements involved in the regulation of YTH genes will facilitate the study of the regulation mechanism and the putative functions of these YTH genes. The 2000 bp region upstream of the start codon of each YTH gene from 4 cotton species was extracted as the promoter region. The cis-acting elements, including salt responsive element, ethylene response element, and abscisic acid responsive element, were identified in the promoters of YTH genes of four cotton species. The number of these cis-acting elements in different species was basically consistent ( Figure 5A). Inducible promoters can respond immediately to stimulus signals and be involved in the related gene expression regulation under particular environmental and stimulus conditions. Furthermore, we identified a total of 15 inducible promoters of YTH genes in cotton ( Figure 5B).

Prediction of YTH Genes Targeted by MiRNAs in G. Hirsutum
Unlike cis-regulatory elements, the miRNAs were considered transac targeting the mRNA. The relationships between the miRNAs and YTHs in predicted in the psRNATarget database (Table S3). In total, 18 miRNA targe genes were identified. Among 21 identified GhYTHs, 17 of them were f targeted by miRNAs, except for GhYTH10, GhYTH16, GhYTH17, and GhYTH GhYTH2, GhYTH3, GhYTH7, GhYTH13, and GhYTH18 were targeted by thr respectively. GhYTH12 and GhYTH19 were targeted by two miRNAs, respec of the other nine GhYTHs was only targeted by only one miRNA. Almost targeted a single or two GhYTHs. ghr-7510R and ghr-7494 targeted fiv GhYTHs, respectively ( Figure 6).

Prediction of YTH Genes Targeted by MiRNAs in G. Hirsutum
Unlike cis-regulatory elements, the miRNAs were considered transacting factors targeting the mRNA. The relationships between the miRNAs and YTHs in cotton were predicted in the psRNATarget database (Table S3). In total, 18 miRNA targets of GhYTH genes were identified. Among 21 identified GhYTHs, 17 of them were found to be targeted by miRNAs, except for GhYTH10, GhYTH16, GhYTH17, and GhYTH21. GhYTH1, GhYTH2, GhYTH3, GhYTH7, GhYTH13, and GhYTH18 were targeted by three miRNAs, respectively. GhYTH12 and GhYTH19 were targeted by two miRNAs, respectively. Each of the other nine GhYTHs was only targeted by only one miRNA. Almost all miRNAs targeted a single or two GhYTHs. ghr-7510R and ghr-7494 targeted five and three GhYTHs, respectively ( Figure 6).

Prediction of YTH Genes Targeted by MiRNAs in G. Hirsutum
Unlike cis-regulatory elements, the miRNAs were considered transacting factors targeting the mRNA. The relationships between the miRNAs and YTHs in cotton were predicted in the psRNATarget database (Table S3). In total, 18 miRNA targets of GhYTH genes were identified. Among 21 identified GhYTHs, 17 of them were found to be targeted by miRNAs, except for GhYTH10, GhYTH16, GhYTH17, and GhYTH21. GhYTH1, GhYTH2, GhYTH3, GhYTH7, GhYTH13, and GhYTH18 were targeted by three miRNAs, respectively. GhYTH12 and GhYTH19 were targeted by two miRNAs, respectively. Each of the other nine GhYTHs was only targeted by only one miRNA. Almost all miRNAs targeted a single or two GhYTHs. ghr-7510R and ghr-7494 targeted five and three GhYTHs, respectively ( Figure 6).

Tissue Expression Patterns of YTH Genes in G. hirsutum
Gene expression patterns are closely related to their functions. To explore the spatialtemporal expression variations in YTH genes across various anatomical tissues, we analyzed the transcriptome data of G. hirsutum during the development of ovule and fiber (Table S4). In ovules, the expression levels of GhYTH2/5/8/10/15/19/21 were higher from 5 DPA to 20 DPA ( Figure 7A). During fiber development, the expression levels of GhYTH5 and GhYTH15 increased substantially from 10 DPA to 25 DPA. The expression levels of GhYTH2/12/14/16, GhYTH13, and GhYTH3/18 were relatively high at 10 DPA, 15 DPA, and 20 DPA, respectively. The expression of GhYTH4/6/7/8/10/12/17 was relatively high at 10 DPA and 15 DPA ( Figure 7B). The results indicated that each GhYTH gene might have a specific function at different time points, including the initiation and elongation stages of cotton fibers.
In addition, the expression profiles of GhYTH genes from ovule, fiber, torus, stem, sepal, root, pistil, petal, leaf, filament, bract, and anther were also studied. GhYTH genes displayed dynamic expression profiles in these tissues and organs (Table S4). Briefly, the expression of GhYTHs was relatively low in leaves, filaments, petals, and bracts. GhYTH2/3/12/13 in anther, GhYTH4/14/18 in root, GhYTH8/19 in sepal, GhYTH16/18 in stem, and GhYTH7/8/17/19 in torus were highly expressed ( Figures 7C and 8). The results showed that GhYTH genes might be involved in the growth and development processes of different organs in cotton.

Tissue Expression Patterns of YTH Genes in G. hirsutum
Gene expression patterns are closely related to their functions. To explore spatial-temporal expression variations in YTH genes across various anatomical tissu we analyzed the transcriptome data of G. hirsutum during the development of ovule a fiber (Table S4). In ovules, the expression levels of GhYTH2/5/8/10/15/19/21 were hig from 5 DPA to 20 DPA ( Figure 7A). During fiber development, the expression levels GhYTH5 and GhYTH15 increased substantially from 10 DPA to 25 DPA. The express levels of GhYTH2/12/14/16, GhYTH13, and GhYTH3/18 were relatively high at 10 DPA, DPA, and 20 DPA, respectively. The expression of GhYTH4/6/7/8/10/12/17 was relativ high at 10 DPA and 15 DPA ( Figure 7B). The results indicated that each GhYTH ge might have a specific function at different time points, including the initiation a elongation stages of cotton fibers.

Expression Patterns of YTH Genes in G. hirsutum under Stress Conditions
We analyzed the transcriptome data of G. hirsutum to uncover the GhYTH participation in response to cold, hot, salt, drought, and salt stress conditions (Table S5). Under cold stress, the expression levels of GhYTH11/18 reached a peak at 6 h after treatment and then decreased, the expression of GhYTH3/5/6/13/14/15/16/19/20 peaked at 12 h after treatment then decreased, and the expression of GhYTH1/2/4/7/8/9/10/12/17/21 was peaked at 24 h post-treatment ( Figure 9A). Under hot stress, the expression of GhYTH genes peaked at 1 h and then decreased except for GhYTH6/12/14/16. However, the expression trend of GhYTH6/14/16 continuously increased and peaked at 24 h post treatment ( Figure 9B). Under drought stress, GhYTH genes exhibited higher expression in the first 6 h after treatment and then showed decreased expression except for GhYTH1/4/5/11/13/14/15/20. The expression level of GhYTH2/6/8/9/10/12/16/19 gradually increased post-drought treatment, indicating their possible roles in response to drought stress. (Figure 9C). Under salt stress, GhYTH6/16/20 were highly expressed at 1 h post-treatment ( Figure 9D). In brief, our data showed that the expression level of most YTH genes was induced in response to multiple stresses.

Expression Patterns of YTH Genes in G. hirsutum under Stress Conditions
We analyzed the transcriptome data of G. hirsutum to uncover the GhYTH participation in response to cold, hot, salt, drought, and salt stress conditions (Table S5). Under cold stress, the expression levels of GhYTH11/18 reached a peak at 6 h after treatment and then decreased, the expression of GhYTH3/5/6/13/14/15/16/19/20 peaked at 12 h after treatment then decreased, and the expression of GhYTH1/2/4/7/8/9/10/12/17/21 was peaked at 24 h post-treatment ( Figure 9A). Under hot stress, the expression of GhYTH genes peaked at 1 h and then decreased except for GhYTH6/12/14/16. However, the expression trend of GhYTH6/14/16 continuously increased and peaked at 24 h post treatment ( Figure 9B). Under drought stress, GhYTH genes exhibited higher expression in the first 6 h after treatment and then showed decreased expression except for GhYTH1/4/5/11/13/14/15/20. The expression level of GhYTH2/6/8/9/10/12/16/19 gradually increased post-drought treatment, indicating their possible roles in response to drought stress. (Figure 9C). Under salt stress, GhYTH6/16/20 were highly expressed at 1 h post-treatment ( Figure 9D). In brief, our data showed that the expression level of most YTH genes was induced in response to multiple stresses.

Subcellular Localization of GhYTH8 and GhYTH16
To investigate the subcellular localization of the GhYTH proteins, we fused the coding sequences of two of the GhYTH proteins (GhYTH8 and GhYTH16) that may be involved in the drought stress response to the GFP vector and transiently expressed in expanded leaves of N. benthamiana. Transient expression assays revealed that GFP signals in the GFP-YTH8 and GFP-GhYTH16 expressed in N. benthamiana leaves were observed in the whole cells, similar to the pattern observed for GFP alone ( Figure 10). These results suggested that the GhYTH8 and GhYTH16 proteins do not localize in specific compartments and are likely targeted ubiquitously in cells, which indicates their functional diversity and complexity. Figure 10. Subcellular analysis of GhYTH8 and GhYTH16 proteins. The GhYTH8::GFP, GhYTH16::GFP, and GFP alone were transiently expressed via agroinfiltration using N. benthamiana leaves. Green fluorescence was observed with confocal laser microscopy (left). The same cells were also observed by transmission microscopy, and both images were merged (from the middle to the right).

Subcellular Localization of GhYTH8 and GhYTH16
To investigate the subcellular localization of the GhYTH proteins, we fused the coding sequences of two of the GhYTH proteins (GhYTH8 and GhYTH16) that may be involved in the drought stress response to the GFP vector and transiently expressed in expanded leaves of N. benthamiana. Transient expression assays revealed that GFP signals in the GFP-YTH8 and GFP-GhYTH16 expressed in N. benthamiana leaves were observed in the whole cells, similar to the pattern observed for GFP alone ( Figure 10). These results suggested that the GhYTH8 and GhYTH16 proteins do not localize in specific compartments and are likely targeted ubiquitously in cells, which indicates their functional diversity and complexity.

Subcellular Localization of GhYTH8 and GhYTH16
To investigate the subcellular localization of the GhYTH proteins, we fused the coding sequences of two of the GhYTH proteins (GhYTH8 and GhYTH16) that may be involved in the drought stress response to the GFP vector and transiently expressed in expanded leaves of N. benthamiana. Transient expression assays revealed that GFP signals in the GFP-YTH8 and GFP-GhYTH16 expressed in N. benthamiana leaves were observed in the whole cells, similar to the pattern observed for GFP alone ( Figure 10). These results suggested that the GhYTH8 and GhYTH16 proteins do not localize in specific compartments and are likely targeted ubiquitously in cells, which indicates their functional diversity and complexity. Figure 10. Subcellular analysis of GhYTH8 and GhYTH16 proteins. The GhYTH8::GFP, GhYTH16::GFP, and GFP alone were transiently expressed via agroinfiltration using N. benthamiana leaves. Green fluorescence was observed with confocal laser microscopy (left). The same cells were also observed by transmission microscopy, and both images were merged (from the middle to the right). Figure 10. Subcellular analysis of GhYTH8 and GhYTH16 proteins. The GhYTH8::GFP, GhYTH16::GFP, and GFP alone were transiently expressed via agroinfiltration using N. benthamiana leaves. Green fluorescence was observed with confocal laser microscopy (left). The same cells were also observed by transmission microscopy, and both images were merged (from the middle to the right).

Functional Verification of GhYTH8 in Cotton under Drought Stress
To further explore the role of YTH genes that may be involved in the drought stress response according to their expression levels under drought stress, we performed a VIGS assay based on the full-length sequence of GhYTH8 in TM-1. The success of the VIGS experiment was verified by the albinism of positive control leaves and the detection of silencing efficiency. When the cotton grew to the three-leaf stage, the cotton infected by pYL156:00 (TRV:00) was used as the control, and the cotton plants were kept without water to simulate drought conditions. After 14 days of drought treatment, we can observe significant differences of the phenotypes of leaves. Compared with the control, the leaves of pYL156:GhYTH8 (TRV:GhYTH8)-infected plants under drought conditions wilted obviously ( Figure 11A). The expression level of GhYTH8 was significantly lower in the silenced plants than in the control (TRV:00) ( Figure 11B). Our data demonstrated that GhYTH8-silenced plants were more sensitive to drought stress, suggesting the functional roles of YTH proteins in the stress response of cotton.

Functional Verification of GhYTH8 in Cotton under Drought Stress
To further explore the role of YTH genes that may be involved in the drought stress response according to their expression levels under drought stress, we performed a VIGS assay based on the full-length sequence of GhYTH8 in TM-1. The success of the VIGS experiment was verified by the albinism of positive control leaves and the detection of silencing efficiency. When the cotton grew to the three-leaf stage, the cotton infected by pYL156:00 (TRV:00) was used as the control, and the cotton plants were kept without water to simulate drought conditions. After 14 days of drought treatment, we can observe significant differences of the phenotypes of leaves. Compared with the control, the leaves of pYL156:GhYTH8 (TRV:GhYTH8)-infected plants under drought conditions wilted obviously ( Figure 11A). The expression level of GhYTH8 was significantly lower in the silenced plants than in the control (TRV:00) ( Figure 11B). Our data demonstrated that GhYTH8-silenced plants were more sensitive to drought stress, suggesting the functional roles of YTH proteins in the stress response of cotton.

Discussion
The YTH domain-containing RBP has been recognized in eukaryotes for more than 20 years. Nevertheless, polyploidization events increase the complexity and make it harder to recognize the role of a gene family. Cotton is an important economic crop around the world, which provides a natural fiber source for the textile industry. Genome and transcriptome sequencing have facilitated the understanding of gene families in cotton. However, the identification of YTH genes and their biological functions in cotton remained elusive. In this present study, a total of 64 YTH genes (10 GaYTHs, 11 GrYTHs, 22 GbYTHs, and 21 GhYTHs) were identified and characterized in four cotton species at the genome-wide scale. The genome size of G. arboreum (1746 Mb) is almost twice of G. raimondii (885 Mb) [24]. The genome size of G. barbadense (~2.22 Gb) is almost the same as that of G. hirsutum (~2.30 Gb) [25]. Previous studies reported that 13,12,9,15, and 5 members of the YTH gene family were identified in Arabidopsis, rice, tomato, apple, and cucumber, respectively [21]. Our results showed that the number of YTH genes is not correlated with the size of the genomes.
The m 6 A modification is essential for modulating gene expression at the post-transcriptional level. To perform its regulatory function, the m 6 A modification complex needs to recruit reader proteins [26]. YTH family members have been shown to

Discussion
The YTH domain-containing RBP has been recognized in eukaryotes for more than 20 years. Nevertheless, polyploidization events increase the complexity and make it harder to recognize the role of a gene family. Cotton is an important economic crop around the world, which provides a natural fiber source for the textile industry. Genome and transcriptome sequencing have facilitated the understanding of gene families in cotton. However, the identification of YTH genes and their biological functions in cotton remained elusive. In this present study, a total of 64 YTH genes (10 GaYTHs, 11 GrYTHs, 22 GbYTHs, and 21 GhYTHs) were identified and characterized in four cotton species at the genomewide scale. The genome size of G. arboreum (1746 Mb) is almost twice of G. raimondii (885 Mb) [24]. The genome size of G. barbadense (~2.22 Gb) is almost the same as that of G. hirsutum (~2.30 Gb) [25]. Previous studies reported that 13,12,9,15, and 5 members of the YTH gene family were identified in Arabidopsis, rice, tomato, apple, and cucumber, respectively [21]. Our results showed that the number of YTH genes is not correlated with the size of the genomes.
The m 6 A modification is essential for modulating gene expression at the post-transcriptional level. To perform its regulatory function, the m 6 A modification complex needs to recruit reader proteins [26]. YTH family members have been shown to act as m 6 A reader proteins [27]. In humans, five YTH proteins, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3, were identified. [28]. Studies have shown that a conserved mechanism was applied by YTHDFs and YTHDCs to identify m 6 A [29]. In this study, 64 YTH genes were clustered into three groups. Arabidopsis YTH genes AtYTH7/9/13/5 (also named AtDF1A-AtDF4A), AtYTH10/4/6/1 (also named AtDF1C-AtDF4C), AtYTH8/2/12 (also named AtDF1B-AtDF1B) were clustered into Groups I-III, respectively. It is worth noting that although AtYTH3 (AtDC1A) and AtYTH11 (AtDC1B) were grouped into Group III, no other cotton YTH genes were in the same branch with these two Arabidopsis YTHDC genes. Furthermore, YTH1 superfamily domains were identified in all three common wheat YTHDC proteins [16]. However, no YTH1 superfamily domains have been found in cotton YTH proteins in this study, suggesting no YTHDC gene exists in cotton species. G. barbadense and G. hirsutum are allotetraploids that originated from the hybridization among A-genome-like ancestral African species and D-genome-like American species [30]. Ancestral diploids, G. arboreum and G. raimondii, contributed to the At and Dt subgenomes for both G. hirsutum and G. barbadense, respectively, through genome-wide duplication events [25,31]. YTHs located on At subgenomes in both G. barbadense (GbYTH1-GbYTH11) and G. hirsutum (GhYTH1-GhYTH10) tend to be clustered with YTHs in G. arboreum (A genome). Similarly, YTHs located on Dt subgenomes in both G. barbadense (GbYTH12-GbYTH22) and G. hirsutum (GhYTH11-GhYTH21) tend to be clustered with YTHs in G. raimondii (D genome). Our findings also confirmed the sources of At subgenomes and Dt subgenomes in two tetraploid species. YTH genes of G. hirsutum and G. barbadense tend to be clustered in the same branch. Our results indicated that the genetic relationship between two allotetraploids was closer than that between allotetraploids and diploids. It might be due to the coevolution of G. hirsutum and G. barbadense.
GaYTH and GrYTH genes are dispersed across 7 out of 13 chromosomes. Moreover, GbYTH and GhYTH genes are distributed on 7 out of 13 chromosomes across both subgenomes. YTH genes showed an uneven distribution on cotton chromosomes. This chromosomal distribution pattern of YTH genes was also found in Arabidopsis, rice, and cucumber [13,21]. It is worth noting that GbYTH8 anchored on chromosome A10 in G. barbadense was not identified in G. hirsutum, indicating that duplication and deletion events may appear during evolution.
Polyploidy is one of the major mechanisms of plant formation and environmental adaptation [32]. Gene duplications are regarded as one of the key driving forces during the evolution of genomes and genetic systems. Duplicated genes offer raw material for the formation of novel genes, which may facilitate the generation of novel functions [33]. Tandem duplication, segmental duplication, and transposition events are three principal evolutionary patterns, of which tandem duplications and segmental duplications were found to be the main causes of gene family expansion in plants [34,35]. Furthermore, segmental duplications have been identified to play important roles in expanding the gene family members in plants [36,37]. In this study, no tandem duplication was detected in the YTH gene family, though the genes may have undergone independent evolutionary processes in two allotetraploids cotton according to the synteny analysis.
Environmental stress responses involve various physiological, cellular, and molecular adaptations. Plants can produce and accumulate phytohormones, such as abscisic acid (ABA), Jasmonic acid (JA), ethylene (ET), and salicylic acid (SA), in response to stresses [38]. Under stressed conditions, molecular level changes in plants are principally caused by transcription factors binding to the cis-regulatory elements upstream to the stress-responsive genes [39]. The pattern of existence of these cis-elements reveals the mechanism of stressresponsive upregulation of downstream genes [40]. In this study, many inducible elements, such as ABA, ET, and SA responsive elements were recognized in the promoter region of YTH genes, suggesting functions of cotton YTH genes in stress responses. Further research is required to decipher their roles in stress responses by recognizing their gene structures, protein structures, and target RNAs.
In general, the abundant or increased expression of a gene in a tissue, a developmental stage, or a stress condition may suggest its roles correlated to developmental and stress responses. Publicly available transcriptomic data were analyzed to explore the possible function of the GhYTH genes in development and stress responses [41]. The production and quality of cotton are closely related to the development of ovule and fiber. In our study, a set of GhYTH genes were expressed at higher levels during the development of the ovule and fiber of cotton. These highly expressed YTH genes in wheat spikes might function in the development and architecture of spikes [16]. Furthermore, the expressions of the GhYTH gene among different tissues were different, suggesting that they may play numerous roles in different cotton tissues. The YTH genes are known to play essential roles in stress responses in several plant species. In this present study, expressions of some GhYTH genes were induced after stress treatments. Similar expression patterns were also observed in Arabidopsis, rice, cucumber, and apple [14,20,21,42].
sRNAs (small RNAs) ranging from 21 to 24 nt have obtained much attention for their biological roles in plant growth, development, and biotic and abiotic stress responses [43]. MiRNAs are a sort of sRNA, which could suppress the expression of target genes via binding to the mRNAs [44]. Drought, one of the most severe abiotic stresses, limits crop growth and yield by changing metabolic activity and biological functions [45]. MiRNAs have been found to play important roles in the drought stress response in cotton. In this study, 18 miRNA targets of GhYTH genes were detected. GhYTH3 and GhYTH13 could be targeted by ghr-miR156a, ghr-miR156b, and ghr-miR156d. One of these targets, ghr-miR156a was expressed differently between control and drought treatment in a previous report [46]. In addition, the miR156 family was identified to serve as the main regulator in response to high-temperature stress during G. hirsutum anther development [47]. These predicted miRNA targets could provide valuable candidates for the experimental study of the stress response of GhYTH genes.
In humans, YTH domain-containing RNA-binding proteins were extensively studied [48]. YTHDF1/2/3 could bind to N6-methyladenosine and thus control the stability of mRNA [27]. YT521 (YTHDC1) is characterized by alternately spliced isoforms with regulatory impact on cancer [49]. In contrast, the functional research of the YTH proteins in plants is very limited. VIGS was widely used to explore the gene functions in cotton for the difficulties inherent in cotton transformation [50]. The effect of GhYTH silencing by VIGS was evaluated, and the results showed that the silencing of GhYTH8 decreased the resistance of the plants to droughts. Our findings verified the GhYTH had a vital function in the resistance of G. hirsutum to droughts.

Identification and Physiochemical Properties Analysis of YTH Gene Family Members in Four Gossypium Species
The gene annotation and genome files of two diploid species, G. arboreum L. (CRI) and G. raimondii L. (JGI), and two tetraploid species, G. hirsutum L. (CRI) and G. barbadense L. (ZJU), were retrieved from the Cotton Functional Genomics Database (https://cottonfgd. org/(accessed on 11 January 2023)) [41]. The protein sequences of YTH genes in Arabidopsis were obtained from the TAIR10 database (https://www.arabidopsis.org/ (accessed on 11 January 2023)) [51]. The query sequences of candidate YTH genes against the protein sequences of four cotton species were aligned by BLASTp. The candidate genes were analyzed in the NCBI-CDD (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb. cgi (accessed on 11 January 2023)) database to identify the YTH domain. The identified YTH genes were further confirmed by HMM search function in HMMER3.0 [52]. The physiochemical properties, relative molecular mass (MW), and theoretical isoelectric point (pI) were determined through the ExPASy proteomics server (https://web.expasy.org/ compute_pi/ (accessed on 11 January 2023)) [53]. The genomic locations of YTH genes in the chromosome were extracted from the annotation files and were visualized by MapChart (https://www.wur.nl/en/show/Mapchart/ (accessed on 11 January 2023)).

Evolutionary, Gene Structure, and Conserved Motif Analysis of YTH Genes
All the identified YTH protein sequences from four cotton species and Arabidopsis were aligned using MUSCLE algorithm program in MEGA-X [54]. The maximum likelihood (ML) method of MEGA-X was exploited to construct the phylogenetic tree with the LG + G model and 1000 bootstrap replications, followed by embellishing the phylogenetic tree with Evolview (https://evolgenius.info/ (accessed on 11 January 2023)) [55].
The conserved motifs of cotton YTH protein sequences were analyzed using the MEME program (http://meme-suite.org/ (accessed on 11 January 2023)) [56]. The gene structure information of the cotton YTH genes was retrieved from Cotton Functional Genomics Database, and the gene structures were graphically visualized using TBtools [57].

Chromosomal Distribution and Synteny Analysis of YTH Genes
CottonFGD genome annotation files were used to determine the chromosomal distribution of YTH genes in 4 cotton species. Duplication events and syntenic blocks of cotton YTH genes were found by McScanX [58], followed by visualizing the orthologous YTH genes between four cotton species through CIRCOS. Ka/Ks analysis was performed as previously described [59].

Regulatory Elements Analysis of YTH Genes
The 2000 bp upstream region of coding sequence (CDS) of YTH genes extracted using TBtools was regarded as the promoter region. The cis-acting regulatory elements of YTH genes were identified using PlantCare (http://bioinformatics.psb.ugent.be/webtools/ plantcare/html/ (accessed on 11 January 2023)) [60] and then visualized by R.

Virus-induced Gene Silencing (VIGS) Experiment and Drought Treatment
We performed VIGS experiments to further confirm whether YTH genes affect drought stress tolerance in cotton. Total RNA was extracted from the leaves and then reverse transcribed to cDNA. A 313 bp fragment of GhYTH8 ORF was amplified and cloned into the pYL156 to construct the fusion vector pYL156:GhYTH8. The vectors (pYL156:GhYTH8 and pYL156:CLA1) were transformed into A. tumefaciens strain GV3101, followed by cultivating overnight in L.B. medium with kanamycin 50 mg/L and rifampicin 25 mg/L at 30°C. pYL156:00 (TRV:00), pYL156:CLA1, and pYL156:GhYTH8 were mixed with a 1:1 ratio, followed by agroinfiltrating into fully unfolded cotton cotyledons with a needlefree syringe. The standard variety TM-1 (G. hirsutum) was used as the experimental plant. After a night of darkness, the plants were transferred to a greenhouse at 20 • C with a 16 h light/8 h dark photoperiod. When TRV:CLA1 plants showed phenotype, TRV:GhYTH8 and TRV:00 plants were kept without water to simulate drought conditions.

Quantitative Real-Time PCR (qPCR)
Three sets of samples from different plants were collected as three biological replicates. The leaves of TRV:00 and TRV:GhYTH8 cotton seedlings were used for total RNA extraction with an EASYspin Plus Plant RNA Kit (Aidlab, Beijing, China) according to the manufacturer's instructions. The RNA reverse-transcribed used Trans Script One-Step gDNA Removal and cDNA Synthesis SuperMix (Trans Gen, Beijing, China). Cotton Histon3 (GhHiston3) was used as an endogenous standard control. qPCR assays were performed on a QuantStudio 6 Flex thermocycler (Applied Biosystems, Foster City, CA, USA) using Trans Start Top Green qPCR SuperMix (TransGen, Beijing, China) and repeated three times.

Conclusions
Above all, our results highlight the potential role of the YTH genes that are involved in the drought stress response. Our results provide new insights to further study the biological roles of the YTH genes in stress response.