Conservation and Diversity of miR166 Family Members From Highbush Blueberry (Vaccinium corymbosum) and Their Potential Functions in Abiotic Stress

MicroRNA166 (miR166) is highly conserved and has diverse functions across plant species. The highbush blueberry (Vaccinium corymbosum) genome is thought to harbor 10 miRNA166 loci (Vco-miR166), but the extent of their evolutionary conservation or functional diversification remains unknown. In this study, we identified six additional Vco-miR166 loci based on conserved features of the miR166 family. Phylogenetic analyses showed that mature Vco-miR166s and their precursor cluster in several clades are evolutionary conserved with diverse species. The cis-regulatory elements in the Vco-miR166 promoters indicated functions related to different phytohormones and defense responses. We also identified putative targets of vco-miR166s, which targeted the same gene families, suggesting the functional conservation and diversification of Vco-miR166 family members. Furthermore, we examined the accumulation patterns of six mature Vco-miR166s in response to abiotic stresses by stem-loop reverse RT-qPCR, which revealed their upregulation under freezing, cold, and heat stress, while they were downregulated by drought compared to control growth conditions. However, Vco-miR166 members showed different expression patterns when exposed to salt stress. These results showed that conserved Vco-miR166 family members display functional diversification but also coordinately influence plant responses to abiotic stress.

Introduction

MicroRNAs (miRNAs) are small [20–24 nucleotides (nt)] non-coding regulatory RNAs derived from longer precursor (pre-miRNA) molecules with stem-loop structures. The pre-miRNA is processed by DICER LIKE1 (DCL1) into the mature miRNA, which then interacts with its target transcripts. Mature miRNAs mediate almost all plant cellular and metabolic processes by regulating gene expression at the posttranscriptional level (Thakur et al., 2011; Fang and Wang, 2021).

In plants, miR166 is a highly conserved miRNA implicated in a wide range of cellular and physiological processes by regulating its cognate target genes encoding members from the homeodomain leucine zipper (HD-ZIP III) family of transcription factors (Miyashima et al., 2013; Singh et al., 2014; Yadav et al., 2021). Recent evidence showed that the miR166 family plays crucial roles in response to abiotic stress. In mulberry (Morus multicaulis), miR166f targets the transcripts encoding an HD-ZIP III transcription factor whose expression is induced by drought stress and is a positive regulator of drought stress tolerance (Li R. et al., 2018). In Arabidopsis (Arabidopsis thaliana), multiple miR166 members accumulate when seedlings are exposed to low temperature (Zhou et al., 2008). In potato (Solanum tuberosum), the levels of miR166 family members were downregulated by 4°C treatment but were upregulated by salt stress (Kitazumi et al., 2015; Ou et al., 2015). In barley (Hordeum vulgare), miR166 was induced in leaves upon dehydration stress (Kantar et al., 2010). In wild emmer wheat (Triticum turgidum), miR166 was downregulated by drought stress (Kantar et al., 2011). In soybean (Glycine max), miR166 members play an important role in response to cold, drought, and salinity stress (Li et al., 2017).

Blueberries (Vaccinium spp.) are an economically important small fruit crop. The main distribution of commercial blueberries from northern to southern China comprises lowbush blueberry (V. angustifolium), highbush blueberry (V. corymbosum), and rabbiteye blueberry (V. ashei). The complement of miRNAs has been identified in blueberries using high-throughput sequencing. For example, 412 conserved miRNAs belonging to 29 families were identified in fruits from rabbiteye blueberry (Yue et al., 2017); similarly, 127 known and 101 novel miRNAs were identified from flowers, young fruits, and ripe fruits of rabbiteye blueberry (Li G. et al., 2018). In highbush blueberry, 84 known miRNAs and 16 novel miRNAs were identified from fruits, with 10 miRNAs from miR166 family (Hou et al., 2017). However, the function of miR166s in highbush blueberry fruits and the genomic sequence of miR166 loci are unknown. In particular, the response of the highbush blueberry miR166s to abiotic stress has not been reported.

As mature miRNAs are highly conserved within the same family, sequence prediction allows the identification of multiple members of a given miRNA family, including miR166 (Li et al., 2017). The sequence and assembly of the highbush blueberry genome were recently completed, opening the doors to a systematic identification of miRNAs and their precursors in this species (Colle et al., 2019). In this study, we predicted the sequences of miRNA166 family members, their precursors, and their encoding genes by mining the highbush blueberry genome sequence. We demonstrated the evolutionary conservation and functional diversification of this family based on a sequence analysis of mature Vco-miR166s, Vco-miR166 loci, and Vco-miR166 promoters and predicted vco-miR166 target genes. We also examined the accumulation patterns of six Vco-miR166s in response to abiotic stresses. The results help elucidate the functions of the miR166 family in highbush blueberry.

Materials and Methods

Plant Materials and Stress Treatments

One-year-old plants from the highbush blueberry cultivar ‘Bluecrop’, grown in a growth chamber, were used for salt, drought, freezing, cold, and heat treatments. Plants were also maintained in tissue culture on modified woody plant medium with Murashige and Skoog (MS) vitamins, 1.0 mg/L trans-zeatin, and 3% (w/v) sucrose. Plants were subcultured every 28 d and were used for drought stress. All plants used for stress treatments were grown at 25°C under a 16-h-light/8-h-dark photoperiod. Plants were exposed to −5°C for 0, 0.5, 1, 2, 4, or 6 h for freezing treatment; 4°C for 0, 2, 4, 6, 9, 12, for 24 h as cold treatment; and 40°C for 0, 0.5, 1, 1.5, 2, and 3 h as heat treatment. Plants were watered with a 200 mM NaCl solution for 0, 1, 2, 6, 9, 12, or 24 h as salt treatment. For drought treatment, plants in tissue culture were transferred to medium (overlaid with 20% (w/v) PEG-6000 for 1 day, and then PEG-6000 was removed) for 0, 1, 2, 6, 9, 12, and 24 h. For all treatments, the 0-h treatment served as the control. Leaves were collected and immediately frozen in liquid nitrogen for microRNA extraction.

Identification of Novel MiR166s in Highbush Blueberry

The precursor sequences of miR166 from sweet orange (Citrus sinensis), grape (Vitis vinifera), and apple (Malus domestica) were downloaded from the miRbase database (https://www.mirbase.org) and used as a query against the V. corymbosum cv. Draper v1.0 genome scaffolds database (https://www.vaccinium.org/). Homologs were screened as previously described (Li et al., 2017; Hou et al., 2017). The secondary structures and the negative minimal folding free energy (MFE) of each miRNA were assessed and computed using the RNAfold WebServer (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) (Mathews et al., 2004; Gruber et al., 2008). The minimal folding free energy index (MFEI) was calculated with the following equation: MFEI = [(MFE/length of the RNA sequence) *100]/(G + C)%. The candidates were then confirmed by comparison with miR166s from other species using the NCBI database (https://www.ncbi.nlm.nih.gov). The mature miR166s of highbush blueberry were predicted by searching against the miRbase database using the precursor sequences of miR166 of highbush blueberry.

Phylogenetic Analysis and Sequences Alignment

Phylogenetic trees were separately constructed using the UPGMA method with the MEGA X program from miR166 precursors and mature miR166s from highbush blueberry (Vco-miR166s), grape (Vvi-miR166s), sweet orange (Csi-miR166s), and apple (Mdm-miR166s) (Kumar et al., 2018). Sequence alignments were performed using DNAMAN software, version 8.192.

Promoter Analysis of Vco-miR166s

The upstream sequences of Vco-miR166s were retrieved from the highbush blueberry genome (https://www.vaccinium.org/). The transcription start sites (TSSs) and TATA-box of Vco-miR166s were predicted by TSSP software in Softberry (http://linux1.softberry.com/berry.phtml?topic=tssp&group=programs&subgroup=promoter). About 1,600 bp of sequence upstream of the TSS was analyzed for the presence of known cis-acting regulatory sequences held in the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Prediction and Functional Analysis of Vco-miR166 Targets

The target sequences of highbush blueberry Vco-miR166s were predicted with the psRNATarget program (http://plantgrn.noble.org/psRNATarget/) against American cranberry (Vaccinium macrocarpon) transcripts (Polashock et al., 2014; Dai et al., 2018). The obtained target sequences were then used to query target genes in the V. corymbosum cv. Draper v1.0 genome transcript database (Colle et al., 2019). The functions of Vco-miR166s target genes were predicted by BLAST against the Swiss-Prot database (https://ww.uniprot.org). Genetic similarity and genetic distance were calculated by DICE and Nei’s similarity coefficient, respectively. The phylogenetic tree was constructed using the UPGMA method with NTSYSpc-2.11 F software.

Reverse-Transcription Quantitative PCR (RT-qPCR)

The microRNAs were extracted using an EASYspin Plant microRNA rapid extraction kit (Aidlab, Beijing, China) from control and abiotic stress samples. Reverse transcription and qPCR were performed with the Mir-X miRNA qRT-PCR TB Green Kit (Takara Inc., Dalian, China) using specific stem-loop RT and forward RT-qPCR primers. qPCR was conducted using the ABI StepOnePlus Real-Time PCR System using U6 as the reference. Primer sequences are listed in Supplementary Table S1.

Results

Identification and Characterization of Vco-miR166s in Highbush Blueberry

Based on the conserved features shared by pre-miR166s sequences from sweet orange (Citrus sinensis), grape (Vitis vinifera), and apple (Malus domestica), we performed a BLAST search of the highbush blueberry genome. We identified three new pre-miR166s, designated here Vco-miR166g, h, and i, following a previous report (Hou et al., 2017) (Table 1). The length of these three new pre-miRNAs ranged from 107 to 186 nt, with three to four predicted mismatches between the 5′ and 3′ arms of the miRNA (Supplementary Figure S1). The minimal folding free energy (MFE) of these new mature Vco-miRNAs was lower (from −67.9 to −43.30) and their MFE index (MFEI) (from 0.85 to 1.08) was higher than those of transfer RNAs (tRNAs), ribosomal RNAs (rRNA), and messengers RNAs (mRNAs). These results suggest that these predicted precursors may be bona fide precursors of highbush blueberry miRNAs.

TABLE 1

TABLE 1. Characteristics of three newly identified miR166 precursors in highbush blueberry.

BLAST analysis determined that Vco-miR166g, Vco-miR166h, and Vco-miR166i precursors share 90%, 90%, and 63% identity with Csi-miR166g, Vvi-MIR166g, and Sly-miR166 b precursors, respectively. We thus concluded that these precursors likely belong to the miR166 family of miRNAs (Table 1). At the same time, these three precursors from highbush blueberry showed conserved features in the 3′ arm and differ by several bases in the 5′ arm with other plant species (Supplementary Figure S2).

From the three pre-miR166s above, we predicted six putative mature miR166 sequences (Table 2), located along the 3′ or 5′ arm of the secondary stem-loop structure of each pre-miR166 sequence by searching against the miRbase database. We designated these candidate miRNAs as Vco-miR166g-3p, Vco-miR166g-5p, Vco-miR166h-3p, Vco-miR166h-5p, Vco-miR166i-3p, and Vco-miR166i-5p. We noticed that the sequence of Vco-miR166i-3p was identical to that of Vco-miR166b-3p, while Vco-miR166i-5p was identical to Vco-miR166d-5p. In addition, BLAST analysis showed that Vco-miR166g-3p, Vco-miR166h-3p, and Vco-miR166i-3p are identical in sequence to known mature miR166s from other species, with the number of mismatches ranging between one (for Vco-miR166g-5p and Vco-miR166h-5p) and two (for Vco-miR166i-5p). These results thus indicated that these new Vco-miR166 candidates do in fact belong to the miR166 family and that the 3′ arm sequences are more conserved than the 5′ arm.

TABLE 2

TABLE 2. Characteristics of six newly identified mature miR166s in highbush blueberry.

Conservation and Diversification of the Vco-miR166 Family

We also aligned the sequences of pre-miR166s from highbush blueberry and other major fruit trees (sweet orange, grape, and apple) and produced the corresponding phylogenetic tree (Figures 1, 2). We observed the highest degree of conservation along the arms of the secondary stem-loop structure, especially on the 3′ arm and in the core region of the mature miRNA (GGA​CCA​GGC​TTC​aTT​CC), with more variation flanking the mature miRNA (Figure 1). The nine Vco-miR166s roughly clustered into clades Ⅰ and Ⅱ but largely did not cluster together within each clade (Figure 2A). These results supported the evolutionary conservation of miR166 family members across various plant species.

FIGURE 1

FIGURE 1. Sequence alignment of miR166 precursors. Vco, Vaccinium corymbosum; Csi, Citrus sinensis; Vvi, Vitis vinifera; Mdm, Malus domestic. Pink, 100% conservation; green, ≥75% conservation; yellow, ≥50% conservation. The core region of the mature miRNA is indicated in the black box.

FIGURE 2

FIGURE 2. Phylogenetic analyses of miR166 loci and their predicted 3′ and 5′ arms. Sequences were aligned in MEGA X, and the phylogenetic tree was constructed with the UPGMA method (A) Phylogenetic analysis of Vco-miR166 precursors and related sequences from Csi, Citrus sinensis; Vvi, Vitis vinifera ; and Mdm, Malus domestica. Csi-MIR167 b served as the outgroup (B,C) Phylogenetic analysis of mature Vco-miR166s and Csi-miR166s for the 3′ arm (B) and 5′ arm (C). Black and red dots indicate mature miR166s or their precursors of highbush blueberry from a previous study (Hou et al., 2017) and this study, respectively. The numbers at the nodes indicate the percentage of support from 1,000 bootstrap replicates. The scale indicates the average number of amino acid substitutions per site.

We repeated the phylogenetic analysis separately on the 3′ and 5′ arms of nine mature highbush blueberry miR166 sequences along with related sequences from sweet orange (Figures 2B,C). For the analysis of the 3′ arm, Vco-miR166i-3p and Vco-miR166b-3p were identical to Csi-miR166a, c, d, and e-3p, while Vco-miR166g-3p was identical to Csi-miR166b, g, and f-3p in clade I. Eight Vco-miR166s and seven Csi-miR166s clustered into clade I. For the 5′ arm of miR166 sequences, eight Vco-miR166s and eight Csi-miR166s formed four clades, with only clade IV not including a Vco-miR166 member. These results supported the notion that the miR166 family is highly conserved in diverse species. Notably, Vco-miR166e-3p appeared to be distant from other Vco-miR166s based on its position in the phylogenetic tree, indicating the evolutionary diversification of Vco-miR166 members.

Putative cis-Regulatory Elements of Vco-miR166 Promoters

To explore the potential functions of Vco-miR166s, we analyzed the putative cis-regulatory elements contained within 1,600 bp of the Vco-miR166 promoter sequence upstream of the TSS (+1 bp) (Supplementary Table S2). We identified core promoter elements (TATA-box), common cis-acting elements, enhancer regions (CAAT-box), as well as light-, phytohormone-, and defense-related elements in all Vco-miR166 promoter regions. Several Vco-miR166 promoters also contained elements related to the cell cycle and meristem expression, seed development, synthesis, and metabolism (Supplementary Table S3). Table 3 summarizes the phytohormone- and defense-related elements in Vco-miR166 promoters. For plant hormone–related elements, 77.8% (7/9) of Vco-miR166 promoters harbored gibberellin- and abscisic acid–responsive elements, with 66.7% (6/9) of Vco-miR166 promoters also presenting elements involved in methyl jasmonate and auxin responsiveness. By contrast, only Vco-miR166a and Vco-miR166h promoters contained elements associated with salicylic acid responsiveness. For defense-related elements, we detected anaerobic induction elements in all Vco-miR166 promoter regions, while drought inducibility, low-temperature, and defense- and stress-responsive elements were present in four, three, and two of the nine Vco-miR166 promoters, respectively. These results indicated that Vco-miR166 loci may be involved in responses to different phytohormones and defense, but also reflected their potential functional diversification.

TABLE 3

TABLE 3. Putative cis-regulatory elements in upstream sequences of blueberry Vco-miR166s.

Prediction of Vco-miR166 Target Genes

To obtain a better understanding of Vco-miR166 functions, we predicted their targets with the psRNATarget program against the transcript database for American cranberry (V. macrocarpon) transcripts, another Vaccinium species, because there is no highbush blueberry database in the psRNATarget program (Supplementary Table S4). We then used the American cranberry targets to search for related sequences in the highbush blueberry (V. corymbosum cv. Draper v1.0) transcripts (Supplementary Table S5). All Vco-miR166s were predicted to target the transcripts of a gene encoding a pentatricopeptide repeat–containing protein (Supplementary Table S6). In addition, all Vco-miR166s except Vco-miR166h-3p were predicted to target the transcripts encoding a leucine-rich repeat (LRR) receptor-like serine/threonine-protein kinase. Many Vco-miR166s-3p were predicted to regulate genes encoding homeobox-leucine zipper proteins (HD-ZIP III), while many Vco-miR166s-5p were predicted to target the transcripts for genes encoding phosphatidylinositol/phosphatidylcholine transfer proteins and patellin. These results supported the functional conservation of Vco-miR166 family members. Several Vco-miR166s involved in abiotic stresses, including heat, stress response, and cold, were predicted to target transcripts for cold-regulated proteins, while other Vco-miR166s appeared to target the transcripts encoding many protein families. For example, Vco-miR166a,b,i-3p targeted ethylene-responsive and basic helix-loop-helix (bHLH) transcription factors, Vco-miR166c-3p targeted WD repeat–containing protein, Vco-miR166d-3p targeted bHLH transcription factors, and Vco-miR166e,g-3p targeted WRKY transcription factors. These results illustrated the functional diversification of Vco-miR166 family members in highbush blueberry.

To explore the diversity of Vco-miR166 functions, we calculated the genetic similarity and genetic distance according to the presence or absence of functions of their target genes in Supplementary Table S6. The genetic similarity was the highest (0.88) and the genetic distance was the closest (0.12) between Vco-miR166a-3p and Vco-miR166b/i-3p, whereas Vco-miR166b/i-3p and Vco-miR166b-5p exhibited the lowest genetic similarity (0.05) and the farthest genetic distance (3.03) (Table 4). The genetic similarity between Vco-miR166h-3p and Vco-miR166i-3p was also low at only 0.12. In fact, a phylogenetic tree according to the genetic similarity between miR166 targets showed that Vco-miR166s are divided into two clades, with Vco-miR166s derived from the same arm clustering together (Supplementary Figure S3). These results indicated that the functions of Vco-miR166s are diverse and that Vco-miR166s from the same arm exert similar regulatory functions.

TABLE 4

TABLE 4. Genetic distance and genetic similarity between Vco-miR166s from highbush blueberry according to the functions of target genes.

Accumulation of Vco-miR166s Under Abiotic Stress

To explore how Vco-miR166s are related to abiotic stress, we selected six Vco-miR166s to examine their abundance in response to abiotic stresses (Figure 3). The six Vco-miR166s had similar accumulation patterns under freezing stress, as they all rapidly increased in abundance, with Vco-miR166h-3p and Vco-miR166i-3p increasing in abundance by 18.6- and 28.6-fold upon 2 h of freezing stress compared to the control at 0 h. The abundance of the six Vco-miR166 transcripts sharply declined and showed no accumulation after 6 h of cold exposure in 1-year-old blueberry plants, which died after 6 h of –5°C treatment (Figure 3A). These results suggested that these six Vco-miR166s, especially Vco-miR166h-3p and Vco-miR166i-3p, may be involved in freezing response and as negative regulators.

FIGURE 3

FIGURE 3. RT-qPCR analysis of Vco-miR166 abundance under various stress conditions. Highbush blueberry plants were exposed to freezing (A), cold (B), heat (C), drought (D), and salt (E) stress. Data represent the means ± standard deviation (SD) of three independent biological replicates, each with three technical replicates. Different letters indicate significant differences (p < 0.05) between samples, as determined by Tukey’s test.

In blueberry plants exposed to cold stress, the abundance of Vco-miR166a, Vco-miR166e-3p, and Vco-miR166i-3p reached their maximum levels after 9 h of freezing stress, corresponding to an increase of 4.1-, 10.2-, and 3.9-fold compared to the control at 0 h, respectively, followed by a sharp decrease to control levels for Vco-miR166a and Vco-miR166e-3p or to lower levels than control for Vcc-miR166i-3p. The abundance of Vco-miR166d-5p, Vco-miR166f-5p, and Vco-miR166h-3p was upregulated and reached its highest levels at 24, 4, and 12 h during freezing stress, increasing 12.8-, 15.4-, and 8.5-fold compared to the control samples, respectively (Figure 3B). Thus, the Vco-miR166 family appeared to be involved in cold stress responses.

When blueberry plants were exposed to heat stress, Vco-miR166a precursor abundance was upregulated and increased by 38.1-fold at 2 h compared to the control samples at 0 h, followed by a sharp decrease at 3 h. Vco-miR166d-5p and Vco-miR166f-5p showed similar accumulation patterns, with their levels reaching their peak at 2 h and their lowest levels at 1.5 and 3 h. The accumulation of Vco-miR166e-3p, Vco-miR166h-3p, and Vco-miR166i-3p was also upregulated and reached its highest levels at 0.5 h for Vco-miR166h-3p, at 1 and 1.5 h for Vco-miR166e-3p, and at 0.5 and 2 h for Vco-miR166i-3p. All six miRNAs showed a rapid decrease in abundance at 3 h (Figure 3C). Therefore, Vco-miR166 family members might respond during early stages of heat stress.

Under drought stress, six Vco-miR166 precursors showed similar trends in their abundance, with a sharp decrease and reaching their lowest levels at 12 and 72 h and a modest increase in their levels at 24 and 48 h (Figure 3D). Vco-miR166 members might thus promote the accumulation of their target transcripts and play a positive role in stress adaptation.

In the context of salt stress, Vco-miR166a transcript levels were upregulated after 1 h and then rapidly decreased, with the lowest abundance after 12 h. We observed similar patterns for Vco-miR166d-5p, whose transcript levels gradually increased by 6.9- and 6.0-fold at 6 and 9 h compared to the 0 h control and then rapidly increased until 24 h. Vco-miR166h-3p also increased in abundance before decreasing, but increased rapidly at 24 h. On the contrary, the transcript levels of Vco-miR166e-3p, Vco-miR166f-5p, and Vco-miR166i-3p were downregulated compared to the control under salt stress, especially Vco-miR166e-3p, whose transcript levels rapidly decreased at 1 h and then remained low after salt stress (Figure 3E). These results suggested that these Vco-miR166 members might respond to salt stress at an early stage.

Discussion

miRNAs play various roles in regulating almost all processes of plant growth and in metabolism. Several conserved miRNAs were reported in highbush blueberry based on an analysis of expressed sequence tags (ESTs) (Li et al., 2014), leading to the identification of 10 miR166 family members using high-throughput sequencing of small RNAs and reference sequences from the American cranberry genome, and ESTs and transcriptome data from blueberries during fruit maturation (Hou et al., 2017). The completion of the highbush blueberry genome (Gupta et al., 2015) allowed us to identify three new Vco-pre-miR166 members and determine their mature miR166s based on the V. corymbosum cv. Draper v1.0 genome. We determined here that the newly identified pre-miR166s exhibit low MFE and high MFEI values and form stable secondary hairpin structures (Table 1 and Supplementary Figure S1). In addition, we showed that the mature Vco-miR166 sequences are highly identical to other miR166 family members from other species (Table 2). We thus concluded that the three newly identified Vco-miR166 loci belong to the miR166 family.

The conserved nature of miR166 family members has been demonstrated in many plant species (Li et al., 2017; Barik et al., 2014). In this study, we showed that the newly identified mature Vco-miR166s derived from the 3′ and 5′ arms have an identical sequence and differ by one or two bases with mature miR166s from other plant species (Table 2). Phylogenetic analysis indicated that some Vco-miR166s are identical to mature Csi-miR166s, while the precursor sequences of miR166s from blueberry, sweet orange, grape, and apple also showed a high degree of conservation along the 3′ arm (Figure 1). Thus, miR166 family members are highly conserved across diverse plant species, especially along the 3′ arm (Barik et al., 2014). Mature miR166 members are known to regulate the abundance of transcripts encoding HD-ZIP III transcription factors in many species (Boualem et al., 2008; Li R. et al., 2018; Yadav et al., 2021). In this study, we predicted that all Vco-miR166s derived from the 3′ arm target HD-ZIP III transcripts, supporting the functional conservation of miR166 family members in various species (Supplementary Table S6).

The multiple miR166 family members likely diverged from a common ancestor, as reflected by the evolutionary diversification of the miR166 precursor sequences (Barik et al., 2014; Li et al., 2017). In addition, an analysis of the Vco-miR166 promoters and the putative Vco-miR166 target genes indicated that Vco-miR166 exert diverse functions (Li et al., 2017). Indeed, we detected different defense-related elements and anaerobic induction elements in the promoters of Vco-miR166 loci (Table 3). The prediction of Vco-miR166 target genes showed that most Vco-miR166s target different genes (Table 4). These results suggest the functional diversification of Vco-miR166s, in agreement with similar results obtained in other plant species (Barik et al., 2014; Li et al., 2017).

Temperature is an important environmental factor that affects plant growth, yield, and geographical distribution. Plants are constantly subjected to freezing, cold, and heat stress from their surrounding natural environment, with miRNAs playing an active role in temperature stress responses (Sun et al., 2019; Fang and Wang, 2021). In rice (Oryza sativa), Osa-miR166 m and Osa-miR166 k are involved in cold stress (Devi et al., 2013). Similarly, Ath-miR166 abundance is upregulated in Arabidopsis seedlings exposed to cold stress (Zhou et al., 2008), while Stu-miR166 (from tomato) and miR166a-e [from almond (Prunus dulcis)] abundance decreased under cold stress (Ou et al., 2015; Karimi et al., 2016). In tomato, upregulation of Sly-miR166 and the concomitant downregulation of its target gene encoding an HD-ZIP III showed that Sly-miR166 plays a critical role in regulating cold stress responses in this species (Valiollahi et al., 2014). Under chilling conditions, miR166e-3p abundance in field mustard (Brassica rapa) was upregulated, while miR166a levels were downregulated; in almond, the abundance of miR166a-e decreased upon exposure to cold (Zeng et al., 2018; Karimi et al., 2016). In the context of heat exposure, miR166 abundance decreased in radish (Raphanus sativus) but increased in rice (Wang et al., 2015; Li et al., 2014a). However, miR166a/b/c-3p abundance increased while that of miR166c/f-5p decreased in Betula luminifera under heat stress (Pan et al., 2017). The abundance of almost all Vco-miR166s characterized here transiently increased before dropping under temperature stress. In particular, the levels of Vco-miR166s rapidly increased within 1 or 2 h under freezing stress (Figure 3A). At the same time, the promoters of several Vco-miR166s (Vco-miR166a, Vco-miR166g, and Vco-miR166h) harbored a low-temperature response element; notably, Vco-miR166s also targeted transcripts encoding HD-Zip III transcription factors, which are involved in temperature stress (Valiollahi et al., 2014). The functions of predicted target genes showed that Vco-miR166f-5p and other Vco-miR166s may target heat and cold shock proteins and directly regulate temperature stress responses (Supplementary Table S6). Our results indicated that Vco-miR166s may be early positive regulators of the response to temperature stress.

Drought is another major abiotic stress that limits plant growth and development. Many miRNAs regulate drought responses in plants (Li et al., 2013; Budak et al., 2015; Pagliarani et al., 2017). Indeed, miR166 expression was shown to be downregulated in barley and in emmer wheat, while being upregulated in rubber tree (Hevea brasiliensis) in response to drought (Kantar et al., 2010; Kuruvilla et al., 2019; Kantar et al., 2011). In rice, manipulating the levels of miR166 demonstrated its contribution to drought resistance (Zhang et al., 2018). In white mulberry (Morus alba multicaulis), miR166f might function as a positive regulator of drought stress (Li R. et al., 2018). In this study, we showed that the abundance of all Vco-miR166s decreases under drought stress. Thus, the Vco-miR166 family might act as a negative regulator of drought stress.

Salt stress is a major abiotic stress that affects plant growth and development. Many miRNAs, including miR166 family members, play important roles in response to salt stress in plants (Si et al., 2014; Sun et al., 2019; Fang and Wang, 2021). The abundance of miR166a-5p and miR166j-3p was upregulated in soybean roots, while that of miR166a-5p was downregulated in barley roots (Wang et al., 2020; Kuang et al., 2021). However, miR166 family members were upregulated in roots and downregulated in leaves of sweet potato (Ipomoea batatas) (Yang et al., 2020), while Pgu-miR166t was also downregulated in the leaves of common guava (Psidium guajava) (Sharma et al., 2020). In our study, the levels of Vco-miR166e-3p, Vco-miR166f-5p, and Vco-miR166i-3p were downregulated in highbush blueberry leaves, while those of Vco-miR166d-5p and Vco-miR166h-3p were upregulated upon exposure to salt stress relative to controls. However, Vco-miR166a levels rose at 1 h and dropped at 12 h following exposure to salt stress (Figure 3E). These results indicated that different members of the Vco-miR166 family exhibit distinct expression patterns and form a complex regulatory network under salt stress.

Here, we identified three new miR166 loci in the highbush blueberry genome, which increased the size of the highbush blueberry miR166 family. The resulting 16 Vco-miR166 members of highbush blueberry exhibited evolutionary conservation and functional diversification based on sequence analysis, promoter dissection, and target gene predictions. The analysis of their upstream sequences and stem-loop RT-qPCR showed that Vco-miR166 family members are involved in various abiotic stresses, possibly by positively and/or negatively regulating target genes.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

CZ conceived and designed the project. YL, XW, QG, XZ, LZ, and YZ participated in the experiments and data analysis. YL and XW drafted the manuscript. CZ modified the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by the Natural Science Foundation of Jilin Province, China (20210101006JC).

Conflict of Interest

YZ was employed by Helong Forestry Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.919856/full#supplementary-material

References

Barik, S., SarkarDas, S., Singh, A., Gautam, V., Kumar, P., Majee, M., et al. (2014). Phylogenetic Analysis Reveals Conservation and Diversification of Micro RNA166 Genes Among Diverse Plant Speciesfication of Micro RNA166 Genes Among Diverse Plant Species. Genomics 103, 114–121. doi:10.1016/j.ygeno.2013.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Boualem, A., Laporte, P., Jovanovic, M., Laffont, C., Plet, J., Combier, J.-P., et al. (2008). MicroRNA166 Controls Root and Nodule Development in Medicago Truncatula. Plant J. 54, 876–887. doi:10.1111/j.1365-313x.2008.03448.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Budak, H., Kantar, M., Bulut, R., and Akpinar, B. A. (2015). Stress Responsive miRNAs and isomiRs in Cereals. Plant Sci. 235, 1–13. doi:10.1016/j.plantsci.2015.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Colle, M., Leisner, C. P., Wai, C. M., Ou, S., Bird, K. A., Wang, J., et al. (2019). Haplotype-phased Genome and Evolution of Phytonutrient Pathways of Tetraploid Blueberry. GigaScience 8, 1–15. doi:10.1093/gigascience/giz012

CrossRef Full Text | Google Scholar

Dai, X., Zhuang, Z., and Zhao, P. X. (2018). PsRNATarget: a Plant Small RNA Target Analysis Server (2017 Release). Nucleic Acids Res. 46, W49–W54. doi:10.1093/nar/gky316

PubMed Abstract | CrossRef Full Text | Google Scholar

Devi, S. J. S. R., Madhav, M. S., Kumar, G. R., Goel, A. K., Umakanth, B., Jahnavi, B., et al. (2013). Identification of Abiotic Stress miRNA Transcription Factor Binding Motifs (TFBMs) in Rice. Gene 531, 15–22. doi:10.1016/j.gene.2013.08.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, L., and Wang, Y. (2021). MicroRNAs in Woody Plants. Front. Plant Sci. 12, 686831. doi:10.3389/fpls.2021.686831

PubMed Abstract | CrossRef Full Text | Google Scholar

Gruber, A. R., Lorenz, R., Bernhart, S. H., Neuböck, R., and Hofacker, I. L. (2008). The Vienna RNA Websuite. Nucleic Acids Res. 36, W70–W74. doi:10.1093/nar/gkn188

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, V., Estrada, A. D., Blakley, I., Reid, R., Patel, K., Meyer, M. D., et al. (2015). RNA-seq Analysis and Annotation of a Draft Blueberry Genome Assembly Identifies Candidate Genes Involved in Fruit Ripening, Biosynthesis of Bioactive Compounds, and Stage-specific Alternative Splicing. GigaSci 4, 5. doi:10.1186/s13742-015-0046-9

CrossRef Full Text | Google Scholar

Hou, Y., Zhai, L., Li, X., Xue, Y., Wang, J., Yang, P., et al. (2017). Comparative Analysis of Fruit Ripening-Related miRNAs and Their Targets in Blueberry Using Small RNA and Degradome Sequencing. Ijms 18, 2767. doi:10.3390/ijms18122767

PubMed Abstract | CrossRef Full Text | Google Scholar

Kantar, M., Unver, T., and Budak, H. (2010). Regulation of Barley miRNAs upon Dehydration Stress Correlated with Target Gene Expression. Funct. Integr. Genomics 10, 493–507. doi:10.1007/s10142-010-0181-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kantar, M., Lucas, S. J., and Budak, H. (2011). MiRNA Expression Patterns of Triticum Dicoccoides in Response to Shock Drought Stress. Planta 233, 471–484. doi:10.1007/s00425-010-1309-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Karimi, M., Ghazanfari, F., Fadaei, A., Ahmadi, L., Shiran, B., Rabei, M., et al. (2016). The Small-RNA Profiles of Almond (Prunus Ducis Mill.) Reproductive Tissues in Response to Cold Stress. Plos One 11, e0156519. doi:10.1371/journal.pone.0156519

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitazumi, A., KawaharaOnda, Y. T. S., Onda, T. S., De Koeyer, D., and de los Reyes, B. G. (2015). Implications of miR166 and miR159 Induction to the Basal Response Mechanisms of an Andigena Potato (Solanum tuberosum Subsp. Andigena) to Salinity Stress, Predicted from Network Models in Arabidopsis. Genome 58, 13–24. doi:10.1139/gen-2015-0011

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuang, L., Yu, J., Shen, Q., Fu, L., and Wu, L. (2021). Identification of microRNAs Responding to Aluminium, Cadmium and Salt Stresses in Barley Roots. Plants 10, 2754. doi:10.3390/plants10122754

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S., Stecher, G., Li, M., Knyaz, C., and Tamura, K. (2018). MEGA X: Molecular Evolutionary Genetics Analysis across Computing Platforms. Mol. Biol. Evol. 35, 1547–1549. doi:10.1093/molbev/msy096

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuruvilla, L., Sathik, M., Luke, L. P., and Thomas, M. (2019). Identification and Validation of Drought-Responsive microRNAs from Hevea Brasiliensis. Acta Physiol. Plant. 41, 14. doi:10.1007/s11738-018-2803-8

CrossRef Full Text | Google Scholar

Li, J.-s., Fu, F.-l., An, M., Zhou, S.-f., She, Y.-h., and Li, W.-c. (2013). Differential Expression of microRNAs in Response to Drought Stress in Maize. J. Integr. Agric. 12, 1414–1422. doi:10.1016/S2095-3119(13)60311-1

CrossRef Full Text | Google Scholar

Li, X., Xie, X., Li, J., Cui, Y., Hou, Y., Zhai, L., et al. (2017). Conservation and Diversification of the miR166 Family in Soybean and Potential Roles of Newly Identified miR166s. BMC Plant Biol. 17, 32. doi:10.1186/s12870-017-0983-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, G., Wang, Y., Lou, X., Li, H., and Zhang, C. (2018a). Identification of Blueberry miRNAs and Their Targets Based on High-Throughput Sequencing and Degradome Analyses. Ijms 19, 983. doi:10.3390/ijms19040983

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Fan, T., Wang, T., Dominic, K., Hu, F., Liu, L., et al. (2018b). Characterization and Functional Analysis of miR166f in Drought Stress Tolerance in Mulberry (Morus Multicaulis). Mol. Breed. 38, 132. doi:10.1007/s11032-018-0886-y

CrossRef Full Text | Google Scholar

Li, X., Hou, Y., Zhang, L., Zhang, W., Quan, C., Cui, Y., et al. (2014a). Computation Identification of Conserved microRNAs and Their Targets from Expression Sequence Tags of Blueberry (Vaccinium Corybosum). Plant Signal. Behav. 9, e29462. doi:10.4161/psb.29462

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Wu, L.-Q., Zheng, W.-Y., Wang, R.-F., and Yang, L.-X. (2014b). Genome-wide Identification of microRNAs Responsive to High Temperature in Rice (Oryza Sativa) by High-Throughput Deep Sequencing. J. Agro Crop Sci. 201, 379–388. doi:10.1111/jac.12114

CrossRef Full Text | Google Scholar

Mathews, D. H., Disney, M. D., Childs, J. L., Schroeder, S. J., Zuker, M., and Turner, D. H. (2004). Incorporating Chemical Modification Constraints into a Dynamic Programming Algorithm for Prediction of RNA Secondary Structure. Proc. Natl. Acad. Sci. U.S.A. 101, 7287–7292. doi:10.1073/PANS.040179910110.1073/pnas.0401799101

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyashima, S., Honda, M., Hashimoto, K., Tatematsu, K., Hashimoto, T., Sato-Nara, K., et al. (2013). A Comprehensive Expression Analysis of the Arabidopsis MICRORNA165/6 Gene Family during Embryogenesis Reveals a Conserved Role in Meristem Specification and a Non-cell-autonomous Function. Plant Cell Physiol. 54, 375–384. doi:10.1093/pcp/pcs188

PubMed Abstract | CrossRef Full Text | Google Scholar

Ou, Y., Liu, X., Xie, C., Zhang, H., Lin, Y., Li, M., et al. (2015). Genome-Wide Identification of microRNAs and Their Targets in Cold-Stored Potato Tubers by Deep Sequencing and Degradome Analysis. Plant Mol. Biol. Rep. 33, 584–597. doi:10.1007/s11105-014-0771-8

CrossRef Full Text | Google Scholar

Pagliarani, C., Vitali, M., Ferrero, M., Vitulo, N., Incarbone, M., Lovisolo, C., et al. (2017). The Accumulation of miRNAs Differentially Modulated by Drought Stress Is Affected by Grafting in Grapevine. Plant Physiol. 173, 2180–2195. doi:10.1104/pp.16.01119

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, Y., Niu, M., Liang, J., Lin, E., Tong, Z., and Zhang, J. (2017). Idenfification of Heat-Responsive miRNAs to Reverl the miRNA-Mediated Regulatory Network of Heat Stress Response in Betula Luminifera. Trees 31, 1635–1652. doi:10.1007/s00468-017-1575-x

CrossRef Full Text | Google Scholar

Polashock, J., Zelzion, E., Fajardo, D., Zalapa, J., Georgi, L., Bhattacharya, D., et al. (2014). The American Cranberry: First Insights into the Whole Genome of a Species Adapted to Bog Habitat. BMC Plant Biol. 14, 165. doi:10.1186/1471-2229-14-165

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, A., Ruiz-Manriquez, L. M., Serrano-Cano, F. I., Reyes-Pérez, P. R., Tovar Alfaro, C. K., Barrón Andrade, Y. E., et al. (2020). Idenfification of microRNA and Their Expression in Leaf Tissues of Guava (Psidium Guajava L.) under Salinity Stress. Agronomy 10, 1920. doi:10.3390/agronomy10121920

CrossRef Full Text | Google Scholar

Si, J., Zhou, T., Bo, W., Xu, F., and Wu, R. (2014). Genome-Wide Analysis of Salt-Responsive and Novel microRNAs in Populus Euphratica by Deep Sequencing. BMC Genet. 15, S6. doi:10.1186/1471-2156-15-S1-S6

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, A., Singh, S., Panigrahi, K. C. S., Reski, R., and Sarkar, A. K. (2014). Balanced Activity of microRNA166/165 and its Target Transcripts from the Class III Homeodomain-Leucine Zipper Family Regulates Root Growth in Arabidopsis thaliana. Plant Cell Rep. 33, 945–953. doi:10.1007/s00299-014-1573-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Lin, L., and Sui, N. (2019). Regulation Mechanism of microRNA in Plant Response to Abiotic Stress and Breeding. Mol. Biol. Rep. 46, 1447–1457. doi:10.1007/s11033-018-4511-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Thakur, V., Wanchana, S., Xu, M., Bruskiewich, R., Quick, W. P., Mosig, A., et al. (2011). Characterization of Statistical Features for Plant microRNA Prediction. BMC Genomics 12, 108. doi:10.1186/1471-2164-12-108

PubMed Abstract | CrossRef Full Text | Google Scholar

Valiollahi, E., Farsi, M., and Kakhki, A. M. (2014). Sly-miR166 and Sly-miR319 Are Components of the Cold Stress Response in Solanum lycopersicum. Plant Biotechnol. Rep. 8, 349–356. doi:10.1007/s11816-014-0326-3

CrossRef Full Text | Google Scholar

Wang, R., Xu, L., Zhu, X., Zhai, L., Wang, Y., Yu, R., et al. (2015). Transcriptome-Wide Characterization of Novel and Heat-Stress-Responsive microRNAs in Radish (Raphanus Sativus L.) Using Next-Generation Sequencing. Plant Mol. Biol. Rep. 33, 867–880. doi:10.1007/s11105-014-0786-1

CrossRef Full Text | Google Scholar

Wang, Q., Yang, Y., Lu, G., Sun, X., Feng, Y., Yan, S., et al. (2020). Genome-wide Identification of microRNAs and Phased siRNAs in Soybean Roots under Long-Term Salt Stress. Genes Genom. 42, 1239–1249. doi:10.1007/s13258-020-00990-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, A., Kumar, S., Verma, R., Lata, C., Sanyal, I., and Rai, S. P. (2021). MicroRNA 166: an Evolutionarily Conserved Stress Biomarker in Land Plants Targeting HD-ZIP Family. Physiol. Mol. Biol. Plants 27, 2471–2485. doi:10.1007/s12298-021-01096-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z., Zhu, P., Kang, H., Liu, L., Cao, Q., Sun, J., et al. (2020). High-throughput Deep Sequencing Reveals the Important Role that microRNAs Play in the Salt Response in Sweet Potato (Ipomoea Batatas L.). BMC Genomics 21, 164. doi:10.1186/s12864-020-6567-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yue, J., Lu, X., Zhang, H., Ge, J., Gao, X., and Liu, Y. (2017). Identification of Conserved and Novel microRNAs in Blueberry. Front. Plant Sci. 8, 1155. doi:10.3389/fpls.2017.01155

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, X., Xu, Y., Jiang, J., Zhang, F., Ma, L., Wu, D., et al. (2018). Identification of Cold Stress Responsive microRNAs in Two Winter Turnip Rape (Brassica Rapa L.) by High Throughput Sequencing. BMC Plant Biol. 18, 52. doi:10.1186/s12870-018-1242-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Zhang, H., Srivastava, A. K., Pan, Y., Bai, J., Fang, J., et al. (2018). Knockdown of Rice microRNA166 Confers Drought Resistance by Causing Leaf Rolling and Altering Stem Xylem Development. Plant Physiol. 176, 2082–2094. doi:10.1104/pp.17.01432

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Wang, G., Sutoh, K., Zhu, J., and Zhang, W. (2008). Identification of Cold-Inducible microRNAs in Plants by Transcriptome Analysis. Biochim. Biophys. Acta (BBA) - Gene Regul. Mech. 1779, 780–788. doi:10.1016/j.bbagrm.2008.04.005

CrossRef Full Text | Google Scholar