Identification of miRNAs and Their Targets Involved in Flower and Fruit Development across Domesticated and Wild Capsicum Species

MicroRNAs (miRNAs) are regulators of the post-transcription stage of gene activity documented to play central roles in flower and fruit development in model plant species. However, little is known about their roles and differences in domesticated and wild Capsicum species. In this study, we used high-throughput sequencing to analyze the miRNA content at three developmental stages (flower, small fruit, and middle fruit) from two cultivated (C. baccatum and C. annuum) and two wild (C. chacoense and C. eximium) pepper species. This analysis revealed 22 known and 27 novel miRNAs differentially expressed across species and tissues. A number of stage- and species-specific miRNAs were identified, and Gene Ontology terms were assigned to 138 genes targeted by the miRNAs. Most Gene Ontology terms were for the categories “genetic information processing”, “signaling and cellular processes”, “amino acid metabolism”, and “carbohydrate metabolism”. Enriched KEGG analysis revealed the pathways amino acids, sugar and nucleotide metabolism, starch and sucrose metabolism, and fructose-mannose metabolism among the principal ones regulated by miRNAs during pepper fruit ripening. We predicted miRNA–target gene interactions regulating flowering time and fruit development, including miR156/157 with SPL genes, miR159 with GaMYB proteins, miR160 with ARF genes, miR172 with AP2-like transcription factors, and miR408 with CLAVATA1 gene across the different Capsicum species. In addition, novel miRNAs play an important role in regulating interactions potentially controlling plant pathogen defense and fruit quality via fructokinase, alpha-L-arabinofuranosidase, and aromatic and neutral amino acid transporter. Overall, the small RNA-sequencing results from this study represent valuable information that provides a solid foundation for uncovering the miRNA-mediated mechanisms of flower and fruit development between domesticated and wild Capsicum species.


Introduction
MicroRNAs (miRNAs) are a specific class of 20-to 24-nt endogenous small nonprotein coding RNAs involved in post-transcriptional and translational gene expression regulation in plants and animals [1,2]. Mature miRNA, generated from longer pri-RNA via nuclease cleavage processes [3], negatively regulate gene expression by recognition and complementary binding to the open reading frame or untranslated regions (UTRs) of target genes. The expression of multiple genes can be regulated by a single miRNA, and multiple miRNAs can control a single gene expression [4]. In plants, the complementarity between miRNA and their targets is very high, which results in RNA-induced silencing

Analysis of Small RNAs in Pepper
For comprehensive analysis of the Capsicum miRNAome and to identify putative miRNAs associated with fruit development, we sequenced 12 small-RNA libraries derived from flower, small fruit (6 days post anthesis) and medium fruit (25 days post anthesis) developmental stages of four Capsicum species (Figure 1) by using the Illumina HiSeq 2500 instrument. A total of 162,841,574 raw reads were obtained from the 12 libraries. After removing adaptors, low-quality reads, poly-A sequences, and reads of <18 and >30 nt in length, the number of clean reads ranged from~3 million in C. annuum small fruit, tõ 8 million in C. baccatum small fruit (Table S1). The small RNA clean reads were further classified into different categories by BLAST searches against Rfam and Repbase databases, and noncoding RNAs including rRNA, tRNA, snoRNA, snRNA sequences were removed. The remaining sequences were then mapped to the pepper genome to determine whether they could be candidate miRNAs and were selected based on strict criteria for annotation of plant miRNA. Furthermore, parameters such as length distribution and 5 -end of candidate miRNAs were considered by using the clean reads from all species and tissues ( Figure 2). Most of these small RNAs had a 5 -end terminal U or A, which indicates canonical small RNAs [24]. Likewise, the length distribution ranged from 20 to 24 nt, the typical length of canonical miRNAs. Although the size distribution of all small RNAs was diverse and varied across species and tissues, the dominant length was 24 nt followed by 23,22, and 21 nt in C. baccatum and C. eximium (Figure 2a,d). However, this length pattern was not observed in C. chacoense and C. annuum species, whose predominant lengths were 23 nt followed by 24 nt in C. annuum; 21 nt was the predominant length in C. chacoense at middle fruit stage. Harvested fruit at different developmental stages from days after pollination of (a) C. baccatum, (b) C. annuum, (c) C. chacoense, and (d) C. eximium. Red arrows indicate small fruit (6-dpa) and middle fruit (25-dpa) stages that were collected for small RNA sequencing.

Identification of Known and Novel miRNAs in Pepper
To identify conserved miRNAs in each library, unique sRNA reads were BLASTN searched against currently known mature miRNAs in miRBase (v22.1) allowing one or two mismatches between sequences. Consequently, 22 known miRNAs from 21 conserved families were identified ( Table 1). The only miRNA family that showed variants was miR166, with two members (miR166a and miR166b). To predict novel miRNAs in Capsicum species, sRNA reads that were unmatched to the miRbase were aligned with the genome sequences of C. annuum by using miRDeep2. Ultimately, 27 novel miRNA candidates with stable hairpin structures were identified and designated miR01 to miR27. In addition, precursor sequences from the pepper genome database were obtained (Table S2). Notably, candidate precursors for all the 27 predicted novel miRNAs were suggested. Furthermore, the hairpin structures of these novel small RNAs were predicted by using the RNAfold web server. All the novel miRNAs precursors in pepper possessed typical stem-loop structures and had negative folding free energies ranging from -62.6 to -15.9 according to RNAfold. The predicted secondary structures of six randomly selected novel miRNAs candidates are in Figure 3. In addition, we determined whether these novel miRNA candidates contained miRNA* sequences in the 5p or 3p terminus, which resulted in the identification of the miRNA* sequence in most of the candidate miRNAs, with exception of 6 ( Table S2). The first cleavage position is critical to determine the mature miRNA sequence and resulting target specificity [25,26]. Therefore, base composition of miRNAs plays another important role because it may affect secondary structures and biological properties. Nucleotide bias analysis revealed that 21-to 22-nt miRNAs more frequently contained G and U at the first position; however, 23-and 24-nt miRNAs had a strong preference for G across tissues and species ( Figure S1). Although the nucleotide bias at each position had wide fluctuation, they followed similar patterns through all species and tissues. This was especially observed at the 10th nucleotide, where the predominant base was G followed by C and U, and at the 11th nucleotide, where the four bases were present with dominancy of base C ( Figure S2). Previous studies suggested that the first nucleotide is important for miRNA sorting [26] and that the 10th and 11th nucleotides are responsible for guiding the miRNA to cleave the target mRNA [27]. Altogether, these results suggest that novel miRNAs might be involved in regulating similar physiological and biological processes across Capsicum species.

Expression Pattern of Known and Novel miRNAs in Pepper
We analyzed the sequence frequencies of known and novel miRNAs from the 12 libraries to estimate the pattern expression of the miRNAs found at a specific developmental stage or tissue and to infer their possible roles in the different Capsicum species. Seven of 22 known miRNAs were expressed at different densities in all stages and species (Table 1 and Figure 4): 9, 11 and 12 known miRNAs were expressed at flower, small fruit, and medium fruit, respectively, among all Capsicum species. From the miRNA expression at flower stage, C. baccatum, C. annuum, and C. chacoense were clustered together; however, miR157, 166b, 172 and miR319 showed higher expression in C. annuum. Furthermore, the expression of 13 miRNAs (miR156, 159, 160, 162, 166a, 167, 168, 171, 390, 394, 396, 403, and 6478) at small and middle fruit stages was higher in C. baccatum than the other species. Similarly, miR164, 165 and 398 showed high expression in C. chacoense and C. eximium at middle fruit stage, whereas miR408 was highly expressed in C. eximium at small fruit stage. For novel miRNAs, 8 of 27 were expressed in all tissues and species at different levels (Table 2, Figure 5). Additionally, 15, 11, and 10 were shared between Capsicum species at flower, small fruit, and middle fruit stage, respectively (Figure 5b,c,d). Similar to known miRNAs, some novel miRNAs also showed species-specific expression; for instance, 14 miR-NAs including miR01, 02, 03, 04, 05, 06, 08, 09, 10, 11, 12, 16, 18, and 27 were upregulated in C. baccatum at small and middle fruit stages. Likewise, at the same stages, miR014 showed high expression in C. annuum, whereas miR07 and 20 were highly expressed in C. chacoense. Moreover, miR17, 22, 23, and 24 showed high abundance at flower stage in C. eximium, but miR25 was highly expressed at the same stage in C. chacoense. Finally, miR13, 15, 19, and 21 showed low expression patterns across all species and stages. Base pair probability is represented by color: red, high-probability, and purple, low-probability.

Prediction of Putative Targets for Known and Novel miRNAs in Pepper
To better understand the function of the identified miRNAs and assuming that plant miRNAs have a nearly perfect match to their target mRNAs, putative target genes were predicted by using the psRNA target program. Additionally, we performed 'microRNA:mRNA' seed sequence similarity analysis between known and novel miRNAs and their target mR-NAs because the seed sequence is essential for binding to take place (Mullany et al. 2016) (Table S3). We identified 78 and 60 potential candidate targets for known and novel miR-NAs, respectively (Table S4). The number of putative target genes for a single miRNA ranged from 1 to 6 for known miRNAs and from 1 to 4 for novel miRNAs. The set of predicted target proteins was functionally characterized by using the BLASTKOALA sequence similarity tool. BLASTKOALA annotated 54 (39.1%) of 138 proteins identified. A group of 21 genes was thought to be involved in genetic information processing, but other genes were predicted to participate in processes such as signaling and cellular processes (11), metabolism (7), environmental information processing (3), amino acid metabolism (2), energy metabolism (2), xenobiotics biodegradation and metabolism (2), carbohydrate metabolism (2), lipid metabolism (1), metabolism of terpenoids and polyketides (1), metabolism of cofactors and vitamins (1), and glycan biosynthesis and metabolism (1) ( Figure S3).  Moreover, to gain insights into a global overview of the regulatory functions of the miRNAs across Capsicum species, we analyzed the GO terms for 116 of the 138 identified targets. GO analysis suggested the putative participation of miRNA targets in multiple biological processes, molecular functions, and cellular component ( Figure S4). The major biological processes predicted for these GO-defined target genes were regulation of transcription, DNA-replication (20), oxidation-reduction process (10), protein phosphorylation (10), auxin-activated signaling pathway (5), and hormone-mediated signaling pathway (4). The molecular functions were mostly classified as DNA binding (25), ATP binding (20), protein serine/threonine kinase activity (13), DNA-binding transcription factor (TF) activity (10), and zinc ion binding (8). For cellular components, most genes were related to nucleus (35), membrane (35), plasma membrane (17), cytoplasm (8), and cytosol (8). Furthermore, enriched KEGG pathway analysis of target genes of differentially expressed miRNAs revealed that miRNAs regulate genes involved in amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism, and fructose-mannose metabolism pathways ( Figure S5).

Expression Profile of miRNA Target Genes in C. annuum
To investigate the expression profile of miRNA predicted target genes across different tissues, including leaf, steam, and placenta tissue at 6, 16, and 25 dpa, we used publicly available RNA-seq data for C. annuum cv. CM33428. A heat map was used to visualize tissue-specific expression patterns of target genes ( Figure 6). Overall, 112 of 138 predicted genes were identified, and most genes exhibited unique expression profiles. Genes showing exclusively high expression in leaf tissue encoded proteins for ABC transporter (CA06g14420), cellulose synthase (CA10g10190), FGGY carbohydrate kinase (CA07g00430), or class III HD-Zip protein 3 (CA02g10530) or were involved in the PHO system (CA02g17380). Placenta tissue at 6-dpa showed high expression of genes involved in regulation of transcription such as SPL (squamosa promoter binding protein) domain class TF (CA03g12170), TCP TF (CA08g04030), and NAC (NAM, ATAF, and CUC) domain TF (CA06g18770), thus indicating the participation of these genes in fruit development.
The same tissue at 16-dpa showed high expression of genes associated with plant hormone signaling: auxin response factor (CA04g09370) and serine/threonine-protein kinase (CA04g16720). Likewise, genes exhibiting the same profile encoded proteins such as class III HD-Zip protein (CA12g13110), protein SUPPRESSOR OF GENE SILENCING 3 (CA03g17450), and fructokinase activity (CA02g24240). Meanwhile, genes upregu-lated in placenta tissue at 25 dpa were related to carbohydrate metabolism (CA01g27530, CA03g19360, and CA03g29340), aromatic and neutral amino acid transporter (CA04g05160), leucine rich repeat receptor protein kinase CLAVATA1 (CA02g24570), and TFs with key roles in ripening such as AP2 (APETALA2) and NAC TFs (CA11g14070, CA12g13470). Genes showing high expression during pepper fruit development and ripening at 16 and 25-dpa were the SPL TF (CA02g15200) and the response regulator ARR12-like protein (CA07g02260).

Validation of miRNA and Target Gene Expression
To validate the expression profiles of miRNAs in Capsicum species at flower, small fruit (6-dpa), and middle fruit (25-dpa) stages, we investigated 6 known (miR159, 162, 166a, 319, 396, and 6478) and six novel miRNAs (miR01, 05, 10, 13, 16, and 23) by using stem-loop qRT-PCR based on their sequencing frequencies (Figure 7). In general, expression of most of these miRNAs agreed with sequencing data, with slight variation. The expression profiles of miR162, 166a, 319, 396, 6478, 01, 05, and 16 were consistent with the results from sRNAsequencing. However, the remaining four miRNAs tested showed inconsistent results between stem-loop qRT-PCR and sRNA sequencing in all or specific Capsicum species. For instance, miR159 showed a discrepancy in C. annuum and C. eximium species, miR10 showed differences in only C. eximium, and miR23 showed a different expression trend in C.
annuum. To validate the expression profiles of protein-coding genes targeted by miRNAs, we used qRT-PCR for all Capsicum species at different stages. The expression profiles of auxin-induced protein, class III HD-Zip protein, TCP TF, and serine carboxypeptidase protein, regulated by miR162, 166a, 319, and 016, respectively, showed opposite expression profiles from their target genes, as expected for miRNA targets ( Figure 8). As expected, the expression of the remaining predicted targets was not opposite to their miRNA profiles in at least one species, which suggests that the activity of these genes could be determined by translational repression and/or by multiple miRNAs.

Discussion
Small RNAs, including miRNAs, are key regulators of biological processes such as biotic and abiotic stress tolerance, plant growth and development, metabolic pathways, and morphogenesis [28,29]. In plants, miRNAs regulate gene expression mainly by targeting mRNAs for cleavage and/or by translation inhibition of the target mRNAs during or after transcription [30]. Thus, miRNAs usually negatively regulate the accumulation of mRNAs and show an inverse correlated expression pattern in the same plant cells [31]. Such changes in the expression profile of plant genes play an important role in establishing specific phenotypes between plant species [32]. However, the molecular mechanisms underlying these changes are largely unknown.
Although miRNAs have been identified in pepper fruits by using small RNA-seq [22,23], we lack a systemic comparison study of miRNA expression in pep-per between domesticated and wild species, despite the worldwide importance of these fruits. In this study, we investigated the contribution of gene expression regulation by miRNAs in cultivated (C. baccatum and C. annuum) and wild (C. chacoense and C. eximium) Capsicum species at flower and two fruit development stages (6-and 25-dpa fruits, namely small and medium fruit stages). Our study identified 22 known and 27 novel miRNAs differentially expressed across Capsicum species that may be closely involved in biological processes controlling flower and fruit ripening, including fruit development, morphogenesis, pigmentation, and quality in wild and cultivated pepper species.

Role of miRNAs and Their Regulators in Flowering across Capsicum Species
Flowering is an important biological process for plants that represents the phase transition from vegetative growth to reproductive growth, ensuring success in plant reproduction [33,34]. Several miRNAs have been reported to be a main regulator of the floral phase transition in many plants, including Arabidopsis [35], tomato [36], and apple [37]. For instance, miR156/157 and 172 are the two main key members of the aging pathway regulating SPL and AP2-like TFs [34]. These miRNAs play antagonistic but related roles in Arabidopsis flower induction, whereas a high level of miR156 extends the juvenile phase and delays flowering, and miR172 accumulation leads to early flowering [38,39]. As shown in Figure 4a, miR172 was highly accumulated in flowers in domesticated species, especially C. annuum, so cultivated species may exhibit early flowering as compared with wild species, which could be beneficial for pepper productivity. Similar to miR172, miR319 regulating TCP TFs showed high expression in C. annuum. Although TCP TFs have been more associated with leaf development and petal growth [40], recently Li et al. [41] described the function of these TFs as key regulators of flower development. Likewise, Wang et al. [42] reported that overexpression of TCP8 significantly delayed flowering in Arabidopsis under long-and short-day conditions; additionally, high expression of TCP4 promoted pistil abortion in Prunus mume [43]. Our results suggest that overexpression of miR319 in C. annuum negatively regulates TCP TF expression, which results in early flowering in this cultivated species and needs further confirmation.
Besides miR156/157, 172, and 319, many other miRNA families involved in the control of plant flowering time include miR159, 165/166, 167, 169, 171, 319, 390, and 399 [44]. miR159 regulates GAMYB TFs involved in the gibberellic acid (GA) pathway in different land plants as well as meristem formation and seed development [45,46]. Our results revealed a higher expression of miR159 than other miRNAs across species and stages; however, specifically at flowering stage, all species with the exception of C. chacoense showed similar expression. These results may suggest that GA content is an important player in flower initiation and development across Capsicum species. Furthermore, miR166b and 390 showed high expression in C. annuum and C. baccatum, respectively, and regulated LRR receptor-like serine/threonine-protein kinase (LRR-RLK). Previous studies of Arabidopsis reported the role of these proteins in floral organ abscission [47], in the main abscisic acid-mediated (ABA) signaling pathway and in early ABA perception [48]. Thus, miR166b and 390 may participate as negative regulators of flower organ development through the ABA pathway in cultivated species.
Along with the known miRNAs mentioned above, novel miRNAs were upregulated during flower stage, mainly in wild species. For instance, miR17, 22, 23, and 24 showed high expression in C. eximium, whereas miR25 was highly expressed in C. chacoense. These miRNAs regulate genes coding for disease resistance proteins, suppressor of gene silencing 3-like (SGS3) proteins, and nonsense-mediated mRNA decay proteins, involved in plant pathogen defense and destruction of aberrant mRNAs [49,50]. The novel miRNAs we identified may be linked to plant pathogen protection against virus infection during flowering in wild Capsicum species that represent valuable germplasm resources for crop improvement.

Role of miRNAs and Their Regulators in Pepper Fruit Development
Capsicum species are highly diverse, and fruit attributes are one of the principal phenotypic differences among accessions [51]. Domestication events and continuous selection have increased the variability of fruit features such as shape, size, color, and aroma. Domesticated fruit-bearing crop species have largely increased their fruit size compared with those normally found in progenitor wild species [52,53]. As shown in Figure 1, domesticated species used in this study had large fruits, C. annuum having the largest. However, wild species had smaller fruits, with a round and elongated fruit shape for C. chacoense and C. eximium, respectively. miRNAs orchestrate different fruit development process such as fruit set, formation, shape, size, ripening, and quality in multiple plants including tomato [36,54,55], rice [56], cassava [57], Lycium barbarum [58], sweet potato [59], orange [60], blueberry [61], and diverse cucurbit species [62].
The hierarchical clustering in Figures 4 and 5 for known and novel miRNAs, respectively, revealed that domesticated species were grouped, sharing similar miRNA expression patterns at fruit stages; however, some miRNAs showed high expression in C. baccatum versus other Capsicum species (Figure 5c,d). In addition to their role in plant vegetative versus reproductive phase change, miR156 and 172 have also been associated with fruit size in different species depending on the fruit type [63,64]. For instance, overexpression of miR156a reduced the fruit size and yield in tomato [65], whereas overexpression of miR172 in apple led to reduced fruit size and weight [66]. In our study, miR172 was uniquely expressed in C. baccatum at small fruit stage, whereas miR156 was highly abundant in C. baccatum at both fruit stages. Likewise, miR159 displayed different levels of upregulation in C. baccatum at both fruit stages. Overexpression of miR159 has been associated with abnormal ovule development that affects the initial fruit set, precocious fruit initiation, and seedless fruits in tomato [67]; in contrast, in strawberry, miR159 repression inhibited receptacle ripening and color formation [68]. The increased expression of miR156 and 159 in C. baccatum may interfere somehow in the fruit set and size of this species, which remains to be addressed.
Furthermore, Nimmakayala et al. [69] reported an Ankyrin protein with acyltransferase activity associated with fruit weight in C. annuum. Similarly, Qin et al. [20] reported an acyltransferase involved in the selective sweep signals in the cultivated C. annuum L. (Zunla-1) versus its wild progenitor Chiltepin (C. annuum var. glabriusculum). Our results showed that miR12 with high expression in C. baccatum at middle fruit stage was related to regulation of an acyltransferase protein, which may be involved in the differences in fruit size between cultivated and wild pepper species. As well, we identified a CLAVATA1 gene regulated by miR408. It has been reported in tomato that a mutation in the CLAVATA3 gene involved in the CLAVATA-WUSCHEL pathway promoted stem-proliferation, increasing the meristem size, and resulting in the development of extra organs in flowers and larger fruits [70,71]. These results suggest that CLAVATA genes are implicated in pepper fruit size and are regulated by novel miRNAs. Moreover, Qin et al. [20] also reported an ABC transporter and Pleiotropic drug resistance protein 1 (PDR1) genes likely associated with the morphological and physiological differences between wild and domesticated pepper species. Likewise, it has been reported that the ABC transporter is a candidate gene at the fw3.2 locus associated with fruit weight in tomato [72]. In pepper, ABC transporter has also been associated with fruit weight [69] and capsaicinoid content [73]. In our study, we identified two novel miRNAs, miR20 and miR21, regulating an ABC transporter and PDR4 proteins, suggesting their importance for differences in pepper fruit development between wild and domesticated species.
In addition to their involvement in fruit size, miRNAs are also important for their involvement in the auxin signaling pathway regulating auxin responsive factors (ARFs) [74]. Auxin is a key hormone implicated in most plant organ developmental processes [75] and has been reported to affect fruit shape and development [32]. miR160 regulates ARFs in different plants such as Arabidopsis [76], cotton [77], and tomato [78]; however, in our study, miR162 was also found to regulate an auxin-induced protein. miR160 was downregulated across Capsicum species and fruit stages, whereas miR162 showed high expression at middle fruit stage in C. baccatum. Previous studies have reported that upregulation of miR160 targeting ARF10 affects fruit shape in tomato [79], whereas downregulation alters the abscission of petal and anther, having an effect on tomato fruit set [78]. Together, miR160 and 162 may play an important role in the auxin signaling pathway in Capsicum species and thus are involved in regulating fruit shape and development.
Quality and flavor are two important attributes of pepper fruits. Both features are strongly related to sugar and acid composition [23]. KEGG pathway analysis revealed that some of the miRNAs regulate major candidate genes involved in amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism, and fructose-mannose metabolism pathways. Carbohydrate metabolism-related genes such as cellulose synthase, FGGY carbohydrate kinase, SlArf/Xyl2 and fructokinase proteins were regulated by miR10, 11, 13, and 15, respectively, which suggests that these novel miRNAs may be involved in modulating pepper fruit quality. Particularly, miR13 and 15 were downregulated in all species, whereas miR10 and 11 showed high expression at both fruit stages in C. baccatum. Likewise, novel miR06, which targets a methylenetetrahydrofolate reductase and is involved in folate metabolism, was upregulated in C. baccatum at both fruit stages, so novel miRNAs may also be involved in the concentration of the B vitamin group across Capsicum species. Furthermore, the novel miRNAs miR23 and 24 regulated a gene encoding 1-aminocyclopropane-1-carboxylate oxidase enzyme (ACO), which is directly involved in ethylene biosynthesis [80]. ACO is the rate-limiting step in ethylene production during certain dedicated processes such as fruit ripening in tomato [81,82]. ACO is regulated by the known miR396 in citrus, which inhibits ethylene biosynthesis [83]. Although pepper fruit ripening is widely classified as non-climacteric, patterns of ethylene production and respiratory rates vary across pepper species due to their wide diversity [84,85]. For instance, hot pepper (C. frutescens) fruit ripening is regulated by ethylene and ABA biosynthesis [86]. In addition, some varieties such as 'Camelot', 'King Arthur', and 'Tabasco' exhibit characteristics intermediate between climacteric and non-climacteric fruit ripening [87]. Another miRNA related to fruit ripening in tomato and Lycium barbarum involving the ethylene biosynthesis is miR164, regulating NAC TFs [58]. In tomato, SlNAC4 RNAi-knockout plants showed delayed fruit ripening and a reduction of 30% of total carotenoid content, which suggests that SlNAC4 is a positive regulator of ripening and carotenoid accumulation [88]. In our study, miR164 showed high accumulation at middle fruit stage in C. chacoense, which suggests its role as a negative regulator of fruit ripening in this species. Moreover, fruit ripening and softness are influenced by alpha-N-acetylglucosaminidase and hydrolase (hydrolyzing O-glycosyl compounds) [89,90], which are regulated by miRNAs such as the identified miR168 and 6478.
Narrow genetic diversity which typifies cultivars of domesticated crops has increased the susceptibility of domesticated plants to major diseases [91]. Wild relatives of cultivated crops are recognized as a rich source of genes for disease resistance and stress tolerance. Our results showed different known and novel miRNAs regulating disease resistant genes. For instance, miRNAs such as miR165 and 398 were upregulated in the wild pepper relative C. eximium; these miRNAs were found to target aspartic proteinases CDR1 (constitutive disease resistance 1). Along with their roles in plant development, aspartic proteinases have been associated with immune defense in pepper plants [92]. Additionally, we identified miR394 and miR403 that regulate an F-box protein. F-box genes are well known to be involved in hormone signaling and in response to abiotic stress in pepper [93]. Lim et al. [94] reported a novel F-box protein, CaDIF1 (C. annuum Drought-Induced Fbox Protein 1) in pepper. CaDIF1-silenced pepper plants exhibited a drought-sensitive phenotype, whereas CaDIF1-overexpressing plants exhibited ABA-sensitive and droughttolerant phenotypes. MiR394 and miR403 were upregulated at both fruit stages in C. baccatum, which may suggest a lower accumulation of F-box proteins in this domesticated pepper species and therefore a high susceptibility of C. baccatum to drought stress. In addition, we identified miR403 which regulates an Argonaute (AGO2) protein. AGO proteins bind to small-interfering (si)RNAs and micro (mi)RNAs to target RNA silencing against viruses [95], suggesting that miR403 is directly involved in regulation of pepper defense response to virus infection. Furthermore, miR390 and miR22 associated with a leucine-rich repeat receptor kinase (LRR receptor-like protein) and miR396 associated with a late blight resistance protein were identified. Qin et al. [20] reported 10 disease resistance genes such as LRR and late blight resistance proteins showing a strong selective sweep signals in the cultivated peppers, indicating that these genes seemed to have been affected by selection during domestication. Our results support the involvement of miR390, miR396, and miR22 in pepper domestication. Overall, our results suggest that miRNAs identified, and their target genes, likely serve crucial regulatory roles in pepper related to vegetative to reproductive phase changes, fruit development and quality, and disease resistance across domesticated and wild pepper species. Further research to characterize the microRNAome in Capsicum will be required to confirm these inferences.

Plant Materials
C. baccatum cv. Lemon drop, C. annuum cv. Cayenne, C. chacoense cv. PI 669106, and C. eximium cv. PI 645681 were grown in the greenhouse at West Virginia State University with the following conditions: 16 h of light at 26 • C and 8 h of darkness at 20 • C, relative humidity of 75 ± 2%, and plants were watered daily using manual irrigation. Flowers were uniformly sampled for unopened and fully matured buds from all cultivars, while fruit samples were collected at 6 and 25 days post-anthesis (dpa) from all four cultivars. All samples were harvested in three biological replications and frozen immediately in liquid nitrogen and stored at −80 • C for RNA isolation.

Construction and Sequencing of Small RNA Libraries
Total RNA was isolated by using TRIzol reagent (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions. RNA purity was achieved by using the RNA clean and concentrator kit (Zymo Research, Irvine, CA, USA) along with on-column DNaseI digestion to remove genomic DNA. Total RNA quality and quantity were measured by using the Agilent 2100 Bioanalyzer and Qubit 4 Fluorometer (Invitrogen, Carlsbad, CA, USA), respectively. Total RNA from three biological replicates was pooled for each sample before small RNA-Seq library preparation. At least 1.5 µg high-quality total RNA was used to construct each small RNA library and for sequencing by the Beijing Genomics Institute (Hong Kong) with the Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA). The small RNA-seq dataset was deposited in the NCBI SRA database under the accession number PRJNA718887.

Bioinformatics Analysis of miRNAs
After sequencing, the quality of raw reads was ascertained by checking the adapter, GC distribution, average base content, and quality score of the distribution by using the fastqc program. The Cutadapt toolkit was used for removal of poor-quality reads, sequences with a poly-A tail, low-quality reads with ambiguous bases ("N"), reads with <18 or >30 bases, and adapter sequences [96]. Furthermore, the clean reads were BLAST-searched against the Rfam database (v14.1) [97] and Repbase [98] to identify and exclude known noncoding RNAs including rRNAs, scRNA, snoRNAs, snRNAs, and tRNAs. The remaining clean reads were searched against the miRBase database (v22.1, http://www.mirbase.org/, accessed on 1 August 2020) to identify known putative miRNAs [17]. The final miRNAs dataset underwent sequence length distribution and nucleotide preference analysis at each position. Meanwhile, the remaining unannotated sRNA sequences that were not mapped to any pre-miRNAs in miRbase were analyzed by using miRDeep2 [99] to predict potential novel miRNAs using the C. annuum cv. CM334 reference genome [19]. The prediction of novel miRNAs was based on previously reported criteria including (a) no more than four mismatches between the small RNA and the target (G-U bases count as 0.5 mismatches), (b) no more than two adjacent mismatches in the miRNA/target duplex, (c) no adjacent mismatches in positions 2 to 12 of the miRNA/target duplex (5 of miRNAs), (d) no mismatches in positions 10 to 11 of the miRNA/target duplex, (e) no more than 2.5 mismatches in positions 1 to 12 of the miRNA/target duplex (5 of miRNAs), and (f) a minimum free energy (MFE) of the miRNA/target > 75% [100,101]. Additionally, the putative precursor sequence was folded for each candidate novel miRNA by using RNAfold from the Vienna RNA software package [102]. Similarly, as for known miRNAs, for all the predicted novel miRNAs, properties including miRNA count, length, and nucleotide bias at each position were determined.

Analysis of Differentially Expressed miRNAs
We used the R package DEGseq [103] to identify differentially expressed known and novel miRNAs across species and stages. The fold changes of miRNA expression were normalized by transcripts per million according to the formula: Normalized expression (NE) = Actual miRNA reads count/Total count of clean reads × 1,000,000 (Zhou et al., 2010).

Prediction of miRNA Targets and Enrichment Analyses
The putative target genes of identified miRNAs were predicted by using psRNATarget [104] with a pepper transcriptome dataset [20]. Genes with a final sequence score ≤ 5 were considered potential candidate miRNA targets. The target genes were functionally annotated by using the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/, accessed on 1 September 2020) pathway database with the BlastKOALA sequence similarity tool [105]. Pathways with false discovery rate (FDR) ≤ 0.5 were considered significantly enriched. Gene Ontology (GO) analysis of these genes involved using Blast2GO with a cutoff E-value of 10 −5 (http://www.blast2go.com, accessed on 4 May 2021) [106]. The Reads Per Kilobase of transcript expression values from leaf, stem, and placenta tissues (6, 16, 25 dpa) for the identified target genes were retrieved from published RNA-seq data [19] and used to generate a heatmap with the ClustVis web tool (https://biit.cs.ut.ee/clustvis/, accessed on 1 September 2020).

Quantitative RT-qPCR and Stem-Loop RT-qPCR
To validate our results from the bioinformatics-based analysis, we used stem-loop RT-qPCR for miRNAs and RT-qPCR for target genes in flower, small fruit (6-dpa) and medium fruit (25-dpa) tissues with three biological replications across the Capsicum species analyzed. The RT-qPCR primers were designed from miRNA sequences as described [107] (Table S5). Total RNA used for small RNA sequencing was used for reverse transcription using Superscript IV reverse transcriptase. RT-qPCR analysis was achieved with the Universal SYBR Green Master Mix (Thermo Fisher Scientific, USA) on the StepOne Plus Real-Time PCR System (Applied Biosystems, USA) with the following conditions: 95 • C for 10 min, followed by 40 cycles of 95 • C for 15 s, 56 • C for 30s, and 72 • C for 15 s, and melting curve analysis from 65 to 95 • C. To normalize gene expression, small nuclear RNA (snRNA) U6 and β-tubulin were internal controls for miRNA and target genes, respectively. All RT-qPCR reactions involved three technical replicates, and the relative gene expression of miRNAs and targets was estimated by the 2 -∆∆Ct method [108]. The oligonucleotide primers corresponding to the predicted target genes are in Table S6.

Conclusions
Considering the global importance for pepper worldwide, we used high-throughput sequencing and identified 22 known miRNAs and 27 novel miRNAs differentially expressed in flower and pepper fruits at different developmental stages across 4 Capsicum species. Analysis of differential expression patterns combined with target prediction suggested key roles for these miRNAs in controlling flower time and pepper fruit development in cultivated and wild species. The results expand the study of miRNAs in plants by providing a better understanding of their essential roles in miRNA-based regulation processes in pepper. Similarly, our results provide insight into the biology, evolution, and domestication process of Capsicum species, accelerating the agricultural applications of the miRNAs, the genes of their biogenesis pathway and providing targets for future investigation in pepper and other plants.