Abstract

Background

   ‘Fengzao’ is an early-ripening bud mutant of ‘Kyoho’, which matures
   nearly 30 days earlier than ‘Kyoho’. To gain a better understanding of
   the regulatory role of miRNAs in early-ripening of grape berry,
   high-throughput sequencing approach and quantitative RT-PCR validation
   were employed to identify miRNAs at the genome-wide level and profile
   the expression patterns of the miRNAs during berry development in
   ‘Kyho’ and ‘Fengzao’, respectively.

Results

   Nine independent small RNA libraries were constructed and sequenced in
   two varieties from key berry development stages. A total of 108 known
   miRNAs and 61 novel miRNAs were identified. Among that, 159 miRNAs
   identified in ‘Fengzao’ all completely expressed in ‘Kyoho’ and there
   were 10 miRNAs specifically expressed in ‘Kyoho’. The expression
   profiles of known and novel miRNAs were quite similar between two
   varieties. As the major differentially expressed miRNAs, novel_144,
   vvi-miR3626-3p and vvi-miR3626-5p only expressed in ‘Kyoho’,
   vvi-miR399b and vvi-miR399e were down-regulated in ‘Fengzao’, while
   vvi-miR477b-3p up-regulated in ‘Fengzao’. According to the expression
   analysis and previous reports, miR169-NF-Y subunit, miR398-CSD,
   miR3626-RNA helicase, miR399- phosphate transporter and miR477-GRAS
   transcription factor were selected as the candidates for further
   investigations of miRNA regulation role in the early-ripening of grape.
   The qRT-PCR analyses validated the contrasting expression patterns for
   these miRNAs and their target genes.

Conclusions

   The miRNAome of the grape berry development of ‘Kyoho’, and its
   early-ripening bud mutant, ‘Fengzao’ were compared by high-throughput
   sequencing. The expression pattern of several key miRNAs and their
   target genes during grape berry development and ripening stages was
   examined. Our results provide valuable basis towards understanding the
   regulatory mechanisms of early-ripening of grape berry.

Electronic supplementary material

   The online version of this article (10.1186/s12870-018-1516-x) contains
   supplementary material, which is available to authorized users.

   Keywords: Grape, Early ripening, miRNA, ROS

Background

   Grape (Vitis vinifera L.) is an important fruit crop in the world with
   a high nutritional value for fresh berry consuming and wine-making. The
   proper ripening time point is vital either for table or wine grape.
   Fruit ripening is a complex processes involved enormous physiological
   and molecular changes in the color, flavor, quality, and the
   regulations of gene expressions. At present, most of our knowledge
   about the genetic and molecular control of fruit ripening are from the
   model fruit tomato with the various ripening-related mutants [[33]1,
   [34]2].

   Studies have been conducted with the ripening mutants from different
   fruit crops to explore the molecular mechanisms of fruit ripening, such
   as in citrus [[35]3], pear [[36]4], and grape [[37]5]. We have found an
   early-ripening bud mutant of ‘Kyoho’, ‘Fengzao’, which matures earlier
   30 days than ‘Kyoho’ [[38]6]. The phenological, physiological and
   molecular differences between ‘Fengzao’ and ‘Kyoho’ have been
   investigated in the previous studies [[39]6–[40]10]. Histological and
   molecular analysis showed that their genetic background are highly
   uniform [[41]9]. RNA-Seq analysis showed that the main differentially
   expressed genes between ‘Fengzao’ and ‘Kyoho’ are ROS (Reactive Oxygen
   Species) related genes [[42]8]. More proofs need to be explored to
   clarify the early-ripening mutant mechanism of ‘Fengzao’.

   MicroRNAs (miRNAs) are a class of ~ 22 nt endogenous non-coding RNAs
   and negatively regulate the expressions of target genes through mRNA
   degradation or translational inhibition of target mRNAs [[43]11]. Plant
   miRNAs have been extensively analyzed in tomato [[44]1], sweet orange
   [[45]3] and melon [[46]12], etc. miRNAs were confirmed to be involved
   in multiple processes, such as stress response [[47]13], fruit
   development [[48]14, [49]15] and floral bud development [[50]16], etc.

   Recently, some studies have suggested that miRNAs play an important
   role in various fruit development and ripening process including pear
   [[51]17], Lycium barbarum [[52]18], Arabidopsis [[53]19], sweet orange
   [[54]3] and melon [[55]12], etc. Karlova et al. [[56]20] reported that
   CNR and SIAP2a were actively modulated by miR156/157 and miR172 during
   tomato ripening. Gao et al. [[57]1] identified several miRNAs, which
   differentially expressed during fruit development in tomato RIN mutant.
   Saminathan et al. [[58]15] showed that miR156, miR156a, miR159a,
   miR159b, and miR319b were upregulated during the later stages of
   pomegranate fruit development. There are a relatively large number of
   grape miRNAs deposited in miRBase Release21 and their expression
   profiles have been analyzed in developing berries [[59]21–[60]23]. Some
   reports suggested that the mutations of the mutant have modified the
   regulation of miRNAs expressions between the mutant and wild genotypes
   [[61]24–[62]26]. However, the mechanism underlying the miRNA-mediated
   regulation of early-ripening of grape berry is largely unknown.
   Therefore, identification and characterization of miRNAs in grape early
   ripening bud mutant may provide new insights for grape berry ripening
   process.

   To gain a better understanding of the regulatory role of miRNAs in
   early-ripening of grape berry, the miRNA profiles during the berry
   development of ‘Kyoho’ and its early-ripening bud mutant, ‘Fengzao’
   were investigated using an Illumina HiSeq 2000 platform. The
   differentially expressed miRNAs during the berry development of
   ‘Fengzao’ and ‘Kyoho’ were identified, and the corresponding target
   genes were predicted. In addition, the expression patterns of several
   miRNAs and their target genes were further examined by qRT-PCR
   analyses. The results provided some insights into the regulatory roles
   of miRNAs during grape berry development, and the data presented here
   will lay a foundation for future studies of grape early-ripening.

Results

Small RNA profiles of grape berry development

   Nine independent small RNA libraries were generated from key berry
   development stages of ‘Fengzao’ and ‘Kyoho’ (Table [63]1). The sampling
   time points were selected based on our earlier researches which
   indicated that these berry development stages are key to reveal the
   significant transcriptional and metabolic varitions between ‘Fengzao’
   and ‘Kyoho’ [[64]7, [65]8]. For samples of ‘Fengzao’, 11,524,331 reads
   were obtained on average per sample (Additional file [66]1: Table S1).
   While for ‘Kyoho’, the averages are 12,717,569 reads
   (Additional file [67]1: Table S1). After filtering the low-quality
   reads, adapter contaminates, reads with ploy Ns, clean reads ranging
   from 10,245,269 to 12,987,903 were obtained per sample. Then, reads of
   length < 18 were excluded, resulting in 9,103,202 reads on average.
   Subsequntly, more than 66% of the screened reads were successfully
   mapped to grape reference sequence.

Table 1.

   The sampling time and corresponding development stages of grape berries
   in this study
   Development stage Sampling date Sampling date Stage code
   Characteristics
   Fengzao Code Kyoho Code
   E-L33 6.17 FZ1 6.27 KY1 I hard green berries
   7.4 KY2 II
   E-L34 6.27 FZ2 7.16 KY3 III starting to soften
   E-L35 7.4 FZ3 7.23 KY4 IV véraison
   E-L37 7.11 FZ4 8.3 KY5 V Berries not quite ripe
   [68]Open in a new tab

   Samples were collected in 2016

   The length of the small RNA (sRNA) ranged from 18 to 35 nt were counted
   (Fig. [69]1). The dominant and abundant sRNA length were 21 and 24 nt
   regardless of the development stages and whatever in ‘Kyoho’ or
   ‘Fengzao’ (Fig. [70]1). The 21 nt sRNAs were the most abundant,
   followed by 24 nt sRNAs. The amount of 21 nt and 24 nt accounted for
   60% of all the sRNA reads. However, the number of 24 nt sRNAs is a
   little more than that of 21 nt at the sample of FZ2 and KY2.

Fig. 1.

   [71]Fig. 1
   [72]Open in a new tab

   Length distribution of unique sRNAs in grape berry libraries. Y-axis
   represents the percentages of small RNA identified in each sample.
   X-axis represents the length of small RNA. FZ1–4 and KY1–5 represent
   the samples from different berry development stages as indicated in
   Table [73]1. The sample codes are the same as below

   The size distribution of 20, 22 and 23-nt redundant reads has no
   obvious difference between ‘Kyoho’ and ‘Fengzao’. From Fig. [74]1, it
   could be observed that the length distribution profiles of sRNAs were
   very similar between ‘Kyoho’ and ‘Fengzao’ during the whole berry
   development stages. The common and specific sRNAs were compared between
   any two libraries, showing that more than 60% of the total sRNAs were
   shared.

   The sRNAs were analyzed using the blastn and blastall alignments
   against Genbank and Rfam databases to filter out the known non-coding
   RNAs: snoRNAs, snRNAs, tRNAs, and rRNAs. All of nine libraries showed
   similar compositions of non-conding RNA types (Table [75]2). The number
   of reads differed in each category. Particularly, the proportion of
   repeat and NAT (natural antisense short interfering RNA) specific reads
   was high (10.6 and 10.7%) among the unique reads. Most abundant
   non-coding RNA class was rRNA for unique reads with the mean value of
   8.21%, and other RNAs had quantity with the average of 0.16% (snRNA)
   and 0.58% (snoRNA) (Table [76]2). The remaining sequences were examined
   to identify miRNAs, and the un-annotated sequences were used for novel
   miRNAs identification. Afterwards, the known miRNA is 4.65%, novel
   miRNA is 2.39% on average.

Table 2.

   Summary of small RNA sequencing data and annotation after alignment to
   GenBank and Rfam in nine libraries
   FZ1 FZ2 FZ3 FZ4 KY1 KY2 KY3 KY4 KY5
   Total 5,511,599 6,292,055 6,206,300 6,835,422 7,080,977 6,531,036
   6,498,722 6,944,941 6,876,953
   known_miRNA 293,764
   (5.33%) 256,923
   (4.08%) 195,676
   (3.15%) 240,682
   (3.52%) 304,883
   (4.31%) 220,877
   (3.38%) 288,939
   (4.45%) 468,847
   (6.75%) 473,968
   (6.89%)
   rRNA 401,277
   (7.28%) 461,547
   (7.34%) 819,911
   (13.21%) 555,173
   (8.12%) 490,108
   (6.92%) 408,859
   (6.26%) 614,342
   (9.45%) 530,932
   (7.64%) 526,145
   (7.65%)
   tRNA 2
   (0.00%) 1
   (0.00%) 0
   (0.00%) 1
   (0.00%) 1
   (0.00%) 2
   (0.00%) 1
   (0.00%) 0
   (0.00%) 0
   (0.00%)
   snRNA 6216
   (0.11%) 7435
   (0.12%) 8612
   (0.14%) 7961
   (0.12%) 8801
   (0.12%) 6698
   (0.10%) 6213
   (0.10%) 7397
   (0.11%) 7916
   (0.12%)
   snoRNA 29,043
   (0.53%) 45,448
   (0.72%) 25,971
   (0.42%) 29,690
   (0.43%) 55,917
   (0.79%) 47,378
   (0.73%) 52,260
   (0.80%) 31,061
   (0.45%) 28,381
   (0.41%)
   repeat 705,647
   (12.80%) 696,292
   (11.07%) 510,942
   (8.23%) 605,622
   (8.86%) 759,921
   (10.73%) 750,842
   (11.50%) 743,158
   (11.44%) 715,266
   (10.30%) 723,924
   (10.53%)
   NAT 609,654
   (11.06%) 640,231
   (10.18%) 620,372
   (10.00%) 636,276
   (9.31%) 664,891
   (9.39%) 647,042
   (9.91%) 688,849
   (10.60%) 866,898
   (12.48%) 935,099
   (13.60%)
   novel_miRNA 123,239
   (2.24%) 150,151
   (2.39%) 124,871
   (2.01%) 181,311
   (2.65%) 171,565
   (2.42%) 127,118
   (1.95%) 168,580
   (2.59%) 166,758
   (2.40%) 200,238
   (2.91%)
   TAS 16,101
   (0.29%) 13,801
   (0.22%) 11,967
   (0.19%) 14,884
   (0.22%) 14,451
   (0.20%) 13,756
   (0.21%) 23,823
   (0.37%) 17,509
   (0.25%) 21,585
   (0.31%)
   Exon_sense 124,655
   (2.26%) 112,184
   (1.78%) 96,842
   (1.56%) 110,968
   (1.62%) 125,125
   (1.77%) 101,251
   (1.55%) 154,621
   (2.38%) 152,438
   (2.19%) 160,113
   (2.33%)
   Exon_antisense 113,779
   (2.06%) 84,968
   (1.35%) 65,792
   (1.06%) 79,570
   (1.16%) 101,249
   (1.43%) 80,722
   (1.24%) 146,666
   (2.26%) 124,805
   (1.80%) 125,799
   (1.83%)
   Intron_sense 50,417
   (0.91%) 67,408
   (1.07%) 56,941
   (0.92%) 72,555
   (1.06%) 59,421
   (0.84%) 67,508
   (1.03%) 70,087
   (1.08%) 76,390
   (1.10%) 85,546
   (1.24%)
   Intron_antisense 29,843
   (0.54%) 30,098
   (0.48%) 30,491
   (0.49%) 30,652
   (0.45%) 32,891
   (0.46%) 34,738
   (0.53%) 40,122
   (0.62%) 47,752
   (0.69%) 57,420
   (0.83%)
   Un_annotated 3,007,962
   (54.58%) 3,725,568
   (59.21%) 3,637,912
   (58.62%) 4,270,077
   (62.47%) 4,291,753
   (60.61%) 4,024,245
   (61.62%) 3,501,061
   (53.87%) 3,738,888
   (53.84%) 3,530,819
   (51.34%)
   [77]Open in a new tab

   rRNA/snRNA/snoRNA/tRNA considered; TAS trans-acting small interfering
   RNA, NAT natural antisense short interfering RNA

Identification of known and novel miRNAs in grape berry

   After further sequence analysis, a total of 159 miRNAs were identified
   in ‘Fengzao’, among which, known miRNAs were 100. While for ‘Kyoho’,
   total and known miRNAs were 169 and 108, respectively
   (Additional file [78]2: Table S2). The sequence of all known and novel
   miRNAs were shown in Additional file [79]3: Table S3. For ‘Fengzao’,
   134 sRNAs were common among four development stages (Fig. [80]2). While
   for ‘Kyoho’, it was 141. Among all these miRNAs, 159 were shared
   between ‘Kyoho’ and ‘Fengzao’, 10 was specific for ‘Kyoho’: novel_139,
   novel_144, vvi-miR156a, vvi-miR156e, vvi-miR171h, vvi-miR3626-3p,
   vvi-miR3626-5p, vvi-miR845a, vvi-miR845c, vvi-miR845d. The expression
   levels of most specific miRNAs in ‘Kyoso’ were low except
   vvi-miR3626-3p and vvi-miR3626-5p. The target genes predicted by
   psRobot softwar were shown in Additional file [81]4: Table S4. Most of
   the target genes were related to activities of ATPase, Zinc finger,
   Leucine-rich repeat, NADP-dependent oxidoreductase domain, Protein
   kinase, expecially Serine/threonine-protein kinase, ATP binding sites,
   Phosphorylation, etc.

Fig. 2.

   [82]Fig. 2
   [83]Open in a new tab

   Venn diagram showing miRNAs identified in this study. The left and
   right showed the miRNAs profile between different berry developmental
   stages of ‘Fengzao’ and ‘Kyoho’, respectively. The middle showed the
   miRNAs profile shared between ‘Fengzao’ and ‘Kyoho’

   In this study, these known miRNAs could be classified into 35 conserved
   miRNA families (Additional file [84]5: Figure S1;
   Additional file [85]6: Table S5). Most of the conserved families were
   present across the different plants species (Additional file [86]5:
   Figure S1). The miRNAs composition in each family identified in this
   study were similar to that of Malus domestica (Additional file [87]5:
   Figure S1), indicating high conservation of miRNA families and high
   homology between two species. Some known but less-conserved miRNAs were
   also found in this study (MIR845_2, MIR3630, MIR3631, MIR845_3,
   MIR2950). Among these, MIR169_2 were the largest represented families
   with 11 members, followed by MIR399 and MIR171_1, with six and five
   members, respectively. Of the remaining families, most of them
   comprised two to four members. Furthermore, various family members from
   the same miRNA family showed diverse expression levels. For example,
   the number of MIR396 family members reads ranged from1839–29,792, while
   for MIR169_2 family members, reads ranged from 2 to 987. For
   vvi-miR399d (2–37 reads) and vvi-miR399a (73–645 rerads), both of them
   belong to the same famlily MIR399, but the expression levels showed
   largely variable.

   Nucleotide bias analysis of these known miRNAs revealed that uridine
   (U) is the most common preferred nucleotide at the first position of
   the 5′end (> 85%); while adenine (A; ~ 40%) was the most common
   nucleotide at both the 10th and 11th nucleotide matched to the cleavage
   site of the targets (Additional file [88]7: Figure S2A). For 18- to
   28-nt miRNAs, the base bias at the first position from the 5′ end had a
   strong preference for U. Nucleotides A and G predominately occupied the
   first position base bias for the 29- and 30-nt miRNAs, respectively
   (Additional file [89]7: Figure S2B). Nucleotide U was preferred more
   than 80% of the time as a first base for 18- to 30-nt miRNAs.

   After removing conserved and known miRNAs, the remaining sequences were
   used to predict novel miRNAs according to the sequence similarity and
   the formation of a stable stem-loop structure in the precursor.
   Finally, 59 and 61 potential novel miRNAs were predicted from ‘Fengzao’
   and ‘Kyoho’, respectively (Additional file [90]2: Table S2). Most of
   predicted novel miRNAs expressed both in ‘Fengzao’ and ‘Kyoho’, except
   novel_139 and novel_144, which only expressed in ‘Kyoho’.

Prediction of putative target genes for known and novel miRNAs

   A total of 2162 unigenes were predicted to be targets of 105 miRNAs
   using PsRobot software (Additional file [91]4: Table S4). A relatively
   high proportion of target genes were annotated as transcription factors
   (TFs) and protein-like domains. A number of miRNAs had multiple
   putative target genes, the numbers of targets for each ranging from 1
   (vvi-miR168, vvi-miR3627-3p) to 331 (for vvi-miR3629a-5p), only 24
   miRNAs below 4 targets, suggesting that these miRNAs may have diverse
   biological functions. In contrast, several putative target genes were
   targeted by multiple miRNAs with up to 14 miRNAs for a hypothetical
   protein (VIT_11s0103g00390), on average 1.44. Noticeably, only 1 out of
   the putative 61 novel miRNAs (novel_141) was successfully predicted to
   target 192 unigenes.

   To reflect a global overview of the regulatory functions of miRNAs, the
   GO terms of all target genes were analyzed through a GO annotation.
   Among the 2162 target genes, 1711 target genes had GO terms. As shown
   in Additional file [92]8: Figure S3, target prediction analysis showed
   that the identified targets regulated in a wide spectrum of biological
   processes, cellular components, and molecular functions. The target
   genes with molecular functions in nucleic acid binding presented the
   highest percentage, corresponding to 12.36%, followed by ATP binding,
   10.7% (Additional file [93]8: Figure S3). GO enrichment analysis
   revealed that the top three terms were copper ion binding;
   oxidoreductase activity, oxidizing metal ions; oxidoreductase activity,
   acting on diphenols and related substances as donors (Fig. [94]3).
   Based on a KEGG pathway analysis, the target genes in this study were
   involved in 115 different pathways. The most top 5 enriched pathways
   include protein processing in endoplasmic reticulum, RNA transport,
   spliceosome, ubiquitin mediated proteolysis (Additional file [95]9:
   Figure S4). There were no enrichment results for KEGG analysis.

Fig. 3.

   [96]Fig. 3
   [97]Open in a new tab

   Scattergram of overrepresented GO terms (P < 0.05) in molecular
   function categories from GO enrichment analysis of miRNA targets gene
   of all the miRNAs identified in this study using ClusterProfiler.
   Enrichment term is represented by colored dots (red indicates high
   enrichment and blue indicates low enrichment)

Characterization of expressed miRNAs in berries

   The relative expression levels of miRNAs were estimated as transcripts
   per million (TPM). The TPM values drastically varied among 9 miRNA
   samples. Some miRNAs were highly expressed in both ‘Fengzao’ and
   ‘Kyoho’, there are 36 sRNAs which accumulated at more than 1000 TPM
   (Additional file [98]10: Table S6).

   In order to define the relationships of the 9 different samples from
   different berry stages of ‘Fengzao’ and ‘Kyoho’, a correlation matrix
   using Pearson coefficients (Fig. [99]4) were established to evaluate
   the level of similarity among the samples based on the miRNA expression
   levels. The high correlation was observed with average R^2 ≥ 0.896
   among ‘Fengzao’ samples (Fig. [100]4). Most of the correlation
   coefficients of ‘Kyoho’ samples was above 0.802 except the coefficients
   between KY3 and KY5 (Fig. [101]4).The correlation results showed that
   the berryies from different development periods of ‘Fengzao’ and
   ‘Kyoho’ are moderately related to each other, highlighting the
   similarities among these samples, except the KY3 samples.

Fig. 4.

   Fig. 4
   [102]Open in a new tab

   Correlation matrix of the 9 grape samples used to define the miRNA
   atlas. One minus Pearson correlation was used as a metric distance

   To better understand the overall expression profile of the common
   miRNAs in 9 samples, K-means clustering method were conducted further
   to elucidate the overall expression trends of the miRNAs based on the
   value of Log10(TPM + 1) –transformed. As a result, cluster analysis of
   the expression patterns divided these miRNAs into 16 groups
   (Fig. [103]5). The miRNAs and corresponding subclusters for the 16
   groups are shown in Additional file [104]8: Table S8 and Fig. [105]5.
   As shown in Fig. [106]5, most subclusters showed similar expression
   patterns between ‘Kyoho’ and ‘Fengzao’. Some differences in miRNA
   expression patterns were revealed between ‘Kyoho’ and ‘Fengzao’ in
   subcluster_1, subcluster_9, subcluster_10 (Fig. [107]5).

Fig. 5.

   [108]Fig. 5
   [109]Open in a new tab

   Cluster anlyses of the miRNAs expression patterns across various
   developmental stages in ‘Fengzao’ and ‘Kyoho’. Clustering was performed
   using k-means statistics. X-axis represents different sampling time
   points and Y-axis represents the expression level of miRNAs (log
   (TPM + 1)). The blue line was the fitted curve to match the expression
   change trend along the berry development

   The subcluster_1 comprised miR3626-5p, miR3626-3p, novel_144,
   novel_169, and novel_139; the expression levels of these miRNAs were
   very low in ‘Fengzao’, but gradually increased in ‘Kyoho’during berry
   developmental stages. These miRNAs were annotated as Serine/threonine-
   protein kinase, Sodium/calcium exchanger membrane region,
   phosphopantetheine-binding domain and SEP domain. The subcluster_9
   composed of novel_117, novel_128 and miR3624-5p (WD40 repeat,
   FAD-binding), they gradually decreased in ‘Fengzao’along the berry
   development, but fluctuated in ‘Kyoho’. The subcluster_10 included
   miR398a (Thioredoxin-like fold), miR169b (Actin-related protein),
   novel_182, novel_64 and novel_46, they gradually increased in
   ‘Fengzao’, while they firstly decreased in KY2 of ‘Kyoho’, then slowly
   increased.

Identification and different expression of miRNAs

   The normalized expression levels of miRNAs were compared between the
   samples from ‘Kyoho’ and ‘Fengzao’ to identify differentially expressed
   miRNAs based on the criteria (|fold-change| > =2, q-value < 0.05). The
   comparisons were made between any two adjacent sampling time points
   within ‘Fengzao’ or ‘Kyoho’, respectivley.

   There are few common differentially expressed miRNAs obtained from
   pairwise comparisons of the samples. Considering just for the samples
   of ‘Fengzao’, only one miRNA, miR159a, was differentially expressed
   between any two adjacent sampling time points, i.e., were
   differentially expressed in all stages of ‘Fengzao’. While for ‘Kyoho’,
   there are six miRNAs, miR159a, miR164a, miR2111-3p, miR3625-5p, miR390
   and miR828a. These results indicated that a larger number of miRNAs
   were differentially expressed during all stages of berry development.

   The overall comparison between ‘Kyoho’ and ‘Fengzao’ revealed that
   novel_144, vvi-miR3626-3p, vvi-miR3626-5p, vvi-miR399b, vvi-miR399e,
   vvi-miR477b-3p were differentially expressed (Table [110]3). The
   information about the expression and target genes of the differently
   expressed miRNAs are listed in Table [111]3. Some miRNAs showed
   genotype-specific expression patterns. For example, several miRNAs were
   expressed only in ‘Kyoho’, but completely didn’t express in any berry
   development stages of ‘Fengzao’, such as novel_144, vvi-miR3626-3p and
   vvi-miR3626-5p (Table [112]3), indicating that these miRNAs may have a
   specific role in berry development of grape. Some miRNAs were
   down-regulated in ‘Fengzao’ when compared with ‘Kyoho’; these miRNAs
   are miR399b and miR399e. The target genes of these miRNAs are annotated
   as NADH-ubiquinone oxidoreductase, Acyl carrier protein (ACP), Acyl
   carrier protein (ACP), Serine/threonine- protein kinase, SEP domain,
   inner membrane protein, Sodium/calciumoi exchanger membrane region,
   Phosphate permease, Ubiquitin-related domain, etc. Contrary to this,
   miR477b-3p was upregulated in ‘Fengzao’. The annotation of miR477b-3p
   target genes are Ubiquitin-related domain, Tetratricopeptide repeat,
   methyltransferase, LRR receptor-like serine threonine-protein kinase
   gso1-like, etc.

Table 3.

   Significantly differentially expressed miRNAs identified in ‘Fengzao’
   compared with ‘Kyoho’
   miRNA FZ KY log2FoldChange pval padj Target genes Annotation from GO
   and interPro
   novel_144 0 43.34 −4.0075 2.02E-14 1.12E-12
   vvi-miR3626-3p 0 29.95 −3.582 2.15E-11 7.94E-10 VIT_06s0061g00770
   VIT_11s0016g05110
   VIT_01s0244g00030
   VIT_18s0076g00210
   VIT_03s0088g00450
   VIT_10s0116g00210
   VIT_16s0050g01870
   VIT_14s0083g01170
   VIT_16s0050g01490
   VIT_02s0087g00400
   VIT_18s0001g06780
   VIT_10s0116g00290
   VIT_19s0014g00220
   VIT_12s0028g00300
   VIT_04s0023g00150
   VIT_14s0066g01000
   VIT_14s0171g00260 50s ribosomal protein l9; Ribosomal protein L9/RNase
   H1, N-terminal
   protein brittle- chloroplastic amyloplastic-like; Mitochondrial carrier
   domai
   Acyl carrier protein (ACP); Polyketide synthase,
   phosphopantetheine-binding domain
   embryogenesis-associated protein emb8; Alpha/Beta hydrolase fold
   fkbp12-rapamycin complex-associated protein; Phosphatidylinositol
   3−/4-kinase, PIK-related kinase
   hypothetical protein; UBX||Ubiquitin-related domain||SEP domain
   CAAX amino terminal protease family protein; CAAX amino terminal
   protease
   mitochondrial carrier protein; Mitochondrial substrate/solute carrier
   NADH-ubiquinone oxidoreductase 13 kda-b subunit; ETC complex I subuni
   protein notum homolog; Alpha/Beta hydrolase fold||Pectinacetylesterase
   probable serine threonine-protein kinase at1g54610-like; Protein
   kinase, ATP binding site||Serine/threonine-
   plant UBX domain-containing protein 3; SEP domain
   big map kinase; Serine/threonine- protein kinase ||Mitogen-activated
   protein (MAP) kinase, ATP binding site
   glutamyl-tRNA amidotransferase subunit a-like; Tetratricopeptide
   repeat-containing domain
   inner membrane protein; Tetratricopeptide repeat-containing domain,
   embrane insertase OXA1/ALB3/YidC
   phosphate phosphoenolpyruvate translocator; Drug/metabolite transporter
   phosducin-like protein 3; Phosducin, thioredoxin-like domain
   vvi-miR3626-5p 0 82.19 −4.7967 1.11E-21 1.23E-19 VIT_10s0116g00910
   VIT_10s0003g03500
   VIT_16s0050g02710
   VIT_09s0018g01840
   VIT_08s0007g04400
   VIT_00s0179g00190 receptor protein; Serine/threonine-protein kinase,
   Protein kinase, ATP binding site
   katanin p60 ATPase-containing subunit a-like 2-like; AAA+ ATPase
   domain||ATPase, AAA-type, core
   probable receptor-like protein kinase at1g67000-like; Wall-associated
   receptor kinase galacturonan-binding domain
   vacuolar cation proton exchanger 2; Sodium/calcium exchanger membrane
   region
   dead-box atp-dependent RNA helicase 38-like; P-loop containing
   nucleoside triphosphate hydrolase||RNA helicase
   transcription factor jumonji domain-containing protein; JmjC domain
   vvi-miR399b 38.95 159.83 −1.5213 0.001754 0.032445 VIT_08s0007g02840
   VIT_06s0004g04430
   VIT_13s0067g03280 LETM1-like protein; LETM1-like
   ubiquitin carrier protein; Ubiquitin-conjugating enzyme, active site
   high affinity inorganic phosphate transporter; General substrate
   transporter||Phosphate permease
   vvi-miR399e 18.72 76.74 −1.761 5.11E-06 0.000142 VIT_13s0067g03280
   VIT_08s0007g02840
   VIT_13s0067g03280
   VIT_06s0009g02380 high affinity inorganic phosphate transporter;
   General substrate transporter||Phosphate permease
   LETM1-like protein; LETM1-like
   high affinity inorganic phosphate transporter; General substrate
   transporter||Phosphate permease
   pseudouridine-5 -phosphate glycosidase; Pseudouridine-5′-phosphate
   glycosidase
   vvi-miR477b-3p 780.94 283.54 1.2555 0.001439 0.031955 VIT_18s0001g01600
   VIT_13s0084g00080
   VIT_09s0002g03460
   VIT_06s0004g08220
   VIT_18s0001g03710
   VIT_12s0055g00360
   VIT_18s0041g01430
   VIT_09s0002g02670
   VIT_18s0122g00260
   VIT_13s0158g00050
   VIT_07s0104g00810
   VIT_01s0010g03820
   VIT_01s0010g00380
   VIT_18s0001g05870 ATP binding; Putative
   S-adenosyl-L-methionine-dependent methyltransferase
   65-kda microtubule-associated protein 5-like protein regulator;
   Microtubule-associated protein, MAP65
   clavaminate synthase-like protein; Taurine catabolism dioxygenase
   TauD/TfdA
   26 s protease regulatory subunit 6b homolog; P-loop containing
   nucleoside triphosphate hydrolase ||ATPase
   tmv resistance protein n-like; Leucine-rich repeat|P-loop containing
   nucleoside triphosphate hydrolase
   protein-protein interaction regulator family protein; Pinin/SDK/MemA
   protein
   tmv resistance protein n-like; Leucine-rich repeat, typical subtype
   RNA methyltransferase family protein; RNA methyltransferase TrmA
   superkiller protein 3-like protein; Tetratricopeptide repeat
   hypothetical protein; Peptidase S10, serine
   carboxypeptidase||Alpha/Beta hydrolase fold
   biotin lipoate a b protein ligase family protein; Biotin/lipoate A/B
   protein ligase||Octanoyltransferase
   LRR receptor-like serine threonine-protein kinase gso1-like; Protein
   kinase, ATP binding site
   LRR receptor-like serine threonine-protein kinase gso1-like;
   Leucine-rich repeat, Protein kinase, ATP binding site
   pentatricopeptide repeat-containing protein at2g13600-like;
   Tetratricopeptide-like helical
   [113]Open in a new tab

qRT-PCR assay

   To validate the transcriptome data, six differentially expressed and
   four specific miRNAs were selected for real-time quantitative reverse
   transcription (qRT)-PCR analysis. Most tested miRNAs showed similar
   expression trends in qRT-PCR data as sRNA-seq revealed (Fig. [114]6).
   The expression of novel_144, vvi-miR3626-3p, vvi-miR3626-5p miRNA
   indeed couldn’t be detected in ‘Fengzao’; accordingly, the expression
   of the target gene of vvi-miR3626-3p and vvi-miR3626-5p,
   VIT_08s0007g04400, expressed more in ‘Fengzao’ than in ‘Kyoho’. The
   expression of the target gene of miR399b and miR399e,
   VIT_13s0067g03280, also expressed more in ‘Fengzao’ than in ‘Kyoho’ due
   to the down-regulated expression of the miRNAs in ‘Fengzao’. The
   expression of the target genes of miR398 (VIT_01s0127g00520), miR169
   (VIT_01s0011g05560) and miR 159 (VIT_06s0009g02480) expressed more in
   ‘Fengzao’ than in ‘Kyoho’. This result illustrated that our
   high-throughput data were reliable.

Fig. 6.

   [115]Fig. 6
   [116]Open in a new tab

   Stem-loop RT-qPCR validation of significantly differentially expressed
   miRNAs (vvi-) and Quantitative RT-PCR of their target genes (VIT_) at
   different berry development stages in ‘Fengzao’ and ‘Kyoho’. Relative
   quantity is based on the expression of the reference gene U6 and
   ubiquitin1 gene for miRNAs and target genes, respectively. X-axis
   indicates different stages (as indicated in Table [117]1) and Y-axis
   the expression of miRNA or target genes relative to reference gene.
   Data are mean ± SD from three biological replicates. *, P < 0.05; **,
   P < 0.01 by Student t test

Discussion

Characterization of miRNAs during berry development

   A lot of miRNAs have been identified and their expression atlas were
   established in grape [[118]21–[119]23, [120]27]. However, the miRNAs
   involved in grapes berry ripening process remained largely unknown. To
   characterize the miRNAs during grape berry development, nine libraries
   were constructed to identify miRNAs and profile their expressions
   across the key developmental stages of ‘Kyoho’ and ‘Fengzao’,
   respectively. As shown in Fig. [121]1, 21- to 24-nt sRNAs were dominant
   and occupied more than 80% of the total sRNAs. Especially, 21-nt sRNAs
   were more abundant than 24-nt sRNAs (Fig. [122]1), which was consistent
   with the previous studies in grape [[123]22, [124]23, [125]27]. It
   showed that the expression levels of miRNAs varied at different
   developmental stages, indicating miRNAs play diverse roles during grape
   berry development.

   Sequence analysis of sRNAs resulted in characterization of 108 known
   miRNAs and 61 novel miRNAs (Additional file [126]2: Table S2). The
   known miRNAs belonged to 35 conserved miRNA families and most of them
   were highly conserved among diverse plant species
   (Additional file [127]5, Figure S1). Normally, conserved plant miRNAs
   regulate homologous target genes at identical target sites and that
   these miRNAs may play the similar function in different species
   [[128]28]. Target gene analysis showed that a single miRNA
   simultaneously has multiple target genes (Additional file [129]4: Table
   S4), such as miR156, miR159, miR172, and miR164, which is consistent
   with the previous reports [[130]22]. Meanwhile, some miRNAs from
   different families could have the same target genes. The similar cases
   were documented in Lycium barbarum [[131]18], tomato [[132]29] and
   melon [[133]12].

   This study revealed many commonly expressed miRNAs as the previous
   researches [[134]22, [135]23] despite the employments of very different
   varieties and berry stages, especially for the most highly expressed
   miRNAs such as miR156, miR159, miR172 and miR164. The amount of
   identified known and novel miRNAs were a little more than that of Wu et
   al. [[136]3]. They analyzed small RNA expression profiles in navel
   orange of ‘Fengjie 72–1’ and its spontaneous late-ripening mutant
   ‘Fengwan’ at 170 days after flowering (DAF). They identified 107
   conserved miRNAs and 21 novel miRNAs. It may be due to more fruit
   developmental stages covered in this study.

   In contrast to conserved miRNAs, novel miRNAs are usually
   species-specific, which expressed at low levels and appeared to lack
   targets based on the current criteria or the research limits. In this
   study, only one novel miRNA (novel_141) was successfully predicted to
   target protein- coding genes, suggesting that most novel miRNAs have as
   yet undetermined roles at present.

Specific expressed miRNAs in ‘Kyoho’

   The overall comparison of miRNAs expressions showed that all the miRNAs
   expressed in ‘Fengzao’ also expressed in ‘Kyoho’, but 10 miRNAs
   (novel_139, novel_144, vvi-miR156a, vvi-miR156e, vvi-miR171h,
   vvi-miR3626-3p, vvi-miR3626-5p, vvi-miR845a, vvi-miR845c, vvi-miR845d)
   only presented in ‘Kyoho’. It suggested that the mutation in ‘Fengzao’
   had resulted in the alteration of the corresponding miRNA expressions.
   MicroRNA regulates gene expression through inhibiting translation or
   degrading mRNA at a certain site, i.e., it down-regulates the
   expression of its target gene [[137]30]. This implies that expression
   of miRNAs and their targets were negatively correlated. So, the miRNAs
   only expressed in ‘Kyoho’ means that the target of these miRNAs in
   ‘Fengzao’ were accordingly not suppressed and up-regulated. Gao et al.
   [[138]1] demonstrated the up-regulation of the target genes during the
   ripening of wild-type tomato fruit, while their concurrent
   down-regulation in never-ripe rin mutant. Palumbo et al. [[139]31] also
   observed that some miRNAs are deactivated and the switch genes are
   expressed during the transition from immature to mature growth of grape
   berry. Interestingly, vvi-miR156e, vvi-miR3626-5p, vvi-miR845a and
   vvi-miR845d were as parts of some miRNAs identified by Palumbo et al.
   [[140]31] as miRNA/switch gene pairs for the transition from
   immature-to-mature stage in grape, only expressed in ‘Kyoho’ in this
   study. Pilati et al. [[141]32] also showed that 13 out of 80 candidate
   “ABA-responsive switch genes” were predicted to be regulated
   post-transcriptionally by miRNAs including vi-miR156e, vvi-miR3626-5p
   and vvi-miR845a. These evidences suggested the important roles of these
   miRNAs for the switching of grape berry development. It need to be
   further explored in the future.

   Kullan et al. [[142]22] showed that miR156 and miR164 involving in
   grape berry ripening. miR156 decreased from the vegetative to the
   reproductive phase of Arabidopsis and its targets, SQUAMOSA PROMOTER
   BINDING PROTEIN-LIKE (SPLs), are key participators in regulating plant
   phase transitions [[143]33]. The target gene of miR164, Vv-NAC, has
   been demonstrated a special function for grape berry ripening
   [[144]34]. Several miR156-targeted genes were differentially expressed
   in pre- and post-anthesis ovaries of tomato, indicating its regulation
   role for fruit development [[145]35]. MiR156/157 could affect the
   tomato ripening process by modulating the known ripening regulators CNR
   and SIAP2a [[146]36]. vvi-miR 156a and vvi-miR156e only expressed in
   ‘Kyoho’ as confirmed by qRT-PCR (data not shown), which suggested that
   they may be related to the early-ripening mutation of ‘Fengzao’ due to
   its absence in ‘Fengzao’.

   Specific set of miRNAs of different Vitis varieties and species have
   been revealed [[147]22]. For example, Kullan et al. [[148]22] and Mica
   et al. [[149]21] previously didn’t detect the expression of vvi-miR845
   in grape. But the expression of vvi-miR845a, vvi-miR845c, vvi-miR845d
   were detected in ‘Kyoho’ in this study. Palumbo et al. [[150]31] showed
   that vvi-miR845a and vvi-miR156e had the same target gene, Pyruvate
   dehydrogenase kinase (VIT_214s0060g00420). Downregulation of the
   pyruvate dehydrogenase gene PDHE1α by VIGS had inhibited respiration
   and ATP biosynthesis, whilst promoted the accumulation of sugar, ABA
   and ETH, which then accelerated the ripening of strawberry [[151]37].
   And both of them were absent in ‘Fengzao’, indicated that vvi-miR845
   and vvi-miR156 may have the similar roles for the grape berry ripening
   as in strawberry. In this study, vvi-miR845a, vvi-miR845c and
   vvi-miR845d were all predicted to target Gag-polypeptide of LTR
   copia-type (VIT_209s0002g08320) retrotransposon. Previous studies have
   suggested retrotransposon was closely related to the formation of grape
   bud mutant [[152]38]. It need to be further explored whether the
   early-ripening mutant of ‘Fengzao’ is related to retrotransposon or
   not.

Expression patterns of miRNAs between ‘Fengzao’ and ‘Kyoho’

   Expression profiles based on the k-means method showed most of the
   miRNAs had the similar expression patterns between ‘Fengzao’ and
   ‘Kyoho’ during fruit development (Fig. [153]5), which indicated that
   the mutation in ‘Fengzao’ didn’t alter the expressions of miRNAs very
   much compared to ‘Kyoho’. The results also suggested that miRNAs from
   the same cluster have parallel expression patterns (Fig. [154]5).
   Furthermore, different clusters contained some similar miRNAs (from the
   same miRNA family) (Fig. [155]5). This implied that those miRNAs even
   from the same family had different expression patterns and they
   synergistically regulate the berry development. A similar situation has
   been also reported [[156]17].

   The differential expression profiles of some miRNAs were identified at
   different berry developmental stages between ‘Fengzao’ and ‘Kyoho’
   during berry development (Fig. [157]5). For example, 5 miRNAs in
   subcluster 1 of Fig. [158]5, including miR3626-5p, miR3626-3p,
   novel_144, novel_169, and novel_139, their expression pattern were
   entirely different between ‘Kyoho’ and ‘Fengzao’. They expressed very
   lowly in ‘Fengzao’ and gradually increased in ‘Kyoho’ along the berry
   development. The expression profiles of miRNAs in subcluster 10 of
   Fig. [159]5, miR398a, miR169b, novel_182, novel_46 and novel_64 are
   also different between ‘Fengzao’ and ‘Kyoho’. The difference in
   expression profiles may suggest that the mutations in ‘Fengzao’
   modified the miRNAs expression during berry development.

   The main targets of miR398 is copper/zinc superoxide dismutase
   (Cu/Zn-SOD, CSD), a scavenger enzyme of ROS (reactive oxygen species),
   which is related to the oxidative stress [[160]39, [161]40]. The
   overexpression of miR398 had led to the downregulation of CSD1 and CSD
   2 enzymes in rice [[162]41]. In our previous study [[163]8], RNAseq
   analysis revealed that SOD was one of the most significantly
   differently expressed genes between ‘Fengzao’ and ‘Kyoho’ during berry
   development and the overall expression level of SOD in ‘Fengzao’ was
   lower than that in ‘Kyoho’, except at the stage of veraison. The
   expression of miR398 increased more in ‘Fengzao’ than in ‘Kyoho’
   (subcluster 10 of Fig. [164]5) which also verified our previous results
   [[165]8]. And the activity of SOD enzyme in ‘Fengzao’ indeed was lower
   than that in ‘Kyoho’ except at the stage of veraison [[166]10].

   The target genes of miR169 were Jasmonate ZIM Domain (JAZ) and nuclear
   transcription factor Y subunit A-3 (NFYA-3) in tomato [[167]20,
   [168]41]. Zeng et al. [[169]18] suggested that miR169 involved in the
   fruit development of Lycium barbarum, it expressed at a higher level in
   fruit than in flowers and leaves [[170]42]. Ripening can be considered
   as a stressful process with a progressive increase in oxidation
   [[171]43, [172]44]. Overexpression of miR169 caused significantly
   down-regulation of its target genes and induced the increased drought
   tolerance of tomato [[173]45]. Accompanied the gradually increasing
   expression of miR169 in ‘Fengzao’, and the overall oxidation status
   (ROS levels) in ‘Fengzao’ was indeed higher than that in ‘Kyoho’
   [[174]10] .

Differentially expressed miRNAs between ‘Fengzao’ and ‘Kyoho’

   Only miR159a was significantly differentially expressed between any
   adjacent stages during the berry development either in ‘Fengzao’ or
   ‘Kyoho’, but it was not significantly differentially expressed between
   ‘Fengzao’ and ‘Kyoho’, which suggested that miR159a constitutively
   expressed along the grape berry development. Based on the expression
   analysis of miRNAs, the significantly differentially expressed miRNAs
   between ‘Fengzao’ and ‘Kyoho’ were novel_144, vvi-miR3626-3p,
   vvi-miR3626-5p, vvi-miR399b, vvi-miR399e and vvi-miR477b-3p. Among
   these miRNAs, novel_144, vvi-miR3626-3p, vvi-miR3626-5p only expressed
   in ‘Kyoho’ and their expression patterns were entirely distinct between
   ‘Fengzao’ and ‘Kyoho’ (Fig. [175]5: subcluster1). vvi-miR3626 was
   previously not detected in grape [[176]21, [177]22].

   It is known that miR159 target the GAMYB transcription factors to
   affect flowering time [[178]46]. The miR159GAMYB system is conserved
   for the vegetativetoreproductive phase transition in plants [[179]40].
   Grape targets of miR159 include the GAMYB transcription factors MYB33,
   MYB65, and MYB101 [[180]27], which participate in the signaling process
   induced by ABA accumulation in the presence of stress [[181]47].The
   up-regulation of miR159 could inhibit ABA signaling through
   down-regulation of MYB transcription factors [[182]47]. Abscisic acid
   has key role for the onset of grape berry ripening [[183]48].
   Therefore, it could be comprehensible for the differentially expressed
   of miR159 during the berry development both in ‘Fengzao’ and ‘Kyoho’.
   Similarly, miR159b and miR319a/b/c/d had high expression level at early
   development stages, and were down-regulated at the late stages of fruit
   ripening in melon [[184]12]. Fa-MIR159a transcript reached its highest
   expression level during the green stage and subsequently decreased
   significantly during the white and red stages strawberry [[185]49].

   The targets of vvi-miR3626-3p, vvi-miR3626-5p in this study are
   predicted as DEAD-box ATP-dependent RNA helicase, fkbp12-rapamycin
   complex-associated protein, Serine/threonine- protein kinase and
   phosphate phosphoenolpyruvate translocator, etc. (Table [186]3). The
   DEAD-box RNA helicases are the largest family of RNA helicases, and
   DEAD-box helicase is believed to play crucial roles in plant growth and
   development [[187]50]. Previous studies have shown that the
   overexpression of helicase genes could elevate SOD activity [[188]51,
   [189]52]. In this study, the lacking expression of vvi-miR3626-3p and
   vvi-miR3626-5p in ‘Fengzao’ suggested the SOD activity should decrease
   accordingly. Interestingly, our previous results indeed confirmed this
   [[190]8, [191]10]. SOD is an important ROS scavenging enzyme integral
   to plant stress tolerance [[192]53]. Two DEAD-box helicases function in
   abscisic acid (ABA)-dependent and ABA-independent abiotic stress
   signaling pathways in Arabidopsis thaliana [[193]54]. Abscisic acid is
   a major regulator of grape berry ripening [[194]48] and ABA also
   involves in plant adaptive responses to abiotic and biotic stresses
   processes [[195]55]. The functional involvement of a putative helicase
   in the antioxidative responses in alfalfa has also been reported
   [[196]52]. Wu et al. [[197]3] proposed that stress response process may
   play an important role during citrus fruit ripening when the
   comparative analyses of miRNAs expression between a spontaneous
   late-ripening sweet orange mutant and its wild-type were conducted
   using small RNA sequencing. Pilati et al. [[198]56] demonstrated an
   oxidative burst in ‘Pinot Noir’ at veraison stage and ROS could
   participate to the regulatory network of fruit development in grape
   [[199]57].

   Therefore, taken the above mentioned expression patterns of vvi-miR398,
   vvi-miR169b and vvi-miR159a; the differentially expression of
   vvi-miR3626-3p and vvi-miR3626-5p; and our previous results [[200]8,
   [201]10] into consideration, we deduced that the stress response (ROS
   related) play a significant role during grape berry ripening and
   partially accounted for the early-ripening mutant of ‘Fengzao’.
   Furhtermore, H[2]O[2] as the exogenous ROS stress indeed promoted the
   early-ripening of ‘Kyoho’ [[202]10].

   miR399 involved in plant responses to phosphate starvation by targeting
   an ubiquitin conjugating enzyme (UBC) gene associated with the function
   of MYB transcription factor [[203]58]. UBC24, functions as a repressor
   that prevents excessive accumulation of Pi [[204]59]. In this study,
   vvi-miR399b and vvi-miR399e targeted to high affinity inorganic
   phosphate transporter. Phosphorus (Pi) is an essential nutrient for
   optimal plant growth and development. The presence of miR399 in
   response to phosphorus deficiency has been confirmed in various plants
   [[205]58]. miR399 family members were also differentially expressed in
   different grape genotypes, which showed they were connected to grape
   berry development [[206]23]. Okamoto et al. [[207]60] showed that
   phosphate enhanced ROS production. Jia et al. [[208]61] proposed that
   overexpression of AtMPT3 elevated the Pi concentration in mitochondrial
   matrix, which accelerated the subsequent processes of electron
   transport, ATP biosynthesis, ROS accumulation and PCD. In this study,
   both vvi-miR399b and vvi-miR399e were down-regulated in ‘Fengzao’,
   which should theoretically result in the elevation of Pi content which
   may further facilitate the ROS production to promote the early-ripening
   of ‘Fengzao’ based on the above assumptions.

   miR477 was found to target GRAS family transcription factor in grape
   [[209]21, [210]22]. GRAS involved in many processes of plant growth and
   development, as well as in plant disease resistance and abiotic stress
   responses [[211]62]. Grimplet et al. [[212]62] revealed the possible
   functions of GRAS genes in grape development and stress responses.
   Several grape GRAS genes showed differential expression among different
   berry ripening stages [[213]63]. In this study, vvi-miR477 was
   upregulated in ‘Fengzao’, which suggested its connection with the
   early-ripening of ‘Fengzao’.

   Candidate miRNAs that showed remarkably differential expression levels
   were further confirmed by Quantitative RT-PCR (qRT-PCR). Totally, 9
   known- and 1 novel- miRNAs were chosen for qRT-PCR confirmation; most
   of the results are consistent with the sequencing data (Fig. [214]6).

A possible network may contribute to the regulation of early-ripening of
‘Fengzao’

   Based on the expression profiles and differentially expressed miRNAs
   between ‘Fengzao’ and ‘Kyoho’, miR159-MYBs, miR169-NF-Y subunit,
   miR398-CSD, miR3626- RNA helicase, miR399-phosphate transporter, and
   miR477-GRAS transcription factor were integrated to further explore
   their potential roles for early-ripening of ‘Fengzao’ (Fig. [215]7).

Fig. 7.

   Fig. 7
   [216]Open in a new tab

   The speculated possible regulatory network of miRNAs-TFs-Genes in
   early-ripening of ‘Fengzao’. PDHE: pyruvate dehydrogenase, CSD:
   superoxide dismutase, ROS: reactive oxygen species, NFYA: nuclear
   transcription factor Y subunit A, PHO2: PHOSPHATE2, SPL:
   squamosa-promoter binding protein-like

   miR159 acts as a negative regulator, while AtMYB33 and AtMYB101 act as
   positive regulators, of the ABA response [[217]47]. MYB transcription
   factor was able to activate miR399, which responded to phosphate (Pi)
   starvation in Arabidopsis [[218]58, [219]64]. Upon Pi deprivation, the
   expression of miR156, miR399, miR778, miR827, and miR2111 were induced,
   whereas the expression of miR169, miR395, and miR398 was repressed
   [[220]65]. During Pi starvation, the upregulated miR399 regulates Pi
   homeostasis by suppressing the expression of UBC gene [[221]65].

   In Arabidopsis, SPL3 could directly activate the expression of several
   Pi starvation inducible genes. A miR156-SPL3 module involved in
   negative feedback regulation of miR399-mediated pathways in response to
   low P was proposed [[222]66]. SPL7 is essential for the expression of
   miR397, miR398, miR408, miR857 [[223]67]. miR169 which targets multiple
   members of subunit A of the NF-Y transcription factor complex, was
   downregulated by N or P starvation [[224]68] and was also downregulated
   by ABA [[225]69]. The PAT1 branch of GRAS family was markedly induced
   by ABA [[226]70]. NtGRAS1 from tobacco was strongly induced by various
   stimulants that raise the intracellular reactive oxygen (ROS) levels
   [[227]71].

   Based on the above progress and profiles of miRNAs detected in this
   study, we presumably constructed the network of miRNAs which regulating
   early-ripening of ‘Fengzao’ as Fig. [228]7 showed. The final connecting
   point of these miRNAs could be concentrated on the regulation of ROS
   level. We have shown the effect of ROS for accelerating the
   early-ripening of ‘Kyoho’ [[229]10], which strongly support the
   conclusions of this study.

   These results indicated that the regulation network of miRNAs for fruit
   ripening are quite complicated, and elucidation of the molecular
   mechanisms underlying the interplay between miRNA and their target
   genes involved in early-ripening of grape requires further study.

Conclusions

   In the present study, we performed a comparative analysis of miRNA
   expression profiles between ‘Fengzao’ and ‘Kyoho’ along grape berry
   development. Our results revealed that specific miRNAs were
   differentially regulated during berry development. The expression
   patterns of several key miRNAs and their target genes during grape
   berry development and ripening stages was examined. According to the
   expression analysis and previous reports, miR169-NF-Y subunit,
   miR398-CSD, miR3626- RNA helicase, miR399-phosphate transporter and
   miR477-GRAS transcription factor were selected as the candidates for
   further investigations of miRNA regulation of early-ripening of
   ‘Fengzao’. The results of this study contribute to the understanding of
   the role of miRNAs in grape berry early-ripening and will provide new
   sights for small RNA engineering approach to grape molecular breeding.

Methods

Plant materials

   The berries of ‘Fengzao’ and ‘Kyoho’ at different stages were collected
   in 2016 from the farm of Henan University of Science & Technology,
   Luoyang, China. Berries from 5 individual vines were pooled and
   immediately frozen in liquid nitrogen, then stored at − 80 °C until
   further use. Based on the comprehensive considerations of the grape
   growth stages which Coombe [[230]72] established and our previous
   results [[231]6–[232]9], the berries at periods of E-L 33 (hard green
   berries), E-L 34 (starting to soften), E-L 35 (véraison), E-L 37 (sugar
   and anthocyanins accumulation, and active growth, Berries not quite
   ripe) were sampled. Because the development status of the berry was
   largely different between ‘Fengzao’ and ‘Kyoho’, especially for the
   interval between E-L33 and E-L34 in ‘Kyoho’, the sampling points in
   ‘Kyoho’ are one more than ‘Fengzao’. The specific sampling time and
   corresponding development stages are shown in Table [233]1.

Small RNA library construction and sequencing

   Total RNA was extracted from the berry according to the method of
   Rienth et al. [[234]73]. The RNA purity and integrity were evaluated
   using Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA)
   and quantified by Qubit® 2.0 Flurometer (Life Technologies, CA, USA). A
   total amount of 3 μg total RNA per sample were used as input material
   for the small RNA library construction. Sequencing libraries were
   generated using NEBNext® Multiplex Small RNA Library Prep Set for
   Illumina®(NEB, USA.) following manufacturer’s recommendations, then
   sequencing on an Illumina Hiseq 2500/2000 platform at Novogene
   Corporation (China).

Bioinformatics analysis of the sequencing data

   The raw reads were first cleaned up and filtered by removing adapter
   sequences and low-quality reads containing ploy-N, with adapter
   contaminants or the insert tag, and lengths < 18 nt. High-quality clean
   small RNA tags were aligned to grape reference sequence
   ([235]ftp://ftp.ensemblgenomes.org/pub/plants/release-23/fasta/vitis_vi
   nifera) with Bowtie 0.12.9. Then, the clean sequences were annotated by
   searching against RepeatMasker and Rfam 11.0 databases to determine
   repeat, rRNA, tRNA, snRNA and snoRNA. The remaining sequences were
   mapped to miRBase 20.0 to identify known miRNAs. Both miREvo [[236]74]
   and mirdeep2 [[237]75] were employed to predict novel miRNA through
   exploring the secondary hairpin structure, the Dicer cleavage site and
   the minimum free energy of the small RNA tags unannotated in the former
   steps. At the same time, custom scripts were used to obtain the
   identified miRNA counts as well as base bias on the first position with
   certain length and on each position of all identified miRNA
   respectively. The miRNA editing sites located on the known miRNA
   sequences were analyzed by aligning all the sRNA tags to mature miRNA,
   allowing one mismatch.

miRNA expression and comparison between ‘Kyoho’ and ‘Fengzao’

   miRNA expression levels were normalized by TPM (transcript per million)
   with the criteria: Normalized expression = mapped read count/ Total
   reads*1000000. DESeq package 1.12.0. was used for differential
   expression analysis of miRNAs between ‘Kyoho’ and ‘Fengzao’ at the
   whole level including all the berry development stages. The P-values
   were adjusted using the Benjamini & Hochberg method. Corrected P-value
   of 0.05 and |log[2](foldchange) > 1| were set as the threshold for
   significantly differential expression. The expression change patterns
   of individual miRNAs were characterized by K-means clustering in both
   ‘Kyoho’ and ‘Fengzao’ based on log 10 (1 + TPM) normalization. All
   statistical analyses were performed in R version 3.4.3.

miRNA-target prediction and enrichment analysis

   Predicting the target genes of miRNA was performed in psRobot software
   [[238]76]. The target gene candidates of differentially expressed
   miRNAs were used for Gene Ontology (GO) enrichment analysis.
   ClusterProfiler 3.8.0 was implemented for GO and Kyoto Encyclopedia of
   Genes and Genomes (KEGG) pathway enrichment analysis [[239]77]. GO
   terms or pathways that had a p value less than 0.05 after Bonferroni
   correction were scored as significant.

Validation of miRNA expression with qRT-PCR

   To validate the expression of the differentially expressed miRNAs and
   their role in grape berry development, differently expressed miRNAs
   were chosen to design miRNA-specific stem-loop primers for stem-loop
   qPCR analysis. The miRNA specific forward primers and stem-loop RT
   primers were designed with the primer premier 5.0 software. All primers
   sequences were provided in Additional file [240]11: Table S7. Stem-loop
   qRT-PCR was performed to validate the expressions of miRNAs with three
   biological replicates based on a previous method [[241]3]. The samples
   were newly collected from three individual vines in 2017; each was a
   biological replicate. cDNAs were reverse transcribed from total RNAs
   using SuperScript III Reverse Transcriptase (Invitrogen, USA). U6 and
   ubiquitin1 gene were used as the internal control for qRT-PCR of miRNAs
   and target genes. For each reaction, 0.5 μL cDNA and 5 μL SYBR Green
   PCR Master mix (Toyobo, Osaka, Japan), and 0.5 μM forward and reverse
   primers were mixed. PCR was run in triplicate at 95 °C for 1 min,
   followed by 40 cycles at 95 °C for 15 s and 60 °C for 1 min with CFX96
   Real Time PCR Detection System (Bio-Rad, USA). Meanwhile, the
   expression profiles of nine predicted target genes were also tested by
   qRT-PCR with primers listed in Additional file [242]6. The relative
   expression changes of mature miRNAs and predicted target genes were
   calculated using the 2^-ΔΔCt method.

Additional files

   [243]Additional file 1:^ (9.9KB, xlsx)

   Table S1. The summary of the sequenced reads of different samples.
   (XLSX 9 kb)
   [244]Additional file 2:^ (12.1KB, xlsx)

   Table S2. The composition of each miRNA family identified in this
   study. (XLSX 12 kb)
   [245]Additional file 3:^ (14.8KB, xlsx)

   Table S3. The mature sequenced of all known and novel miRNAs. (XLSX 14
   kb)
   [246]Additional file 4:^ (115.4KB, xlsx)

   Table S4. The identified miRNAs in ‘Fengzao’ and ‘Kyoho’, respectively.
   (XLSX 115 kb)
   [247]Additional file 5:^ (690.2KB, pdf)

   Figure S1. Deeply sequence conserved and previously reported miRNA
   families detected from developing berries of ‘Fengzao’ and ‘Kyoho’.
   miRNA families (columns) are conserved between plants families (rows)
   for plant species represented in miRBase release 21. The bar represents
   the scale of the numbers of the miRNAs after the Z-score
   standardization. (PDF 690 kb)
   [248]Additional file 6:^ (10.5KB, xlsx)

   Table S5. Target gene prediction results of all the miRNAs identified
   in this study. (XLSX 10 kb)
   [249]Additional file 7:^ (7.3MB, pdf)

   Figure S2. miRNA variants and their nucleotide bias position. A: MiRNA
   nucleotide bias at each miRNA position. B: First nucleotide bias for
   the first position of 18- to 30-nt miRNAs. Relative nucleotide bias at
   each miRNA position compared with the total RNA. (PDF 7480 kb)
   [250]Additional file 8:^ (2.7MB, pdf)

   Figure S3. GO categories and distribution of miRNA targets gene of all
   the miRNAs identified in this study. The left-hand-side scale is the
   percent of the targeted gene numbers corresponding to the GO terms.
   (PDF 2795 kb)
   [251]Additional file 9:^ (6.4KB, pdf)

   Figure S4. KEGG analysis of the 20 most enriched pathways. The coloring
   of the q-values indicates the significance of the rich factor. The
   circle indicates the target genes that are involved, and the size is
   proportional to the gene numbers. The x-axis represents name of
   enrichment pathway. The Y-axis represents rich factor. (PDF 6 kb)
   [252]Additional file 10:^ (42.4KB, xlsx)

   Table S6. The expression data (TPM value) of known and novel miRNAs and
   the location of clusters based on the K-means method in this study.
   (XLSX 42 kb)
   [253]Additional file 11:^ (10.6KB, xls)

   Table S7. The primers of miRNAs and genes used for qRT-PCR
   verification. (XLS 10 kb)

Acknowledgements