Comparative Analysis of Fruit Ripening-Related miRNAs and Their Targets in Blueberry Using Small RNA and Degradome Sequencing

MicroRNAs (miRNAs) play vital roles in the regulation of fruit development and ripening. Blueberry is an important small berry fruit crop with economical and nutritional value. However, nothing is known about the miRNAs and their targets involved in blueberry fruit ripening. In this study, using high-throughput sequencing of small RNAs, 84 known miRNAs belonging to 28 families and 16 novel miRNAs were identified in white fruit (WF) and blue fruit (BF) libraries, which represent fruit ripening onset and in progress, respectively. Among them, 41 miRNAs were shown to be differentially expressed during fruit maturation, and 16 miRNAs representing 16 families were further chosen to validate the sRNA sequencing data by stem-loop qRT-PCR. Meanwhile, 178 targets were identified for 41 known and 7 novel miRNAs in WF and BF libraries using degradome sequencing, and targets of miR160 were validated using RLM-RACE (RNA Ligase-Mediated (RLM)-Rapid Amplification of cDNA Ends) approach. Moreover, the expression patterns of 6 miRNAs and their targets were examined during fruit development and ripening. Finally, integrative analysis of miRNAs and their targets revealed a complex miRNA-mRNA regulatory network involving a wide variety of biological processes. The findings will facilitate future investigations of the miRNA-mediated mechanisms that regulate fruit development and ripening in blueberry.


Introduction
Blueberry is an economically and nutritionally important small fruit crop, and commercially harvested blueberry species include lowbush (Vaccinium angustifolium), highbush (V. corymbosum), and rabbiteye bluberry (V. ashei or V. virgatum). Blueberry fruit not only tastes sweet with variable acidity, but also is an especially rich source of polyphenolic anthocyanin pigments [1], which contribute to recognized health-promoting effects including improvement of night vision, prevention of macular degeneration, anti-cancer activity, and reduced risk of heart disease [2]. For further genetic improvement of the crop, it is crucial to elucidate the molecular mechanisms that control fruit development and biosynthesis of anthocyanins in blueberry.
The developmental process of blueberry fruit after fruit ripening onset can be conceptually divided into two sequential phases: expansion and maturation. In the expansion phase, young fruit is hard and dark green with varied sizes. Subsequently, enlarged fruit begins to soften and accumulate red and then blue pigments in the maturation phase, which is further divided into different stages by fruit colors (light green/white, pink and blue stage) [3]. In recent years, RNA sequencing has become a powerful technology to profile the transcriptome and has provided much valuable information for gene identification and their potential roles during blueberry fruit development and maturation [3][4][5][6]. However, additional regulators need to be discovered to better understand the developmental regulatory network that controls fruit development and ripening.
MicroRNAs (miRNAs) are a class of endogenous, 20-24 nucleotides (nt) in length, single-stranded, non-coding RNAs which act as important regulatory molecules found in diverse eukaryotes [7]. Generally, MIR genes can be transcribed into primary transcripts (pri-miRNA), which are subsequently processed into the stem-loop precursor (pre-miRNA) by a DCL (Dicer-like) protein. The pre-miRNA is further processed by DCL1 into a duplex of miRNA and miRNA* [8]. In plants, miRNA can bind to transcripts of target genes based on perfect or near-perfect pairing, leading either to cleavage-induced degradation of the mRNA or translational inhibition [9,10].
Plant miRNAs act as master regulators implicated in regulating various cellular and metabolic processes [7,[11][12][13][14]. Recently, emerging evidence has demonstrated that miRNAs play an important role in the regulation of fruit development and maturation. For example, abnormal fruit morphologies can be generated by over-expression of AtMIR156b in tomato [15], while small tandem target mimic-mediated blockage of miR858 promoted anthocyanin accumulation in tomato [16]. miR156/157 and miR172 were shown to control the ripening process of tomato via targeting the known ripening regulators CNR and AP2a. Similarly, miR159 acts as a ripening regulator by targeting FaGAMYB, which plays crucial roles during the transition of strawberry receptacle from development to ripening [17]. Most recently, it was reported that miR73 was involved in regulation of strawberry fruit ripening via targeting ABI5 transcript to affect ABA (Abscisic acid) signaling pathway [18]. With the application of the high-throughput sequencing technology in miRNA-related research, numerous miRNAs involved in fruit development and maturation have been identified in many fruit-producing species such as apple [19], grape [20,21], peach [22,23], pomegranate [24], date palm [25], sweet orange [26]. Meanwhile, degradome sequencing has been developed for genome-wide identification of miRNA target genes and has been successfully applied to reveal miRNA-mRNA target pairs during fruit development [19,26]. Thus, the combinatorial analysis of miRNAome and degradome has greatly enhanced our understanding of the roles of miRNA during fruit development.
The blueberry genome is large (600 Mb/haploid genome) [27] and genomic information is limited. We previously reported the computational identification of conserved miRNAs and their targets from expression sequence tags in blueberry [28]. Recently, the availability of the genome sequence of cranberry (closely related to the species of blueberry) and blueberry transcriptome data makes it possible to effectively identify miRNAs and their targets in blueberry at genome level. In this study, a combinatorial analysis of miRNAome and degradome sequencing was conducted to identify and characterize the miRNAs and their targets that act during blueberry fruit maturation, a key developmental phase for fruit quality and biosynthesis of anthocyanin. The libraries for miRNAome and degradome profiling were made from fruits at white and blue fruit stages, which represent fruit ripening onset and in progress, respectively. The findings provide new insights into understanding the regulatory network that controls fruit ripening in blueberry.

Global Analysis of Sequences from Small RNA and Degradome Libraries
Fruit maturation is a very important phase in the developmental process of blueberry fruit, and many complex biochemical changes occur during the phase such as softening, sweetening, cell enlargement and pigmentation. To identify miRNAs involved in the ripening of blueberry fruits, two small RNA libraries were constructed for sequencing from samples of white (WF) and blue fruits (BF), which represent ripening onset and progression, respectively. Many biological changes involved in fruit ripening occur between these two stages such as fruit softening, coloration [6]. In total, approximately 13.89 and 31.40 million raw reads representing 3.08 and 5.99 million unique sequences were generated from WF and BF, respectively (Table S1). After removing the low quality reads, 3 adapters, junk reads, and redundant repeats, 4.58 and 9.56 million clean reads were obtained from the two libraries, including 1.53 and 2.38 Rfam RNAs (rRNA, tRNA, snoRNA and snRNA, Table S1), which showed high quality sequencing data in this study. The small RNA sequences with length of 18-25 nucleotides accounted for 66.66% (WF) and 75.09% (BF) of all the clean reads, respectively. Among them, 24 nt small RNAs are the most abundant in both libraries, followed by 21 nt ( Figure S1). Moreover, a statistical analysis of the nucleotide bias for all locations  showed that sRNAs starting with a uridine (U) at the 5 end are commonly present in blueberry fruit, and the U is also most frequently observed at the 24th nucleotide ( Figure S2).
To further demonstrate the biological functions of miRNAs identified in blueberry, degradome sequencing technology was adopted to perform transcriptome-wide analysis of the miRNA-cleaved mRNAs. A total of 10,686,450 and 13,702,581 raw reads were obtained from the two degradome libraries (WF and BF), respectively (Table S1). Approximately 4.49 and 4.96 million unique sequences were successfully mapped to the V. corymbosum reference sequences (transcriptome and expressed sequence tags), accounting for 99.31% of all unique reads (Table S1). Among them, more than 55% reads were mapped to cDNA sequences of blueberry. All mapped sequences were subsequently applied to detect candidate targets of miRNAs.

Identification of Fruit Ripening-Responsive miRNAs in Blueberry
To identify the miRNAs responsive to fruit ripening, the valid sequences of sRNAs from the libraries of white (WF) and blue fruits (BF) were aligned with the known plant miRNAs in miRBase 21.0 or the Vaccinium genome. Following strict screening criteria for annotating miRNA sequence, a total of 84 known miRNAs belonging to 28 families and 16 novel miRNAs were identified in blueberry fruit (Table 1). Among them, 70 known miRNAs (83.3%) and 16 novel miRNAs (100%) were detected in both libraries, whereas 13 miRNAs were detected only in BF library, and one miRNA in WF (Figure 1), suggesting their specific roles in the regulation of fruit ripening. small RNA libraries were constructed for sequencing from samples of white (WF) and blue fruits (BF), which represent ripening onset and progression, respectively. Many biological changes involved in fruit ripening occur between these two stages such as fruit softening, coloration [6]. In total, approximately 13.89 and 31.40 million raw reads representing 3.08 and 5.99 million unique sequences were generated from WF and BF, respectively (Table S1). After removing the low quality reads, 3′ adapters, junk reads, and redundant repeats, 4.58 and 9.56 million clean reads were obtained from the two libraries, including 1.53 and 2.38 Rfam RNAs (rRNA, tRNA, snoRNA and snRNA, Table S1), which showed high quality sequencing data in this study. The small RNA sequences with length of 18-25 nucleotides accounted for 66.66% (WF) and 75.09% (BF) of all the clean reads, respectively. Among them, 24 nt small RNAs are the most abundant in both libraries, followed by 21 nt ( Figure S1). Moreover, a statistical analysis of the nucleotide bias for all locations  showed that sRNAs starting with a uridine (U) at the 5′ end are commonly present in blueberry fruit, and the U is also most frequently observed at the 24th nucleotide ( Figure S2). To further demonstrate the biological functions of miRNAs identified in blueberry, degradome sequencing technology was adopted to perform transcriptome-wide analysis of the miRNA-cleaved mRNAs. A total of 10,686,450 and 13,702,581 raw reads were obtained from the two degradome libraries (WF and BF), respectively (Table S1). Approximately 4.49 and 4.96 million unique sequences were successfully mapped to the V. corymbosum reference sequences (transcriptome and expressed sequence tags), accounting for 99.31% of all unique reads (Table S1). Among them, more than 55% reads were mapped to cDNA sequences of blueberry. All mapped sequences were subsequently applied to detect candidate targets of miRNAs.

Identification of Fruit Ripening-Responsive miRNAs in Blueberry
To identify the miRNAs responsive to fruit ripening, the valid sequences of sRNAs from the libraries of white (WF) and blue fruits (BF) were aligned with the known plant miRNAs in miRBase 21.0 or the Vaccinium genome. Following strict screening criteria for annotating miRNA sequence, a total of 84 known miRNAs belonging to 28 families and 16 novel miRNAs were identified in blueberry fruit (Table 1). Among them, 70 known miRNAs (83.3%) and 16 novel miRNAs (100%) were detected in both libraries, whereas 13 miRNAs were detected only in BF library, and one miRNA in WF (Figure 1), suggesting their specific roles in the regulation of fruit ripening. The majority of the 84 known miRNAs are conserved across plant species (at least 12 species) based on the data from miRBase 21.0 such as miR156/157, miR159, miR160, miR165/166, miR172, miR319, miR396s, miR2111 and miR2118. Also, 7 less-conserved miRNA families were identified in blueberry fruit, including miR845, miR894, miR2603, miR6149, miR6478, miR7122 and miR7693, and the data from miRBase 21.0 indicated that they are present in only 1-5 plant species. Notably, the member numbers of each miRNA family varied in blueberry fruit. For example, miR156/157, miR166 and miR396 family contained 11, 10 and 13 members, respectively; while only one member was detected for 16 miRNAs families, including miR159, miR162, miR164, miR172, miR394, miR395, miR2111 and some less-conserved miRNA families (Table 1). The majority of the 84 known miRNAs are conserved across plant species (at least 12 species) based on the data from miRBase 21.0 such as miR156/157, miR159, miR160, miR165/166, miR172, miR319, miR396s, miR2111 and miR2118. Also, 7 less-conserved miRNA families were identified in blueberry fruit, including miR845, miR894, miR2603, miR6149, miR6478, miR7122 and miR7693, and the data from miRBase 21.0 indicated that they are present in only 1-5 plant species. Notably, the member numbers of each miRNA family varied in blueberry fruit. For example, miR156/157, miR166 and miR396 family contained 11, 10 and 13 members, respectively; while only one member was detected for 16 miRNAs families, including miR159, miR162, miR164, miR172, miR394, miR395, miR2111 and some less-conserved miRNA families (Table 1). Sixteen novel miRNA sequences derived from 14 loci were identified in both libraries, namely vco-miR_n01 to vco-miR_n14. The characteristics of all the 16 novel miRNAs are summarized in Table 1. The length of these novel miRNAs varied from 19 to 24 nt, with the majority being 21 nt long (62.5%), which is consistent with the previous reports in plant [29]. The sequences of all the novel miRNAs were located on the stem-arm of secondary hairpin structure of their precursors. Among them, 5 of the novel miRNAs (miR_n01-02, miR_n08-5p, miR_n09, miR_n12-3p) began with a 5 uridine, a characteristic feature of miRNAs. The precursors of the novel miRNAs ranged from 96 nt to 221 nt in length, which is acceptable since the lengths of miRNA precursors normally range from 55 to 930 nt in plants. They showed perfect or near-perfect stem-loop structures with no more than 2 mismatches in the miRNA/miRNA* duplex (Table 1 and Figure S3). It has been accepted that the minimal folding free energy index (MFEI) is an important criterion to differentiate the miRNA from other types of RNA [30]. In this study, the precursors of all the novel miRNA had perfect MFEIs (1.1-1.7) with an average of 1.4, which were higher than that of tRNAs (0.64), rRNAs (0.59) or mRNAs (0.62-0.66). The high MFEIs meet the requirements to maintain stability of secondary structure of miRNAs.

Phylogenetic Analysis of the Known and Novel miRNAs in Blueberry Fruit
To systematically exploit the evolutionary conservation of the known miRNAs in plant species, a phylogeny tree was constructed using the sequences of mature miRNAs in blueberry and grape, as they have a lot of commonalities in fruit property and developmental process. As shown in Figure 2, 84 blueberry miRNAs and 105 grape miRNAs were grouped into 10 classes, representing similarities and divergence. For example, all the members of the 9 miRNA families (miR395, miR166, miR408, miR2118, miR393, miR403, miR2603, miR167, miR894) from both blueberry and grape were clustered into class X, suggesting they were evolutionarily closer to each other. It was also observed that blueberry miRNAs were basically grouped together with their homologues in grape. For instance, all the vco-miR396-5ps were clustered together with vvi-miR396a/b/c/d in class IX, while all the vco-miRNA156/157-5ps and vvi-miR156a/b/f/g/h/i formed a discrete clade in the phylogenetic tree. These observations indicate that vco-miRNAs have close relationship with their corresponding homologs in grape. Meanwhile, the topological tree of miRNA precursors in blueberry and grape were constructed. In total, 55 and 93 precursors in blueberry and grape, respectively, were clustered into 5 classes, reflecting both conservation and divergence ( Figure S4). As expected, miRNA precursors in blueberry showed close evolutionary relationship with their homologs in grape, further supporting the conservation of blueberry miRNAs.
To analyze the relationship of the newly identified vco-miRNAs with other miRNAs, their sequences were used as queries to search against all the miRNA sequences in miRBase 21 including animals, viruses and plants. As a result, one novel miRNA (vco-miR_n03) showed near-perfect sequence match with miR9202 (5 mismatches) from Eptesicus fuscus. Furthermore, the precursors of the 16 novel miRNAs were utilized to perform a BLASTX against the nucleotide sequence database at NCBI. The analysis revealed that homologs of all the precursors, except pre-vco-miR_n10, could not be found in any plant species other than those in the Vaccinium genus, suggesting that they might be Vaccinium-specific miRNAs in plants. An uncharacterized ncRNA in Jatropha curcas was found to be homologous to vco-miR_n10 with 82% identity.

Accumulation Patterns of Blueberry miRNAs during Fruit Ripening
It is generally accepted that the number of miRNA reads in a library can be considered as an indicator for estimating the relative abundance of miRNAs. In the present study, the normalized read counts of the known miRNAs varied from 0.42 to 1377.66 in BF library and 0 to 1055.95 in WF library. vco-miR166b-3p and vco-miR894 were the most abundant in BF and WF libraries, respectively (Table 1). Interestingly, different members in the same miRNA family also showed considerably variable abundance. In miR396 family, for example, the read number of miR396b-3p was 60.77 in BF library, whereas only 3.33 reads were detected for miR396c-3p. The variable abundance of miRNAs in each family might reflect their functional divergence during fruit maturation.
To provide clues for their roles during fruit maturation, the abundance of the miRNAs in the two libraries were analyzed and compared using the normalized reads from the sRNA sequencing data. As shown in Figure 3, 100 miRNAs were categorized into three groups based on their accumulation patterns. Group II comprised 32 miRNAs with down-regulated accumulation as fruit ripens, whereas the rest of the miRNAs were up-regulated, which were separated into group I (38 known and 3 novel miRNAs) and III (23 known and 4 novel miRNAs) ( Figure 3). It was also observed that different members within each miRNA family were conversely accumulated in the two libraries. Phylogenetic analysis of known miRNA mature sequences of blueberry (Vaccinium corymbosum) and grape (Vitis vinifera) obtained from the sequencing data (this work) and miRbase 21.0, respectively. The miRNAs were clustered into 10 classes (I-X). The miRNAs in light blue background were from blueberry and the ones in dark grey were from grape in the inner ring. The middle ring shows the miRNA families in different colors. The outer ring shows the length of mature miRNAs.

Accumulation Patterns of Blueberry miRNAs during Fruit Ripening
It is generally accepted that the number of miRNA reads in a library can be considered as an indicator for estimating the relative abundance of miRNAs. In the present study, the normalized read counts of the known miRNAs varied from 0.42 to 1377.66 in BF library and 0 to 1055.95 in WF library. vco-miR166b-3p and vco-miR894 were the most abundant in BF and WF libraries, respectively (Table 1). Interestingly, different members in the same miRNA family also showed considerably variable abundance. In miR396 family, for example, the read number of miR396b-3p was 60.77 in BF library, whereas only 3.33 reads were detected for miR396c-3p. The variable abundance of miRNAs in each family might reflect their functional divergence during fruit maturation.
To provide clues for their roles during fruit maturation, the abundance of the miRNAs in the two libraries were analyzed and compared using the normalized reads from the sRNA sequencing data. As shown in Figure 3, 100 miRNAs were categorized into three groups based on their accumulation patterns. Group II comprised 32 miRNAs with down-regulated accumulation as fruit ripens, whereas the rest of the miRNAs were up-regulated, which were separated into group I (38 known and 3 novel miRNAs) and III (23 known and 4 novel miRNAs) ( Figure 3). It was also observed that different members within each miRNA family were conversely accumulated in the two libraries. For example, the accumulation of miR166a in BF library was increased by 3.0 folds, whereas miR166c-3p was decreased by 1.8 folds as compared to the one in WF library. These observations suggest that functional diversity of miRNAs might occur during fruit maturation, even for members within the same miRNA family. For example, the accumulation of miR166a in BF library was increased by 3.0 folds, whereas miR166c-3p was decreased by 1.8 folds as compared to the one in WF library. These observations suggest that functional diversity of miRNAs might occur during fruit maturation, even for members within the same miRNA family.  To gain deeper insight into the roles of miRNAs during fruit maturation, we performed differential expression pattern analysis of known and novel miRNAs between the WF and BF libraries. The miRNAs were considered as up-regulated or down-regulated ones between the two libraries with the criteria of |log 2 (fold change)| ≥ 1 and p-value ≤ 0.05. A total of 35 known miRNAs belonging to 16 miRNA families and 6 novel miRNAs were differentially expressed during fruit maturation ( Table 2). The majority of these miRNAs were up-regulated as fruit ripens. Among them, 11 miRNAs (vco-miR156c-3p/f, vco-miR159a, vco-miR160a-5p, vco-miR171a-3p/a-5p/b, vco-miR394a-5p, vco-miR398a-3p/b-3p/b-5p) were detected only in the BF library. Additionally, 3 miRNAs (vco-miR166b-5p, vco-miR319b and vco-miR7693-5p) was dramatically altered with 64.44, 16.82 and 13.86 fold changes, respectively. The up-regulation of miRNAs suggests that they might have positive correlation with fruit maturation via targeting repressor(s) of fruit ripening. Meanwhile, two miRNAs (vco-miR894 and vco-miR_n12-3p) were found to be down-regulated as fruit ripens, implying that they might be negatively responsive to fruit maturation.  Figure 4, the change tendency for each of the 16 miRNAs during ripening was quite consistent with the data from sRNA sequencing. The expression levels of the 11 up-regulated. miRNAs (according to the sRNA sequencing data) were all increased, while the other 5 miRNAs (miR164a, miR894, miR_n01, miR_n03, miR_n12-5p) were shown to be down-regulated by both qRT-PCR and high-throughput sequencing. These observations indicated that the sRNA sequencing data presented here are highly reliable, and they could represent relative expression levels of the identified miRNAs in both blueberry libraries.
To validate the sRNA sequencing data, 16 miRNAs representing 16 families were randomly selected for analysis of their expression patterns using stem-loop qRT-PCR. The change tendency was calculated with log2(BF/WF), where BF/WF refers to the ratio of miRNA abundance in BF and WF libraries. As shown in Figure 4, the change tendency for each of the 16 miRNAs during ripening was quite consistent with the data from sRNA sequencing. The expression levels of the 11 up-regulated. miRNAs (according to the sRNA sequencing data) were all increased, while the other 5 miRNAs (miR164a, miR894, miR_n01, miR_n03, miR_n12-5p) were shown to be down-regulated by both qRT-PCR and high-throughput sequencing. These observations indicated that the sRNA sequencing data presented here are highly reliable, and they could represent relative expression levels of the identified miRNAs in both blueberry libraries.

Identification of miRNA Targets during Fruit Ripening by Degradome Sequencing
To understand the biological functions of miRNAs, it is essential to identify their targets. Therefore, a high-throughput degradome sequencing technology was employed to perform a transcriptome-wide analysis of the miRNA-guided cleavage of target mRNAs. As a result, a total of 178 sliced targets for 41 known and 7 novel miRNAs were identified in WF and BF libraries according to the CleaveLand pipeline (Tables S2 and S3). Furthermore, a confidence evaluation of the degradome data was conducted as described previously [29], and the sliced targets were grouped into 5 categories (0-4) accordingly. As exemplified in Figure 5, the abundance of miRNA-guided cleavage fragment is equal to the maximum on the transcript in category 0 and 1, while the cleavage fragment was not the most abundant tag, but it still formed a clear peak in category 2 and 3. Category 4 refers to only 1 raw read of miRNA-guided cleavage fragment detected in degradome sequencing.

Identification of miRNA Targets during Fruit Ripening by Degradome Sequencing
To understand the biological functions of miRNAs, it is essential to identify their targets. Therefore, a high-throughput degradome sequencing technology was employed to perform a transcriptome-wide analysis of the miRNA-guided cleavage of target mRNAs. As a result, a total of 178 sliced targets for 41 known and 7 novel miRNAs were identified in WF and BF libraries according to the CleaveLand pipeline (Tables S2 and S3). Furthermore, a confidence evaluation of the degradome data was conducted as described previously [29], and the sliced targets were grouped into 5 categories (0-4) accordingly. As exemplified in Figure 5, the abundance of miRNA-guided cleavage fragment is equal to the maximum on the transcript in category 0 and 1, while the cleavage fragment was not the most abundant tag, but it still formed a clear peak in category 2 and 3. Category 4 refers to only 1 raw read of miRNA-guided cleavage fragment detected in degradome sequencing. Category 0 is the most reliable for the detection of miRNA target genes, and 30 targets were included such as  (Tables S2 and S3). Meanwhile, 4 targets were detected in category 1, 90 targets in category 2, 10 targets in category 3 and 43 targets in category 4. Further observation indicated that the majority of the cleavage sites were distributed at the position between nucleotides 9 and 11 of the miRNA sequences from the 5 to the 3 end. Moreover, two cleavage targets of miR160b were selected for verification by RLM-5 RACE (RNA Ligase-Mediated (RLM)-Rapid Amplification of cDNA Ends). As expected, the same cleavage site was detected for these two targets ( Figure S5A,F), indicating a high reliability of the degradome data. 9 and 11 of the miRNA sequences from the 5′ to the 3′ end. Moreover, two cleavage targets of miR160b were selected for verification by RLM-5′ RACE (RNA Ligase-Mediated (RLM)-Rapid Amplification of cDNA Ends). As expected, the same cleavage site was detected for these two targets ( Figure S5A,F), indicating a high reliability of the degradome data. To better understand the functional roles of the miRNAs during fruit maturation, 92 cleaved transcripts caused by the differentially accumulated miRNAs in the two libraries were subjected to GO (Gene Ontology) functional analysis using Blast2GO. 57 targets were successfully assigned into 129 GO terms at the second level. By filtering out low abundance transcripts (copy numbers less than 2), the targets were shown to be enriched in 21 biological processes, 13 cellular components and 16 molecular functions ( Figure 5). In biological processes, the predominant terms were related to protein folding (GO:0006457), regulation of transcription (GO:0006355), transcription (GO:0006351), and auxinmediated signaling pathway (GO:0010928). Two GO terms, nucleus (GO:0005634) and cytoplasm (GO:0005737), were found to be predominant in cellular component category, while four terms were highly enriched in the molecular function including peptide binding, peptidyl-prolyl cis-trans isomerase activity, DNA binding, ATP binding ( Figure 5). This analysis indicates that the targets identified by degradome sequencing might be involved in diverse biological processes during fruit maturation, especially DNA-dependent transcription, transcription regulation and protein folding.
To gain insights into their functions, the miRNA targets detected by degradome sequencing were subjected to a BLASTX analysis to search for their homologs in plants. It  Furthermore, an integrative analysis of miRNAs and their targets was conducted to reveal the miRNA-mRNA regulatory network involved in blueberry fruit ripening. As shown in Figure 6, a number of To better understand the functional roles of the miRNAs during fruit maturation, 92 cleaved transcripts caused by the differentially accumulated miRNAs in the two libraries were subjected to GO (Gene Ontology) functional analysis using Blast2GO. 57 targets were successfully assigned into 129 GO terms at the second level. By filtering out low abundance transcripts (copy numbers less than 2), the targets were shown to be enriched in 21 biological processes, 13 cellular components and 16 molecular functions ( Figure 5). In biological processes, the predominant terms were related to protein folding (GO:0006457), regulation of transcription (GO:0006355), transcription (GO:0006351), and auxin-mediated signaling pathway (GO:0010928). Two GO terms, nucleus (GO:0005634) and cytoplasm (GO:0005737), were found to be predominant in cellular component category, while four terms were highly enriched in the molecular function including peptide binding, peptidyl-prolyl cis-trans isomerase activity, DNA binding, ATP binding ( Figure 5). This analysis indicates that the targets identified by degradome sequencing might be involved in diverse biological processes during fruit maturation, especially DNA-dependent transcription, transcription regulation and protein folding.
To gain insights into their functions, the miRNA targets detected by degradome sequencing were subjected to a BLASTX analysis to search for their homologs in plants. As shown in Figure 6, a number of conserved interactions between miRNAs and their targets across plant species were observed in blueberry such as vco-miR156-SPLs, vco-miR159-GAMYB, vco-miR172-AP2, vco-miR396-GRFs, vco-miR160-ARFs, vco-miR164-NAC, vco-miR403-AGO2 and vco-miR319-TIR1. Besides, novel interactions were also identified in blueberry fruit. For example, miR159, targeting GAMYB in many plant species, could also target RER4 and MOT2 in blueberry, while NRT1/PRT and GLD (GLUTAREDOXIN-C10) were found to be targeted by the less-conserved miR7693. Intriguingly, 8 novel interactions were observed for the 7 novel miRNAs, which are involved in different biological processes. These novel interactions might be involved in the maintenance or regulation of blueberry fruit development in a species-specific manner. It was also observed that multiple genes were cleaved by a single miRNA. For instance, miR396 family (miR396a-g) can target 33 unique transcripts belonging to 12 different gene families. This observation suggested that the miRNAs might play diverse roles in the regulation of fruit development and maturation via targeting multiple genes in blueberry. conserved interactions between miRNAs and their targets across plant species were observed in blueberry such as vco-miR156-SPLs, vco-miR159-GAMYB, vco-miR172-AP2, vco-miR396-GRFs, vco-miR160-ARFs, vco-miR164-NAC, vco-miR403-AGO2 and vco-miR319-TIR1. Besides, novel interactions were also identified in blueberry fruit. For example, miR159, targeting GAMYB in many plant species, could also target RER4 and MOT2 in blueberry, while NRT1/PRT and GLD (GLUTAREDOXIN-C10) were found to be targeted by the less-conserved miR7693. Intriguingly, 8 novel interactions were observed for the 7 novel miRNAs, which are involved in different biological processes. These novel interactions might be involved in the maintenance or regulation of blueberry fruit development in a species-specific manner. It was also observed that multiple genes were cleaved by a single miRNA. For instance, miR396 family (miR396a-g) can target 33 unique transcripts belonging to 12 different gene families. This observation suggested that the miRNAs might play diverse roles in the regulation of fruit development and maturation via targeting multiple genes in blueberry.

Dynamic Expression Patterns of Fruit-Related miRNAs and Their Targets during Fruit Development and Maturation
To investigate the dynamic expression patterns of the miRNAs at different developmental stages (green pad, green cup, white, pink and blue fruit) in blueberry, five known (vco-miR160b, vco-miR172a-3p, vco-miR319b, vco-miR396a-5p and vco-miR403a) and one novel (vco-miR_n10) vco-miRNAs were selected for stem-loop qRT-PCR analysis. As shown in Figure 7, these six miRNAs showed differential expression patterns during fruit development and maturation. miR160a was gradually increased as fruit develops and ripens, while 3 miRNAs (miR172a-3p, miR396a-5p and vco-miR_n10) were initially down-regulated, and remarkably elevated at blue fruit stage. miR319b and miR403a showed a stable accumulation at the initial four stages, and then a pronounced accumulation was clearly observed at blue fruit stage. Obviously, all the 6 miRNAs exhibited the highest accumulation at blue fruit stage, suggesting they might play crucial roles in the regulation of fruit maturation.

Dynamic Expression Patterns of Fruit-Related miRNAs and Their Targets during Fruit Development and Maturation
To investigate the dynamic expression patterns of the miRNAs at different developmental stages (green pad, green cup, white, pink and blue fruit) in blueberry, five known (vco-miR160b, vco-miR172a-3p, vco-miR319b, vco-miR396a-5p and vco-miR403a) and one novel (vco-miR_n10) vco-miRNAs were selected for stem-loop qRT-PCR analysis. As shown in Figure 7, these six miRNAs showed differential expression patterns during fruit development and maturation. miR160a was gradually increased as fruit develops and ripens, while 3 miRNAs (miR172a-3p, miR396a-5p and vco-miR_n10) were initially down-regulated, and remarkably elevated at blue fruit stage. miR319b and miR403a showed a stable accumulation at the initial four stages, and then a pronounced accumulation was clearly observed at blue fruit stage. Obviously, all the 6 miRNAs exhibited the highest accumulation at blue fruit stage, suggesting they might play crucial roles in the regulation of fruit maturation. , vco-miR172a-3p and its target gene (B); vco-miR319b and its target gene (C); vco-miR396a-5p and 3 target genes (D); vco-miR403a and its target gene (E); vco-miR_n10 and its target gene (F). The expression level in GP was set as 1. U6 and ACTIN genes were used as the internal control for miRNA expression and target genes expression, respectively. Error bars indicate standard error of three biological and technical replicates. The different letters (a-d) indicate significant differences at p < 0.01 according to Duncan's multiple range tests. GP, green pad; GC, green cup; WF, white fruit; PF, pink fruit; BF, blue fruit.
To confirm the regulatory roles of the vco-miRNAs, the expression patterns of 9 target genes corresponding to the 6 miRNAs were also investigated by qRT-PCR at the same developmental stages as above (Figure 7). For most target genes, their expressions were negatively correlated with those of their corresponding miRNAs. As shown in Figure 7, 8 target genes showed high expression at two fruit developmental stages (green pad and green cup stages), and decreased from onset of fruit Figure 7. Expression patterns of the known and novel vco-miRNAs identified in this work and their target genes during fruit development and ripening. The accumulation patterns of vco-miR160b and 2 target genes (A), vco-miR172a-3p and its target gene (B); vco-miR319b and its target gene (C); vco-miR396a-5p and 3 target genes (D); vco-miR403a and its target gene (E); vco-miR_n10 and its target gene (F). The expression level in GP was set as 1. U6 and ACTIN genes were used as the internal control for miRNA expression and target genes expression, respectively. Error bars indicate standard error of three biological and technical replicates. The different letters (a-d) indicate significant differences at p < 0.01 according to Duncan's multiple range tests. GP, green pad; GC, green cup; WF, white fruit; PF, pink fruit; BF, blue fruit.
To confirm the regulatory roles of the vco-miRNAs, the expression patterns of 9 target genes corresponding to the 6 miRNAs were also investigated by qRT-PCR at the same developmental stages as above (Figure 7). For most target genes, their expressions were negatively correlated with those of their corresponding miRNAs. As shown in Figure 7, 8 target genes showed high expression at two fruit developmental stages (green pad and green cup stages), and decreased from onset of fruit ripening (white stage) until blue fruit stage. Especially, the expressions of miR396-targeted GRF1 and miR160-targeted ARF18 were decreased at white stage with 19.4 and 2.6 fold changes ( Figure 7A,D), respectively, as compared to the ones at green cup stage, implying that they might play a role in the initiation of fruit maturation. It was also noteworthy that miR403-targeted AGO2 was decreased at pink and blue fruit stages by 7.9 and 5.9 folds ( Figure 7E), respectively, as compared to the ones at white stage, suggesting that it might be involved in the regulation of fruit maturation. The apparently negative correlation between the miRNAs and their target genes indicates that the miRNAs exert considerable influence on the expression of their target genes. Unexpectedly, ARF17 showed a stable expression pattern during fruit development and maturation except at green cup stage with a slight increase ( Figure 7A). It is possible that ARF17 and miR160 might be differentially distributed in blueberry fruit or cells; alternatively, it could also be possible that the expression of ARF17 could be mainly regulated by other factors, such as competing endogenous RNAs (ceRNAs), transcription factors and epigenetic regulators.

Discussion
Many studies have implicated the potential roles of miRNAs in fruit development and biosynthesis of anthocyanin [15][16][17]. Blueberries are fleshy berry fruits recognized for various phenolic compounds, especially anthocyanin pigments, which are involved in a wide variety of health benefits. However, it is totally unknown about the miRNAs and their targets involved in blueberry fruit development and maturation. In this study, genome-wide identification of miRNAs and their target genes involved in blueberry fruit ripening were conducted using small RNA and degradome sequencing, and integrative analysis suggested the miRNA-mediated regulatory networks governing blueberry fruit ripening.

miRNAs and Their Targets during Blueberry Fruit Maturation
Numerous fruit development-associated miRNAs have been identified in fruit-producing plants such as apple [19], grape [20,21], citrus [26], peach [23]. In this study, 84 known miRNAs belonging to 28 families and 16 novel miRNAs were identified in blueberry fruit. 41 miRNAs were differentially expressed between white fruit stage and blue fruit stage ( Table 2), suggesting that they might be involved in controlling the cellular and developmental processes during fruit maturation in blueberry. By comparing vco-miRNAs and the ones from other fruit-producing plants such as apple [19], grape [20,21], citrus [26], peach [23], tomato [31] and cucurbits [32], it was found that 21 miRNA families are more common across these fruit-producing plants, suggesting that they might play crucial and fundamental roles in fruit development and maturation in diverse plants. However, 7 miRNAs (such as miR2118, miR2603, miR6149, miR6478, miR7122, miR7693 and miR845) were found only in blueberry fruits. Therefore, we speculate that these 7 miRNAs, together with the 16 novel miRNAs, might be species-specific for blueberry fruit maturation.
It has been proposed that plant miRNAs can mediate gene expression mainly via miRNA-guided cleavage of target transcripts [9]. Therefore, the identification of miRNA targets is a key step to reveal the regulatory roles of miRNAs. In this study, 178 targets were identified for 41 known and 7 novel vco-miRNAs (Tables S2 and S3). Increasing studies indicated that miRNAs seemingly prefer targeting transcription factors in plants [33]. Based on GO analysis, the targets might be mainly involved in DNA-dependent transcription, transcription regulation and protein folding during fruit maturation, which provided preliminary clues for further studying the functions of miRNA-target modules. As have been shown in other species [34,35], the targets of the conserved miRNAs in blueberry are clearly conserved. More importantly, some miRNAs were found to have additional or novel targets in blueberry. For example, miR396s are well-known to target GRF family in many plant species [36]. In this study, 33 targets were detected for miR396s in blueberry. Additionally, 8 targets were detected for the 7 novel miRNAs, which are involved in various biological processes. We speculate that these new target genes and the novel miRNAs are likely acting in a species-specific manner during blueberry fruit maturation. Also, it was observed that some miRNA-target pairs showed weak correlation between their expressions. It is not surprising since gene expression is spatio-temporally regulated by a series of complicated mechanisms. The best-known mechanism is that miRNA can sequestered by ceRNAs, therefore leading to deregulation of target gene [37]. Overall, these findings will facilitate our understanding of the miRNA-mediated regulatory network during fruit maturation.

miRNA-Mediated Regulation Network during Blueberry Fruit Maturation
It has been well documented that auxin is not only the most important hormone in the regulation of cell expansion, but also essential for triggering the onset of ripening in tomato and grape berry at low levels [38,39]. In blueberry, it was found that IAA concentration was sharply increased between the S4 and S5 stages, and then began to decline as ripening proceeded [6]. However, little is known about auxin signaling and transport in ripe fruits. In this study, 3 conserved miRNA families (miR160, miR393 and miR396) were identified in blueberry fruit. Based on our degradome data, they may target some auxin regulators such as ARF17, ARF18, AFB2, TIR1 and ELONGATOR PROTEIN 6 (ELP6) (Figure 7). It has been proposed that ARFs serve as transcription factors that regulate the activation or repression of auxin-responsive genes, and their transcriptional activity can be depressed by heterotypic dimerization with AUX/IAA (AUXIN/INDOLEACETIC ACID) proteins [40]. TIR1 and AFB2 act as auxin receptors to facilitate the degradation of AUX/IAA transcriptional repressor protein, therefore releasing ARFs to activate the transcription of auxin-responsive genes [41]. Likewise, ELP6 may be involved in auxin transport through post-transcriptional modulation of PIN protein levels [42]. Our miRNA sequencing data indicated that the accumulations of miR160, miR393 and miR396 were clearly up-regulated as fruit ripens (Table 1), implying that depression of auxin signaling pathway might occur during fruit ripening in blueberry. More importantly, when coordinated expressions of miR160-ARF18 and miR396-ELP6 were investigated during fruit development and ripening, it was observed that the accumulations of miR160b and miR396a-5p were increased as blueberry fruit develops and ripens ( Figure 7A,D), while ARF18 and ELP6 were conversely expressed ( Figure 7A,D). These observations suggest that these miRNAs play roles in controlling the expression of auxin responsive genes during fruit development and maturation in blueberry. Since auxin delays the ripening-related processes such as sugar accumulation, and anthocyanin content [43], we propose that miR160, miR393 and mi396 might promote these ripening-related processes via suppressing the auxin signaling pathway. In addition, players that might connect different hormones and their signaling pathways were also detected in the study. For example, miR159 was detected to target GAMYB, which is not only responsive to GAs, but also involved in ABA and sucrose biosynthesis [44]; and NRT1/PTR, transporter for the plant hormones auxin, ABA, and GA [45], was targeted by a less conserved miRNA (miR7693). These results indicate the possible roles of these miRNAs in controlling fruit development and maturation through mediating hormonal crosstalk.
Plant organs require a proper and precise control of cell proliferation, expansion and differentiation to achieve their final size and shape. Increasing evidences have indicated that miRNAs play pivotal roles in regulating cell proliferation and expansion. For example, miR319-targeted TCPs can promote cell expansion and repress cell proliferation in Arabidopsis [46], while miR396-targeted CsGRFs may play important roles in regulating the growth of fruit and leaf in citrus [47]. Likewise, transposon insertional allele of miRNA172 showed large fruit phenotype, whereas over-expression of miRNA172 significantly reduced fruit size in apple [48]. In this study, the conserved interactions (miR172-AP2, miR319-TCP and miR396-GFR1) were also found in blueberry fruit ( Figure 6 and Table S2). Our study also revealed a novel interaction, miR396-CYP1 (Figure 7 and Table S2), which is very interesting since CYP1 in tomato (LeCYP1) has been reported to be involved in regulation of early fruit development growth [49,50]. Among the 6 stages of fruit development and ripening, the accumulation of all the 3 miRNAs were the highest at blue fruit stage (Table 1 and Figure 7B-D), while their corresponding targets (AP2, TCP, GRF1 and CYP1) were significantly decreased from onset of fruit ripening (white fruit stage), and reached their minimal level at blue fruit stage ( Figure 7B-D). The negative correlation between the expression patterns of these miRNAs and their targets are consistent with a scenario that cell enlargement of blueberry fruit can last until the early stage of fruit maturation. However, the expression patterns between these miRNAs and their targets were not exactly negatively correlated with each other. We speculated that some factors other than miRNAs might also be involved in the regulation of these genes such as transcription factors and epigenetic regulators. Recently, some evidences indicated that miRNAs exquisitely exert dual roles in the regulation of fruit development and maturation. For example, miRNA172 can be involved in fruit growth through negatively modulating cell division and expansion [48], while miR172-targeted SlAP2a is a negative regulator of ethylene production during ripening in tomato [51], suggesting that the interaction of miR172-AP2 might depress fruit development and promote fruit maturation. In this study, AP2, TCP, GRF and CYP1 were significantly depressed from the onset of fruit maturation until ripening ( Figure 7B-D), implying that these genes can promote fruit growth, but delay fruit maturation. We propose that miRNAs possibly act as crucial regulators to balance the expression of these genes, therefore satisfying the requirement of fruit development and maturation.
Many complex molecular changes occur within fruit during maturation. Protein ubiquitination leads to diverse fates of the target proteins, such as protein degradation, alteration in protein trafficking [52]. We identified two ubiquitination-related interactions such as miR396-UCTH24 (UBIQUITIN CARBOXYL-TERMINAL HYDROLASE 24), miR6478-UCTH9 and miR399-UBC24 ( Figure 6 and Table S2). It has been reported that UCTHs play important roles in different cellular processes including cell growth and differentiation through regulating cellular ubiquitin levels [53], while UBC24 is directly involved in the degradation of ubiquitinated proteins [54]. Differential accumulation of these miRNAs indicated that miRNAs might be involved in the regulation of protein degradation via ubiquitination-related pathway during fruit maturation. Meanwhile, some interactions associated with RNA processing events were identified such as miR403-AGO2, miR845-RBP2I (RNA POLYMERASE II SECOND LARGEST SUBUNIT) and miR7122-PPR ( Figure 6 and Table S2). Notably, AGO2, which was found to be targeted by miR403 in several plant species, is associated with small RNA function [55,56], while PPR proteins are known to mediate specific RNA processing events including RNA editing, transcript processing and translation initiation [57]. These results suggest that diverse events might occur, at both protein and RNA levels, in order to meet the requirement of fruit maturation in blueberry.
Cell wall structure is associated with fruit softening during maturation, and the biochemical characterization of cell wall polymers and their progressive modification during ripening has been documented in blueberry fruits [58]. In this study, two miRNAs (miR166 and miR894) appear to be associated with the modification and formation of cell wall structure ( Figure 6 and Table S2). Generally, miR166 can be involved in many developmental processes including seed maturation, leaf polarity determination and secondary wall biosynthesis, via targeting HD-ZIP III family genes [14,59,60]. During blueberry fruit maturation, however, miR166 was detected to target LEUCINE-RICH REPEAT EXTENSIN-LIKE PROTEIN (LRX) ( Figure 6 and Table S2), which might be involved in either the regulation of enzymatic activities in the cell wall or the recruitment of enzymes to the appropriate location in the cell wall, therefore altering the structure of cell wall [61]. Meanwhile, CAFFEOYL SHIKIMATE ESTERASE (CSE) targeted by miR894 ( Figure 6 and Table S2), serves as an enzyme central to the biosynthetic pathway of lignin [62], which is a major component of plant secondary cell walls. These findings might be useful for guiding future work to understand the control of fruit firmness and softening.
Fruit coloration during ripening can be achieved by chlorophyll breakdown and production of color metabolites such as anthocyanins and flavonoids. In this study, several miRNA-target pairs appear to be involved in the regulation of fruit coloration in blueberry. For example, miR156, a highly conserved miRNA family, is implicated in the regulation of anthocyanin biosynthesis and fruit development by targeting the SPL family genes [63]. In this study, 11 miR156/157 species were identified in blueberry fruit, and most of them showed up-regulated accumulation patterns during fruit ripening ( Table 1). As was shown before, SPL family was found to be target of miR156/157 in blueberry fruit (Table S3). The accumulation pattern of miR156s during fruit ripening is consistent with the requirement of fruit coloration in blueberry. Furthermore, miR396 and miR_n10 were detected to target FtsZs and chloroplastic BAG1, respectively (Table S3). FtsZs are a class of key players in chloroplast division. It has been reported that the expression of FtsZ in fruit was consistent with the chlorophyll concentration in tomato hp3 mutant [64]. Chloroplastic BAG1 regulates diverse cellular pathways in chloroplast, such as programmed cell death and stress responses [65]. Therefore, the up-regulation of miR396 and miR_n10 during fruit ripening led us to speculate that these miRNA-target modules might be involved in the regulation of fruit maturation via depression of chloroplast division and/or maintenance, therefore facilitating coloration in blueberry.
In conclusion, numerous miRNAs and their targets during blueberry fruit maturation were identified by next-generation sequencing approach in the study. The interplays between these miRNAs and their targets depict a picture of the regulatory network involved in blueberry fruit development and maturation, and some miRNAs might play vital roles in fruit maturation and coloration via repressing the expression of auxin-responsive genes or the ones responsible for fruit growth and enlargement ( Figure 8). These findings will facilitate future investigation of miRNA-mediated mechanisms underlying blueberry fruit development and maturation. players in chloroplast division. It has been reported that the expression of FtsZ in fruit was consistent with the chlorophyll concentration in tomato hp3 mutant [64]. Chloroplastic BAG1 regulates diverse cellular pathways in chloroplast, such as programmed cell death and stress responses [65]. Therefore, the up-regulation of miR396 and miR_n10 during fruit ripening led us to speculate that these miRNAtarget modules might be involved in the regulation of fruit maturation via depression of chloroplast division and/or maintenance, therefore facilitating coloration in blueberry.
In conclusion, numerous miRNAs and their targets during blueberry fruit maturation were identified by next-generation sequencing approach in the study. The interplays between these miRNAs and their targets depict a picture of the regulatory network involved in blueberry fruit development and maturation, and some miRNAs might play vital roles in fruit maturation and coloration via repressing the expression of auxin-responsive genes or the ones responsible for fruit growth and enlargement ( Figure 8). These findings will facilitate future investigation of miRNAmediated mechanisms underlying blueberry fruit development and maturation.

Plant Materials and RNA Extraction
5-year-old blueberry trees (Vaccinium corymbosum "Northland") from clonal propagation grown at the experimental station in Jilin University (Changchun, China) were used. Fruit samples were randomly collected and mixed from 5 or 6 different plants at 5 developmental stages (green pad, green cup, white, pink and blue fruit), which correspond to the stage S3, S4, S5, S6, S7/8 as described previously [6], respectively. The collected samples were immediately frozen in liquid nitrogen and stored at −80 °C.

Construction and Sequencing of sRNA and Degradome Libraries
Total RNAs were extracted from all the fruit samples using RNAprep Pure Plant Kit (Tiangen Inc., Beijing, China) according to the manufacturer's instructions. Two small RNA (sRNA) libraries and two degradome libraries were constructed from fruits at white stage (WF) and blue stage (BF). Finally, both sRNA and degradome were profiled by single-end sequencing (36 bp) on an Illumina Hiseq 2500. The raw data from the small RNA and degradome sequencing are available at NCBI Short Read Archive (WF-miRNAs, SRX2658700; BF-miRNAs, SRX2676839; WF-degradome, SRX2660622; BF-degradome, SRX2683796).

Plant Materials and RNA Extraction
5-year-old blueberry trees (Vaccinium corymbosum "Northland") from clonal propagation grown at the experimental station in Jilin University (Changchun, China) were used. Fruit samples were randomly collected and mixed from 5 or 6 different plants at 5 developmental stages (green pad, green cup, white, pink and blue fruit), which correspond to the stage S3, S4, S5, S6, S7/8 as described previously [6], respectively. The collected samples were immediately frozen in liquid nitrogen and stored at −80 • C.

Construction and Sequencing of sRNA and Degradome Libraries
Total RNAs were extracted from all the fruit samples using RNAprep Pure Plant Kit (Tiangen Inc., Beijing, China) according to the manufacturer's instructions. Two small RNA (sRNA) libraries and two degradome libraries were constructed from fruits at white stage (WF) and blue stage (BF). Finally, both sRNA and degradome were profiled by single-end sequencing (36 bp) on an Illumina Hiseq 2500. The raw data from the small RNA and degradome sequencing are available at NCBI Short Read Archive (WF-miRNAs, SRX2658700; BF-miRNAs, SRX2676839; WF-degradome, SRX2660622; BF-degradome, SRX2683796).

Analysis of Small RNA Sequencing Data
The raw data were processed to eliminate low quality reads, adapter sequences, junk reads, low complexity sequences and redundant repeats to obtain clean reads Subsequently, the clean reads were subjected to a blast search against Rfam (http://www.sanger.ac.uk/Software/Rfam) and GenBank (http://www.ncbi.nlm.nih.gov/GenBank/) databases to eliminate coding RNAs and Rfam including rRNAs, tRNAs, small nuclear RNAs (snRNAs) and small nucleolar RNAs (snoRNAs). The remaining sequences with length in 18-25 nt were aligned with all the known plant miRNA sequences from miRBase 21.0 (http://www.mirbase.org) to identify candidate miRNAs. Length variation at both 3 and 5 ends and one mismatch inside of the sequence were allowed in the alignment. Furthermore, the candidates were mapped to the reference sequences (cranberry genome, the EST sequences and transcriptome sequences of blueberry) to obtain their precursor sequences using the SOAP2 program [66]. To improve the reliability of miRNA identification, the precursors of the identified miRNAs need to satisfy typical secondary hairpin structure as described previously [28]. Novel miRNAs were identified based on the following criteria: (1) mature miRNAs are present in one arm of the hairpin structures, (2) unpaired residues should be no more than 2 between miRNA and miRNA*, (3) hairpin precursors lack internal loop and bulge, (4) secondary structures of the hairpins are steady with high minimal folding free energy index (MFEI) [28]. The hairpin RNA structures were assessed using the Mfold web server (http://unafold.rna.albany.edu/?q=mfold) with the default parameters.

Phylogenetic Analysis of Known miRNAs and Their Precursors in Blueberry
Sequences of the known miRNAs and their precursors in blueberry and grape (Vitis vinifera) were collected from sequencing data and miRBase 21.0, and then aligned by ClustalW (Version 2.1, http://www.clustal.org). Subsequently, phylogenetic trees were constructed by maximum likelihood method using software MEGA 7 [67] with 1000 replicates of bootstrap value, and visualized by the online Evolview program (http://www.evolgenius.info/evolview/).

Differential Expression Analysis of miRNAs Between WF and BF Libraries
The expression abundance of miRNAs in two libraries was normalized and if the normalized read counts of miRNAs were less than 1 in both two libraries, and then they were removed. The heatmap was generated using the software HemI1.0 [68]. Furthermore, the differential expression of miRNAs between the two libraries was calculated as: Fold change = miRNA normalized reads in BF/miRNA normalized reads in WF. The miRNAs with fold change greater than 2 or less than 0.5, along with p-value ≤ 0.05, were considered as up-regulated or down-regulated ones during maturation, respectively.

Analysis of Degradome Sequencing Data
Degradome analysis and identification of the sliced miRNA targets were performed according to the CleaveLand 3.0 [69]. In brief, the degradome sequences, free from adapter sequences, were mapped to the blueberry transcriptome sequences and ESTs (expressed sequence tags) to generate a degradome density file. Meanwhile, target prediction was performed using TargetFinder program. The alignment between miRNA and mRNA (blueberry transcriptome sequences and ESTs) were generated using 'targetfinder.pl', and the parameters were set as default. Thereafter, the degradome density file and target predictions were compared and significant hits were identified as target plots (t-plots) [26]. The potential target transcripts were classified into five main categories based on the abundance of degradome reads as reported previously [26]. Gene Ontology (GO) annotations were assigned using the Gene Ontology tool (http://www.Geneontology.org/).

qRT-PCR Validation of miRNAs and Their Potential Targets from Degradome Data
Total RNA was isolated from fruit samples at 5 developmental stages (green pad, green cup, white, pink and blue fruit) using RNAPlant Plus Reagent kit (Tiangen Inc., Beijing, China) according to the manufacturer's instruction. Subsequently, polyethylene glycol (PEG) 8000 was used to precipitate of high-molecular-weight (HMW) RNAs. Finally, small RNAs could be easily recovered from supernatant by ethanol precipitation [70]. Reverse transcription was performed with PrimeScript™ RT reagent Kit with gDNA Eraser (Takara Inc., Dalian, China) using specific stem-loop RT primers for miRNAs and the oligo dT primer for targets (Table S4). qRT-PCR was subsequently conducted with an ABI StepOnePlus PCR system and SYBR Premix Ex Taq (Takara Inc.), using ACTIN (CF811156) and U6 as internal control for targets and miRNAs, respectively. Three biological replicates with three technical replicates for each biological replicate were performed for each sample, and data were analyzed by the software ABI StepOnePlus v2.3 (Applied Biosystems, Foster City, CA, USA). The stem-loop primers and target qRT-PCR primers are listed in Tables S4 and S5, respectively.

RNA Ligase-Mediated 5 RACE
RNA ligase-mediated rapid amplification of 5 cDNA ends (RLM-5 RACE) was performed with SMARTer™ RACE cDNA Amplification Kit (Clontech Laboratories Inc., PaloAlto, CA, USA) according to the manufacturer's instructions. Briefly, RNA adapter was ligated to the purified RNAs. The ligation products were reverse transcribed, followed by nested PCR amplifications using universal primers and gene specific primers (Table S6). The RACE products were gel purified, cloned and sequenced.