Genome-Wide Investigation of the MiR166 Family Provides New Insights into Its Involvement in the Drought Stress Responses of Tea Plants ( Camellia sinensis (L.) O. Kuntze)

: MicroRNA166 (miR166) is a highly conserved plant miRNA that plays a crucial role in plant growth and the resistance to various abiotic stresses. However, the miR166s in tea ( Camellia sinensis (L.) O. Kuntze) have not been comprehensively identiﬁed and analyzed. This study identiﬁed 30 mature miR166s and twelve pre-miR166s in tea plants. An evolutionary analysis revealed that csn-miR166s originating from the 3 (cid:48) arm of their precursors were more conserved than the csn-miR166s derived from the 5 (cid:48) arm of their precursors. The twelve pre-miR166s in tea were divided into two groups, with csn-MIR166 Scaffold364-2 separated from the other precursors. The Mfold-based predictions indicated that the twelve csn-MIR166s formed typical and stable structures comprising a stem-loop hairpin, with minimum free energy ranging from − 110.90 to − 71.80 kcal/mol. An analysis of the CsMIR166 promoters detected diverse cis -acting elements, including those related to light responses, biosynthesis and metabolism, abiotic stress defenses, and hormone responses. There was no one-to-one relationship between the csn-miR166s and their targets, but most csn-miR166s targeted HD-Zip III genes. Physiological characterization of tea plants under drought stress showed that leaf water content proportionally decreased with the aggravation of drought stress. In contrast, tea leaves’ malondialdehyde (MDA) content proportionally increased. Moreover, the cleavage site of the ATHB-15-like transcript was identiﬁed according to a modiﬁed 5 (cid:48) RNA ligase-mediated rapid ampliﬁcation of cDNA ends. The RT-qPCR data indicated that the transcription of nine csn-miR166s was negatively correlated with their target gene. crucial the role of the regulatory mechanism to


Introduction
Tea (Camellia sinensis (L.) O. Kuntze) is an important economic crop that originated in southwestern China. The environmental stresses associated with the frequent occurrence of extreme weather events may adversely affect tea plant growth and development. Drought stress, which is one of the major challenges faced by growers in the major tea-producing regions, seriously affects tea plant growth and development and tea yield and quality [1,2]. To cope with the drought stress, tea plants would use a series of mechanisms, including morphological, physiological (e.g., cellular membrane damage, photosynthesis changes, osmotic regulation, etc.) and molecular levels to adapt to changes in environmental [3][4][5]. Therefore, the mechanism mediating tea plant drought resistance should be thoroughly characterized.
14.37%, 9.29%, and 5.34% under CK, T1, T2, and T3, respectively, and treatment lasted for ten days. The soil water holding capacity was 23.45%. Each treatment was analyzed by collecting three biological replicates, and each replicates in at least ten randomly selected leaves. Tender leaves randomly distributed in the canopies were collected and immediately frozen in liquid nitrogen and kept at −80 • C until analyzed.

Identification of Csn-MiR166s in Tea Plant's Genome and Acquisition of Plant MiR166 Sequences
The csn-miR166 and pre-miR166 sequences were obtained based on small RNA sequencing results and previous reports [37][38][39]. The obtained sequences were named by combining the abbreviation for Camellia sinensis (L.) O. Kuntze (csn) and the common name for their highly homologous sequences and the pre-miR166 positions on scaffolds (e.g., csn-miR166a, b, and c and csn-MIR166 Scaffold1344, csn-MIR166 Scaffold1761, and csn-MIR166 Scaffold2148) (Additional file 1. Table S1 and Additional file 2. Table S2). As one of the highly conserved miRNAs, miR166s have been identified in 45 plant species (including monocotyledonous species and dicotyledonous species). To comprehensively understand the sequence and evolution characteristics of miR166 and pre-miR166, the miR166 and pre-miR166 sequences from Arabidopsis thaliana (L.) Heynh, O. sativa, Zea mays L., Vitis vinifera L., and Malus domestica B. were downloaded from the miRBase database (release 22.1) (http://www.mirbase.org/; accessed on 13 November 2020).

Bioinformatic Analyses of Csn-MiR166
The csn-miR166 and pre-miR166 sequences were aligned using DNAMAN (version 6.0, San Ramon, CA, USA). A phylogenetic tree was subsequently constructed according to the neighbor-joining algorithm (1000 bootstrap replicates) using the MEGA (version 7.0, Auckland, New Zealand). The default parameters of the Mfold 3.5 tool were used to assess the secondary structure.
Based on the genome sequences of 'Shuchazao' [29] and 'Yunkang 10 cultivars [27], the csn-miR166 targets were predicted using the psRNATarget (http://plantgrn.noble.org/ psRNATarget/?function=3; accessed on 26 April 2020) online tool. The default parameters were applied, with the exception of the expected value, which was set to 3.0 [25]. Additionally, the Cytoscape (version 3.9.1, Boston, MA, USA) analyzed the potential regulatory relationship of csn-miR166s and their targets in tea plants.

Physiological Characterization of the Tea Plant under Different Drought Degrees
The soil water content was determined according to the method for determining soil water content (NY/T 52-1987). In brief, 20 g of soil samples were taken from the potting and oven-dried at 105 ± 2 • C for 12 h. Once the soil samples have cooled to room temperature, it is weighed. Soil water contents were determined as the difference between the fresh and dry soil mass, expressed as a percentage. The leaf water content was determined according to the method for tea-determination of moisture content (GB/T 8304-2002). The Malondialdehyde (MDA) content of tea leaves was measured using the MDA Assay Kit (Solabol Technology Co., Ltd., Beijing, China).

Total RNA Isolation and cDNA Synthesis
Total RNA was extracted from the samples mentioned above using Transzol UP (TransGen Biotech, Beijing, China). The quality and concentration of the extracted RNA were respectively determined with a 1% agarose gel electrophoresis and the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The ratios of OD 260 /oD 280 between 1.90 and 2.10 indicated the samples were of sufficient quality for reverse transcribed to cDNA for target genes of miR166s and first-strand cDNA of miRNA using the TransScript ® Uni One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen) and the TranScript ® miRNA First-Strand cDNA Synthesis SuperMix (TransGen) kits, respectively.

Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) Analyses
The Primer3 (https://bioinfo.ut.ee/primer3/; accessed on 9 January 2021) online tool was used to design RT-qPCR primers specific for the csn-miR166s, potential targets, and reference genes (Additional file 3. Table S3). The RT-qPCR analysis of the csn-miR166s and their targets were conducted using the TransStart ® Tip Green qPCR SuperMix (TransGen) and the LightCycler 480 system (Roche Applied Sciences, Basel, Switzerland). Relative expression levels were calculated according to the 2 −∆∆Ct method [40]. Each reaction was performed independently in triplicate. Additionally, U6 small nuclear RNA (U6 snRNA) and SAND were selected as the reference genes for the csn-miR166s and their targets, respectively. Equal amounts of the total RNA extracted from the samples mentioned above were mixed and then reverse transcribed for a modified RLM-RACE using the FirstChoice™ RLM-RACE Kit (Invitrogen, CA, USA) [41]. The target gene cleavage primers were designed using the Primer3 (https://bioinfo.ut.ee/primer3/; accessed on 13 November 2021) online tool (Additional file 3. Table S3). Two rounds of nested PCR amplification were performed using the gene-specific and universal primers. The amplified products were analyzed by 1.5% agarose gel electrophoresis. The target fragments were purified from the gel using the EasyPure ® Quick Gel Extraction Kit (TransGen) and cloned into the T1 vector (TransGen). Single colonies of putative transformants were selected for sequencing (BioSune, Fuzhou, China).

Statistical Analysis
All dates were provided herein as the mean ± standard deviation. The SPSS (version 25.0, Armonk, NY, USA) software was used for a one-way analysis of variance followed by Tukey's post-hoc test. Differences were considered significant at p < 0.05.

Analysis of MiR166 Sequences in Plants
To investigate the evolutionary relationships among miR166s, we used the miR166s from C. sinensis, A. thaliana, O. sativa, Z. mays, M. domestica, and V. vinifera to construct a phylogenetic tree (Figure 1a). All miR166s were classified into two groups according to whether they derived from the 5 arm or 3 arm of their precursors. Mutations, deletions, and insertions were detected in the miR166 sequences at nucleotide positions 1, 2, 4, 6, 17, 19, 20, and 21; the remaining nucleotides were highly conserved ( Figure 1b). Further analyses revealed fewer mutations, deletions, and insertions at the 3 end than at the 5 end of the csn-miR166 precursors.

Analysis of Cis-Acting Elements in CsMIR166 Promoter Regions
To identify the cis-acting elements in CsMIR166 promoter sequences, the TSSP tool of Softberry was used to predict the TSS of CsMIR166s and the 1500-bp region upstream of the TSS of CsMIR166 genes (i.e., hypothetical promoter region) was analyzed using the

Analysis of Cis-Acting Elements in CsMIR166 Promoter Regions
To identify the cis-acting elements in CsMIR166 promoter sequences, the TSSP tool of Softberry was used to predict the TSS of CsMIR166s and the 1500-bp region upstream of the TSS of CsMIR166 genes (i.e., hypothetical promoter region) was analyzed using the PlantCARE server. The core promoter element TATA-box and the common CAAT-box element were detected in the promoter and enhancer regions of all CsMIR166 genes. The other identified elements in CsMIR166 promoters were divided into the following five categories: light response-, biosynthesis and metabolism-, abiotic defense-, and hormone response-related cis-acting elements ( Figure 5). The cis-acting elements were unequally distributed among the CsMIR166 genes. The light-responsive elements were the most abundant (all promoters contained photoresponsive elements), followed by the anaerobic responsive elements (66.7%). In contrast, some cis-acting elements, such as meristematic tissue expression elements (CAT-box, 33.3%), endosperm expression elements (GCN4, 25.0%), drought response elements (MBS, 16.7%), gibberellin response elements (GARE and P-box, 16.7%), and abscisic acid response elements (ABRE, 16.7%), were identified in no more than three CsMIR166 genes.

Prediction of csn-miR166 targets
To explore the biological functions of csn-miR166s, based on the tea reference genome, potential target genes of csn-miR166s were predicted using the psRNATarget tool, which resulted in a total of 151 potential target genes were predicted (Figure 6a and Additional file 5. Table S5). The regulatory relationships between the csn-miR166s and their characteristic target genes revealed a lack of a one-to-one relationship between the regulation of csn-miR166s and their respective targets. The HD-Zip III family, the primary gene family controlled by miR166s, was simultaneously targeted by fifteen csn-miR166s. The targets included one ATHB-14-like (CSS049178.1), two REVOLUTA-like (CSS017546.1 and CSS042067.1), two ATHB-15 (CSS002611.1 and CSS043540.1), and two ATHB-15-like (CSS012922.1 and CSS044030.1) genes (Figure 6b and Additional file 6. Table S6). Further analyses indicated that csn-miR166s from the 3′ arm of their precursors could target HD-Zip III genes, whereas csn-miR166s from the 5′ arm of their precursors could target genes that did not belong to the HD-Zip III family.

Prediction of csn-miR166 targets
To explore the biological functions of csn-miR166s, based on the tea reference genome, potential target genes of csn-miR166s were predicted using the psRNATarget tool, which resulted in a total of 151 potential target genes were predicted (Figure 6a and Additional file 5. Table S5). The regulatory relationships between the csn-miR166s and their characteristic target genes revealed a lack of a one-to-one relationship between the regulation of csn-miR166s and their respective targets. The HD-Zip III family, the primary gene family controlled by miR166s, was simultaneously targeted by fifteen csn-miR166s. The targets included one ATHB-14-like (CSS049178.1), two REVOLUTA-like (CSS017546.1 and CSS042067.1), two ATHB-15 (CSS002611.1 and CSS043540.1), and two ATHB-15-like (CSS012922.1 and CSS044030.1) genes (Figure 6b and Additional file 6. Table S6). Further analyses indicated that csn-miR166s from the 3 arm of their precursors could target HD-Zip III genes, whereas csn-miR166s from the 5 arm of their precursors could target genes that did not belong to the HD-Zip III family.
characteristic target genes revealed a lack of a one-to-one relationship between the re lation of csn-miR166s and their respective targets. The HD-Zip III family, the primary g family controlled by miR166s, was simultaneously targeted by fifteen csn-miR166s. targets included one ATHB-14-like (CSS049178.1), two REVOLUTA-like (CSS017546.1 CSS042067.1), two ATHB-15 (CSS002611.1 and CSS043540.1), and two ATHB-15 (CSS012922.1 and CSS044030.1) genes (Figure 6b and Additional file 6. Table S6). Fur analyses indicated that csn-miR166s from the 3′ arm of their precursors could target Zip III genes, whereas csn-miR166s from the 5′ arm of their precursors could target ge that did not belong to the HD-Zip III family.

Analysis of Tea Plants Relative Water Content and Malondialdehyde Content under Drought Stress
To explore the physiological changes in tea plants under the different drought degrees, the relative water content (soil water content and leaf water content) and malondialdehyde content were determined (Figure 7). With the aggravation of drought stress, the soil water content and leaf water content proportionally decreased in T1, T2, and T3 with respect to the normal water supply (CK), whereas the malondialdehyde content proportionally increased.

Analysis of Tea Plants Relative Water Content and Malondialdehyde Content under Drought Stress
To explore the physiological changes in tea plants under the different drought degrees, the relative water content (soil water content and leaf water content) and malondialdehyde content were determined (Figure 7). With the aggravation of drought stress, the soil water content and leaf water content proportionally decreased in T1, T2, and T3 with respect to the normal water supply (CK), whereas the malondialdehyde content proportionally increased.
Subsequently, the cleavage site of the ATHB-15-like transcript was identified based on a modified 5 RLM-RACE analysis (Figure 9). The transcription of ATHB-15-like was regulated via the cleavage of the binding region between nucleotides 12 and 13 from the 5 end of the csn-miR166s. Hence, miR166s can regulate the HD-Zip III family by modulating the transcription of related target genes following exposure to drought stress. Subsequently, the cleavage site of the ATHB-15-like transcript was identified on a modified 5′ RLM-RACE analysis (Figure 9). The transcription of ATHB-15-lik regulated via the cleavage of the binding region between nucleotides 12 and 13 fro 5′ end of the csn-miR166s. Hence, miR166s can regulate the HD-Zip III family by m lating the transcription of related target genes following exposure to drought stress

Evolutionary Characteristics of Csn-MiR166s
The miR166s are highly conserved miRNAs that are widely distributed in 45 species, including angiosperms, gymnosperms, ferns, and mosses [25,26,43]. The mation available in the miRBase database indicates that miR166s are abundant in plant species, including Glycine max (L.) Merr (26), Z. mays (26), O. sativa (24), Brachyp distachyon (L.) Beauv (18), and Populus trichocarpa (Torr & Gray) (17). Based on the miR database, the G. max and Z. mays had the highest number of miR166 members a plants, each containing 26 miR166 members. In the present study, 30 csn-miR166s identified, more than the number of miR166 in G. max and Z. mays to our knowl Increases in the number of csn-miR166 genes may have been facilitated by whole-ge duplication events [30] in the large tea genome (approximately 3 Gb) [27][28][29]33]. tionally, conserved miRNAs reportedly underwent both conserved and diverse e tionary processes [44]. We identified 30 csn-miR166s and divided them into two categ ( Figure 1). The csn-miR166s derived from the precursor 3′ arm were more conserved the csn-miR166s originated from the precursor 5′ arm (Figure 1b). The more cons csn-miR166s had the same putative target genes in the HD-ZIP III family (Figure 6 contrast, the csn-miR166s derived from the precursor 5′ arm targeted many other encoding functionally diverse proteins, including ABC transporter and proteins re sive to stress and auxin (Additional file 5. Table S5). A previous investigation prove miR166s are highly complex miRNAs regarding their evolution. There is no one-t relationship between the number of miR166s in a species and their corresponding Different mature sequences may have originated from the same or multiple loci [26 example, the miR166s of G. max (26), Z. mays (26), O. sativa (24), B. distachyon (18), a trichocarpa (17) are respectively associated with 21, 14, 13, 10, and 13 genomic loci. I 30 csn-miR166s were revealed to be associated with 12 loci. Interestingly, of the 3 miR166s, some csn-miR166s (e.g., miR166b-2, b-3, d-1, e-1-3p, a-4-5p, b-5p, d-5p, g

Evolutionary Characteristics of Csn-MiR166s
The miR166s are highly conserved miRNAs that are widely distributed in 45 plant species, including angiosperms, gymnosperms, ferns, and mosses [25,26,43]. The information available in the miRBase database indicates that miR166s are abundant in many plant species, including Glycine max (L.) Merr (26), Z. mays (26), O. sativa (24), Brachypodium distachyon (L.) Beauv (18), and Populus trichocarpa (Torr & Gray) (17). Based on the miRBase database, the G. max and Z. mays had the highest number of miR166 members among plants, each containing 26 miR166 members. In the present study, 30 csn-miR166s were identified, more than the number of miR166 in G. max and Z. mays to our knowledge. Increases in the number of csn-miR166 genes may have been facilitated by whole-genome duplication events [30] in the large tea genome (approximately 3 Gb) [27][28][29]33]. Additionally, conserved miRNAs reportedly underwent both conserved and diverse evolutionary processes [44]. We identified 30 csn-miR166s and divided them into two categories (Figure 1). The csn-miR166s derived from the precursor 3 arm were more conserved than the csn-miR166s originated from the precursor 5 arm (Figure 1b). The more conserved csn-miR166s had the same putative target genes in the HD-ZIP III family (Figure 6b). In contrast, the csn-miR166s derived from the precursor 5 arm targeted many other genes encoding functionally diverse proteins, including ABC transporter and proteins responsive to stress and auxin (Additional file 5. Table S5). A previous investigation proved that miR166s are highly complex miRNAs regarding their evolution. There is no one-to-one relationship between the number of miR166s in a species and their corresponding loci. Different mature sequences may have originated from the same or multiple loci [26]. For example, the miR166s of G. max (26), Z. mays (26), O. sativa (24), B. distachyon (18), and P. trichocarpa (17) are respectively associated with 21, 14, 13, 10, and 13 genomic loci. In tea, 30 csn-miR166s were revealed to be associated with 12 loci. Interestingly, of the 30 csn-miR166s, some csn-miR166s (e.g., miR166b-2, b-3, d-1, e-1-3p, a-4-5p, b-5p, d-5p, g-1-5p, g-2-5p, and h-5p) did not precisely match the precursor sequences, suggesting that the csn-miR166 family has been affected by mutations during its formation and may have more undiscovered or unconfirmed loci. This hypothesis needs to be confirmed by further experiments. Additionally, csn-MIR166 scaffold364-2 was separated far from eleven pre-miR166s (Figure 4), suggesting that tea pre-miR166s may have complex and distinct evolutionary rates and diverse modes of evolution. Based on that evidence, we speculated that csn-miR166s are both conservation and diversification in evolution and function.

Cleavage of HD-ZIP III Transcripts by Csn-MiR166s Is a Crucial Mechanism Underlying Tea Plant Drought Tolerance
Examining the putative promoter regions is likely to provide clues regarding the csn-miR166 regulatory mechanism. Cis-regulatory elements are crucial for the transcriptional regulation of MIR genes [6]. We observed that putative CsMIR166 promoter regions include many cis-regulatory elements, including MBS and ABRE elements ( Figure 5). The miRNA genes (e.g., MIR167, MIR168, and MIR396) with promoter regions that contain ABRE elements are important for drought stress responses [45,46]. Similarly, miR166 genes with MBS and ABRE elements in their promoter regions might also play crucial roles during plant responses to drought stress. As expected, the transcription patterns of csn-miR166s are affected by drought stress (Figure 8a). Previous studies reported that different miR166 members possibly present a special manner or regulatory role in withstanding the abiotic stresses [2,25]. In this study, the transcription patterns of csn-miR166 members differed, and the regulation of csn-miR166s on their targets displayed diversity under the different drought degrees (Figure 8 and Additional file 7. Table S7). It is indicated that the regulation of csn-miR166 family members on their targets under the different drought degrees may be complex, with the division of roles.
MicroRNAs can regulate the expression of their target genes by cleaving the corresponding mRNA or inhibiting translation according to the perfect or near-perfect complementarity principle [25]. Cleavage is one of the main miRNA modes of action for silencing their target genes [6]. The HD-ZIP III genes, which are critical for plant responses to various abiotic stresses, are the main targets of miR166s [2,47,48]. In this study, the prediction of csn-miR166 targets showed its role in inhibiting the transcription of their putative targets (HD-ZIP III family members) primarily via transcript cleavage (Additional file 6. Table  S6). This was indeed the case since the csn-miR166s can cleave the ATHB-15-like transcript by applying a modified 5 RLM-RACE strategy ( Figure 9). This is consistent with the previous report that miR166 could respond to drought stress by cleaving HD-ZIP III family members [11,25]. Thus, we speculated that csn-miR166-mediated HD-ZIP III cleavage could be an important mechanism for mediating tea plant drought resistance. The HD-ZIP III genes encode transcription factors that regulate the expression of many downstream stress-associated genes [23]. In future studies, a systematic exploration of the specific mechanism by which csn-miR166s control the transcription of downstream genes through their regulatory effects on HD-ZIP III genes to protect tea plants from drought stress will help to improve our understanding of tea plant resistance to drought stress.

Conclusions
In this study, csn-miR166s and csn-pre-miR166s were analyzed to elucidate the characteristics and effects of csn-miR166s on tea plant drought resistance. We revealed the conserved and diverse evolution and functions of csn-miR166s. Moreover, the RT-qPCR analysis indicated that csn-miR166 family members in responding to drought stress might have the division of labor and complementary function. The csn-miR166-target gene pairs confirmed that csn-miR166s regulate the transcription of HD-Zip III family members. More specifically, they control target gene expression levels to mediate drought stress resistance. Our study provides crucial insights into the role of the csn-miR166-mediated regulatory mechanism contributing to tea plant resistance to drought stress.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/f13040628/s1, Table S1. The sequences of csn-miR166, pre-miR166, and the target genes of csn-miR166s, Table S2. Naming information of csn-miR166s and pre-miR166s, Table S3. Primers for RT-qPCR, Table S4. The location information of the pre-miR166s scaffold, Table S5. Prediction of the csn-miR166 target gene, Table S6. Naming information of csn-miR166 target gene, Table S7. The specific functions of different csn-miR166 members under the different drought degrees.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.