Identification of miRNAs-mediated seed and stone-hardening regulatory networks and their signal pathway of GA-induced seedless berries in grapevine (V. vinifera L.)

Background Stone-hardening stage is crucial to the development of grape seed and berry quality. A significant body of evidence supports the important roles of MicroRNAs in grape-berry development, but their specific molecular functions during grape stone-hardening stage remain unclear. Results Here, a total of 161 conserved and 85 species-specific miRNAs/miRNAs* (precursor) were identified in grape berries at stone-hardening stage using Solexa sequencing. Amongst them, 30 VvmiRNAs were stone-hardening stage-specific, whereas 52 exhibited differential expression profiles during berry development, potentially participating in the modulation of berry development as verified by their expression patterns. GO and KEGG pathway analysis showed that 13 VvmiRNAs might be involved in the regulation of embryo development, another 11 in lignin and cellulose biosynthesis, and also 28 in the modulation of hormone signaling, sugar, and proline metabolism. Furthermore, the target genes for 4 novel VvmiRNAs related to berry development were validated using RNA Ligase-Mediated (RLM)-RACE and Poly(A) Polymerase-Mediated (PPM)-RACE methods, and their cleavage mainly occurred at the 9th–11th sites from the 5′ ends of miRNAs at their binding regions. In view of the regulatory roles of GA in seed embryo development and stone-hardening in grape, we investigated the expression modes of VvmiRNAs and their target genes during GA-induced grape seedless-berry development, and we validated that GA induced the expression of VvmiR31-3p and VvmiR8-5p to negatively regulate the expression levels of CAFFEOYL COENZYME A-3-O-METHYLTRANSFERASE (VvCCoAOMT), and DDB1-CUL4 ASSOCIATED FACTOR1 (VvDCAF1). The series of changes might repress grape stone hardening and embryo development, which might be a potential key molecular mechanism in GA-induced grape seedless-berry development. Finally, a schematic model of miRNA-mediated grape seed and stone-hardening development was proposed. Conclusion This work identified 30 stone-hardening stage-specific VvmiRNAs and 52 significant differential expression ones, and preliminary interpreted the potential molecular mechanism of GA-induced grape parthenocarpy. GA negatively manipulate the expression of VvCCoAOMT and VvDCAF1 by up-regulation the expression of VvmiR31-3p and VvmiR8-5p, thereby repressing seed stone and embryo development to produce grape seedless berries. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03188-y.

Conclusion: This work identified 30 stone-hardening stage-specific VvmiRNAs and 52 significant differential expression ones, and preliminary interpreted the potential molecular mechanism of GA-induced grape parthenocarpy. GA negatively manipulate the expression of VvCCoAOMT and VvDCAF1 by up-regulation the expression of VvmiR31-3p and VvmiR8-5p, thereby repressing seed stone and embryo development to produce grape seedless berries.
Keywords: Grape, Seed and stone hardening development, GA, Seedless berry, miRNA, Target genes Background Grape (V. vinifera L.) is a soft pulpy berry with a thin edible outer skin (exocarp) and fleshy edible inner layers (mesocarp and endocarp) of storage tissues containing seeds. During grape stone-hardening stage, the coats surrounding the seed gradually harden to form a lignified seed coat (stone). Stone hardening is an essential strategy for seed protection and dispersal in different plant species, including cherry (Prunuscerasus and P. avium), peach (P. persica), plum (P. salicina) and grape [1][2][3]. Lignin deposition plays a critical role in seed stone-hardening formation, and several transcriptome and proteome studies have demonstrated that the flavonoid and lignin biosynthesis pathways are highly involved in the lignification of seed-coat structures [4][5][6][7]. However, seed coat development might affect the expansion and ripening of fruits, resulting in low-quality fruit [7]. By contrast, seed abortion and seed-coat degradation during the seed-development stage could lead to seedless berry, which is a favourable trait for consumers [8][9][10]. Therefore, an in-depth understanding of the regulatory mechanisms and molecular basis underlying seed stonehardening formation during grape-berry development is essential for the production of high-quality fruit.
With the recent development of high-throughput sequencing technologies, many miRNAs from various plant species have been released on the miRBase21.0 database (http://www.mirbase.org). The miRNA-guided cleavage of target mRNAs and/or the translational repression of development-related genes especially transcription factors (TFs) in different plant species, including grapevine plays an important role in the regulation of fruit development and ripening [11][12][13]. For example, PbrmiR397a regulates the fruit stone-cell lignification by inhibiting the expression of three LAC-CASE (LAC) genes involved in lignin biosynthesis, resulting in decreased lignin content and stone-cell number in Chinese pear (Pyrus bretschneideri) fruit [7]. Arabidopsis AUXIN RESPONSE FACTOR8 (AtARF8) and FRUITFULL (AtFUL) MADS-domains act together to directly activate the expression of MIR172C, a valvespecific AtmiR172-encoding gene, leading to the repression of the flower-patterning gene APETALA2 (AtAP2) and the promotion of fruit valve growth [14]. Similarly, tomato (S. lycopersicum) SlmiR156/157 and SlmiR172 have been reported as important regulators in the ripening process of tomato fruits by inhibiting the expression of ripening regulatory genes COLORLESS NON-RIPENING (CNR) and SlAP2a [15]. Our research group also identified and characterized many known and grape specific miRNAs (VvmiRNAs) during grape-berry development, of which a large number of VvmiRNAs might be involved in the modulation of berry development and seed formation [12,13,[16][17][18][19], and a good example is that VvmiR058 negatively regulates the expression of POLYPHENOL OXIDASE (VvPPO) gene involved in the synthesis of lignin in peel and seeds during berry development [8]. Nevertheless, the research on the regulatory networks and transcriptome dynamics of VvmiRNAs during seed and stone-hardening development in grape remains imperative.
Gibberellin (GA) is a well-known phytohormone involved in diverse biological processes of grape-berry development, leading to the improvement of berry size, weight, and seedless-berry formation [13,20]. Many VvmiRNAs are differentially expressed during the different stages of grape-berry development in response to GA application [21]. Consequently, the seed stone-hardening stagespecific VvmiRNAs responsive to GA signal could be further identified, which is crucial to the molecular breeding of seedless-grape production. Hence, in the present study, we attempted to identify and characterize the expression of seed stone-hardening stage-specific VvmiRNAs and their target genes using Solexa sequencing technology. Moreover, GO and KEGG pathway analyses were performed to explore the regulatory networks of VvmiRNAmediated seed stone-hardening development stage in grape berries. The dynamic regulatory roles of several important miRNA-mediated target genes involved in the GA signal pathway and seed development in GA-induced grape seedless berries were also explored. Our results can provide novel genetic information for the improvement of grape breeding programs to advance the production of new seedless-grape varieties.

Results
High throughput sequencing of small RNA in grape berries of the stone-hardening stage The sRNA library from the 'Wink' cultivar grape stonehardening berries (SB) at 45 days after flowering (DAF) was constructed and sequenced in depth by Solexa technology to explore the role of sRNA in the grape stonehardening process (Fig. 1a). Sequencing result analysis showed that total of 11,456,656 redundant and 3,591,857 unique clean reads were identified (Supplementary Table  S1), of which about 18.40 and 8.79% were annotated, respectively. Out of them, 5.54 and 1.99% of these two reads, respectively, were further successfully mapped into the non-coding RNA of rRNAs, snRNAs, snoRNAs, and tRNAs. While, 5.71 and 0.04% of the redundant and unique reads, respectively, were determined to be putative known miRNAs, and the majority of the redundant (81.6%) and unique (91.21%) reads, were mapped to unannotated regions in the grape genome (Fig. 1b).
The length distribution of sRNAs in the SB library were uneven, and most sRNA reads were 21 and 24 nucleotides (nt), which were the characteristic lengths of miRNAs and siRNAs, respectively, and consistent with the expected size range generated by Dicer [22]. In contrast to two other sRNA libraries of 5 DAF berry (YB) and 90 DAF berry (MB) (data were not shown), the sRNA distribution at 21 and 24 nt at SB was more similar to that at MB than that at YB (Fig. 1c). The comparative expression level using normalized counts per million of miRNA (21 nt) and siRNA (24 nt) at grape YB, SB, and MB stages indicated a gradual increase in the known miRNA (21 nt) reads towards the MB stage. By contrast, the siRNA (24 nt) reads exhibited a reverse trend with the highest level observed at the early YB stage followed by a gradual decrease towards the MB stage (Fig. 1c). The above results suggested that known miRNA (21 nt) reads were highly abundant during the later stage (MB) of grape-berry ripening compared with siRNA (24 nt) reads.
Identification and characterization of known VvmiRNAs in the stone-hardening stage of grape berries A total of 143 known VvmiRNAs and 18 corresponding precursors, namely, VvmiRNAs* belonging to 48 VvmiRNA families, were identified in the stonehardening stage of grape berries (Supplementary Table  S2). Although the number of known VvmiRNA members within each family varied from 1 to 24 during grape-berry development, they had the most variation at the SB stage compared with the other two stages (YB and MB) ( Fig. 1d and Fig. 1e). Amongst 48 VvmiRNA families, VvmiR169 family was highly represented during all three stages with 22-24 members, followed by VvmiR395 family with 13 members, VvmiR156 family with 8-9 members, and VvmiR166 family with 8 members, whereas the remaining miRNA families consisted of 1-7 members (Fig. 1d). This result suggested a diversification in the functions of these VvmiRNA families during grape development and ripening. Interestingly, from the percent numbers of VvmiRNA members at different grape-berry developments, we found that VvmiR3628 family was sequenced only in mature berries (90 DAF), whereas VvmiR393 family was detected only in young berries (5 DAF) (Fig. 1e), indicating the stagespecificity of VvmiRNAs' expression.
We further characterized most known VvmiRNA families with high expression abundances at the stonehardening stage of grape berries. For instance, amongst the 143 known VvmiRNAs, about 68% exhibited high copy read number, including 42 known VvmiRNAs with read number > 1000, whereas 55 known VvmiRNAs had read numbers ranging between 100 and 1000 (Supplementary Table S2). Specifically, VvmiR166, VvmiR168, VvmiR479, VvmiR156, and VvmiR3636 families (not including their VvmiRNA* sequences) possessed more than 10,000 copy reads (Supplementary Table S2). Meanwhile, 83.3% of known VvmiRNA* family members exhibited lower copy number than their corresponding VvmiRNAs ( Fig. 1f and g, and Supplementary Table S2), which may have been due to the fact VvmiRNAs* were easier to degrade than their corresponding VvmiRNAs. Thus, they usually had less copy numbers, similar to a previous report [23]. Nevertheless, two VvmiRNA* families of VvmiR3623* and VvmiR2950* showed higher copy numbers (14,830 and 3029, respectively) than their corresponding mature sequences (4405 and 1773, respectively) ( Fig. 1g and Supplementary Table S2). This finding implied that VvmiR3623* and VvmiR2950* might play potential roles during the development and ripening of grape berries, similar to VvmiRNAs [17,24].
Screening of novel VvmiRNAs at grape stone-hardening stage and their validation using miR-RACE and qRT-PCR According to the annotation criteria of novel miRNAs [25], all un-annotated sRNAs were used to explore the stem-loop structures of their precursors for the prediction of novel miRNAs. A total of 90,352 reads were identified as novel VvmiRNAs, including 72 novel VvmiRNAs and 12 novel VvmiRNAs*, in the grape stone-hardening stage (Table 1). These novel precursors were folded into stable hairpin structures, and their negative minimal folding free energy (MFE) ranged from − 108.2 kcal mol − 1 to − 20.88 kcal mol − 1 (Table 1), which was in line with the criteria of novel VvmiRNAs (MFE < − 20.0 kcal mol − 1 ) as previously reported [25]. The novel VvmiRNAs and VvmiRNAs* were primarily 21 nt in length, accounting for 84.14% (69/82), and the first base with uracil (U) at the 5′-end of their mature sequences reached 56.0%, confirming that these were novel VvmiR-NAs (Fig. 2a). From our datasets, the novel VvmiRNAs were unconserved, species specific, and low abundance; they usually exhibited lower accumulation level than conserved ones, and in agreement with previous results [23,26,27]. Interestingly, we also observed a few novel VvmiRNAs with high abundance and read number > 1000, such as VvmiR10, VvmiR13, VvmiR29, VvmiR30, VvmiR34, VvmiR37, VvmiR43, and VvmiR71 (Table 1). Amongst them, VvmiR37, VvmiR13, and VvmiR30 showed higher copy numbers (37,985,13,748, and 10,199, Fig. 1 The sRNA dataset of grape berries at 45DAF stage. a Morphology of grape berries at the stone-hardening stage (SB). b Pie charts of diverse types of sRNAs derived from grape berry library at seed coat-hardening stage based on the redundant and unique reads, respectively. c Expression level of miRNA (21 nt) and siRNA (24 nt  respectively) than the others (Table 1). Furthermore, their corresponding VvmiR37* and VvmiR13* exhibited high copy numbers with 2588 and 1332, respectively, which implying that VvmiR37/VvmiR37* and VvmiR13/VvmiR13* may play significant roles during the stone-hardening stage of grape berries. The location of novel VvmiRNAs in their precursors showed that 40 novel VvmiRNAs were located in the 5′-arm of their precursors, whereas another 32 ones were located in the 3′-arm of their precursors (Table 1). Similarly, seven novel VvmiRNAs* were located in the 5′-arm of their precursors, whereas another 5 ones were located in the 3′-arm of their precursors (Table 1). These results indicated that the 5′-arm of miRNA precursors may be more efficient in generating miRNAs and miRNA* than the 3′-arm. However, further confirmation of this mechanism is necessary. The identified novel VvmiRNAs exhibited three types of VvmiRNA-3P or VvmiRNA-5P or the first sequences (3p and 5p) ( Table 1). As shown in Table 1, the miRNA with 3p indicated this miRNA sequence originated only from the 3′ arm of its precursor, and the marked 5p denoting the corresponding miRNA originated just from the 5′ arm of its precursors. Conversely, the miRNA with 3p and 5p represented both arms of miRNA precursors that generated two sequences of miRNA and miRNA* ( Table 1). The distribution of all these novel VvmiRNAs varied between the 19 grape Chrs and 1 unknown Chr (Fig. 2b). Amongst them, Chr19 possessed the highest number (10) of novel VvmiRNAs, followed by Chr8 with 9 novel VvmiRNAs and the unknown Chr with 8 ones (Fig. 2b), whereas Chr3, Chr4, and Chr11 did not harbor any novel miRNAs (Fig. 2b). Subsequently, six novel VvmiRNAs closely related to berry development, such as VvmiR8, VvmiR16, VvmiR31, VvmiR38-5p, VvmiR44-3p, and VvmiR53-3p, were further validated by miR-RACE. Their precise sequences were detected by miR-RACE, and the sequences were consistent with those from the high-throughput sequencing dataset (Table 1 and Supplementary Table S2), which further verified the results of high throughput sequencing. Moreover, their qRT-PCR expression profiles during grape-berry development showed differential expression patterns similar to the high-throughput sequencing dataset (Fig. 2c). Therefore, our miR-RACE and qRT-PCR results confirmed the reliability and expression modes of VvmiRNA involved in the modulation of grape-berry development.

Identification of grape stone-hardening stage-specific VvmiRNAs
The identification of grape stone-hardening stagespecific VvmiRNAs is essential to gain insights into the regulatory roles of grape-berry development. Compared with our other two sRNA libraries from 5 and 90 DAF in our other work (Supplementary Tables S3 and S4 and Fig. 2d), 35 VvmiRNAs/VvmiR-NAs* were identified only at the stone-hardening stage of grape berries, including 28 VvmiRNAs (1 known and 27 novel) and 7 VvmiRNAs* (2 known and 5 novel). Among them, there were a large number of stage-specific novel VvmiRNAs indicated that novel miRNAs may play significant roles in the stonehardening stage of grape-berry development. To identify the differential expression VvmiRNAs during grape-berry development, the fold changes log2 (YB/ SB) or log2 (MB/SB) > 1 cut offs were selected, and the filtered VvmiRNA/VvmiRNA* possessed significant expression difference across diverse development stages of grape berries. Here, we discovered that 52 VvmiRNAs/ VvmiRNAs* exhibited significant differences in their expression levels during grape-berry development and ripening, comprising 44VvmiRNAs (34 known and 10 novel ones) and 8 VvmiRNAs* (5 known and 3 novel ones) (Supplementary Tables S3 and S4 and Fig. 2d). This finding indicated that they may possess dynamic regulatory roles of grape-berry development and of them, more known VvmiRNAs may have dynamic variation in their regulatory roles than novel ones. . Each bar indicated the mean ± SE of triplicate assays. H: high throughput sequencing; Q: qRT-PCR; each assay was conducted by using three repeats. There is no significant difference between high throughput sequencing and qRT-PCR results in P < 0.05 by Student's t-test. d Comparison of specific and significant VvmiRNAs in grape sRNA library. VvmiRNAs with specific expressions at the seed coathardening stage, Spec-VvmiRNAs: VvmiRNAs with the absolute log 2 Fold (YB/SB) or log 2 Fold (MB/SB) > 1, Significant difference VvmiRNAs SNPs and their edit types of known VvmiRNAs/ VvmiRNAs* from grape berries at stone-hardening stage Numerous SNP variations of known VvmiRNAs/VvmiR-NAs* and their Edit types were detected in our datasets (Fig. 3), which were consistent with our previous work in 'Amur' grape [23]. Identifying the characteristics of VvmiRNA SNP helps to recognition of the evolution of VvmiRNA and its overall roles in the process ofthe stone-hardening stage of grape berries. Amongst 161 types of VvmiRNAs/VvmiRNAs*, 71 types of VvmiRNA SNPs and corresponding Edit types were identified, while the mature sequences of the remaining 90 types of VvmiRNA remained unchanged ( Table 2 and Supplementary Table S5). Furthermore, several VvmiRNA families exhibited high SNPs amongst their members. For example, VvmiR166, VvmiR156, and VvmiR167 families exhibited the highest SNPs amongst their members. In contrast, VvmiR169 family possessed the lowest SNP amongst eight members ( Table 2). This finding suggested that divergence of conservation in the sequences amongst various VvmiRNA families ( Table 2 and Fig. 3), similar to the previous report [23].
Interestingly, we observed that diverse VvmiRNA families had diverse SNP variations in their Edit types and numbers (Fig. 3), supported by the Amur grape report [23]. To depict this phenomenon clearly, all VvmiRNA families with SNPs were further classified into several groups (Fig. 3). In group I, most members of each VvmiRNA family had an SNP, and each member with an SNP possessed multiple Edit types of SNP. For example, the VvmiR166 family (VvmiR166s) had 8 members and 282 Edit types of SNV (8,282), followed by VvmiR156s (9, 118), VvmiR167s (5,40), VvmiR164s (4,17), and VvmiR535s (3,9), respectively. From these VvmiRNA families, diverse members with sequence variations obviously exhibited divergent Edit types of SNP. Although the diverse members of one miRNA family with various precursors possessed same mature sequences (e.g., VvmiR166b/c/d/e/f/g/h, VvmiR156b/c/d, VvmiR156f/g/i, VvmiR167b/c/d/e, VvmiR164a/c/d, and VvmiR535a/b/c), they had various Edit types. For instance, VvmiR166b and VvmiR166c/d/e/f/g/h had the same mature sequences, but they possessed 35 and 43 Edit types (35 and 43, respectively), resembling VvmiR156b, VvmiR156c and VvmiR156d (10, 10, and 11, respectively); VvmiR167b, VvmiR167c, VvmiR167d and VvmiR167e (12,9,5, and 12, respectively); and VvmiR164a, VvmiR164b, VvmiR164c, and VvmiR164d  implied that VvmiRNA families with single member exhibited drastic divergence and may thus be active factors during VvmiRNA sequence evolution. In group III, some VvmiRNA families such as VvmiR169b/c/g/h/i/l/r/u, VvmiR160a/b/c/d, VvmiR399a/b/c/d, and VvmiR3629 were also revealed only one Edit type even though they had multiple members with SNPs, indicating that they may possess relatively high conservation during VvmiRNA sequence evolution. In the final group, the remaining VvmiRNA families had fewer members and Edit types (Fig. 3). All these results confirmed the diversification of VvmiRNA families in the evolution of their mature sequences. The total read number of VvmiRNAs with SNP was further revealed to reach 77,141,827, and diverse VvmiRNA families and their various members had conspicuous divergence in the number of sequences with SNP. Amongst them, VvmiR166 and VvmiR156 families had considerably more reads with SNP than the other families (Table 2). Generally, the number of SNPs in VvmiRNA families was less than that of the normal sequences. However, the SNP sequences of 21 VvmiRNAs had more read numbers than miRNAs themselves, including VvmiR156a, VvmiR156e, VvmiR156h, VvmiR160a, VvmiR160b, VvmiR169b, VvmiR169f, VvmiR169g, VvmiR169h, VvmiR169r, VvmiR169u, VvmiR171e, VvmiR3629a, VvmiR3629b, VvmiR3629c, VvmiR3631b*, VvmiR396b, VvmiR396c, VvmiR396d, VvmiR399b, and VvmiR399c, suggesting that these VvmiRNAs had stronger evolution than the others. Interestingly, compared with homologous VvmiRNAs from the grape cv. 'Pinot Noir' in miRBase 21.0 (http:// www.mirbase.org/summary.shtml?org=vvi), some VvmiRNAs could not be identified in this work. However, their SNP sequences such as VvmiR164b, VvmiR169i, and VvmiR828b were identified (Supplementary Table S5; bold and italic words). Thus, SNPs may explain the generation of new members of VvmiRNA family in miRNA evolution.
Functional annotation of specific-VvmiRNA targets during grape stone-hardening stage To further recognize the roles of VvmiRNA during grape stone-hardening stage, PsRNATarget software ( h t t p : / / p l a n t g r n . n o b l e . o r g / p s R N A T a r g e t / result?sessionid=1503987414486479) was utilized to predict the potential miRNA target genes on the basis of our previous RNA-seq data (GEO Accession: GSE77218) by using mature miRNA sequences as queries. A total of 2124 targets for known VvmiRNAs and 885 targets for novel VvmiRNAs were predicted in this work, amongst which 1639 and 1635 target genes may be targeted by 24 stone-hardening-specific VvmiRNAs and VvmiRNAs with significant expressional difference, respectively. GO and KEGG analyses were performed with these predicted target gene sequences to improve our understanding of their functions in grape stone-hardening stage. Totally, 13 VvmiRNAs might be involved in the regulation of embryo development, 11 in lignin and cellulose biosynthesis, and 28 in the modulation of hormone signaling, sugar, and proline metabolism (Supplementary Table S6). From the Fig. 4a, the top 20 of GO enrichment of the stone-hardening-specific VvmiRNA targets included 18 GO enrichment of biological process and 2 of molecular function, whereas those of VvmiRNAs with significant difference had 6 GO enrichment of biological process, 13 of molecular function, and 1 of cell component (Fig. 4b). Further comparison of the biological processes between Fig. 4a and Fig.4b, the target genes for stone-hardening-specific VvmiRNAs primarily focused on the lipid, phosphatelipid biosynthesis/metabolic process (GO: 0008610; GO: 0008654; GO: 0006644), phosphatidylinositol biosynthesis/metabolic process (GO: 0006661; GO: 0046488), and glycerophospholipid/glycerolipid biosynthetic/metabolic process (GO: 0046474; GO: 0006650; GO: 0046486). These GO pathways were closely related to seed development, supporting the modulation of the stone-hardening-specific VvmiRNAs on grape seed development. Whereas, those of VvmiRNAs with significant difference largely participated in cellular metabolic process (GO: 004237), ligin metabolic process (GO: 0009808), secondary/aromatic metabolic process (GO: 0019748; GO: 0019439), cell-wall organisation or biogenesis (GO:0071555; GO: 0071554), and peptidylproline modification (GO: 0018208), thereby providing evidence of supporting their potential regulatory roles in grape-berry development and quality formation.
Verification of target genes for novel VvmiRNAs related to berry development at the stone-hardening stage of grape berries Based on the expression profiles of novel VvmiRNAs (Fig.  2c), together with their potential functional annotation, we found four novel VvmiRNAs closely related to the berry development of VvmiR8-5p, VvmiR31-3p, VvmiR38-5p, and VvmiR53-3p. Their corresponding target genes involved in embryo and seed stone development [VIT_ 204s0008g03060, DDB1-and CUL4-ASSOCIATED FAC-TOR HOMOLOG 1 (VvDCAF1)] [ [31], which were further selected to verify their roles at the stone-hardening stage of grape berries. The cleavage interactions of these four VvmiRNAs on their target genes above at the berries of grape stone-hardening stage were verified with our modified RNA ligase-mediated 5′-rapid amplifification of cDNA ends (RLM-RACE) and developed poly(A) polymerase-mediated 3′-rapid amplifification of cDNA ends (PPM-RACE) procedures [32,33], which are high effeciency, low cost strategies for validating the true target genes for miRNAs by sequencing the 3′-/ 5′-end cleavage products of their target genes.
First, using RLM-RACE, the sequencing of the amplified 3′-end products confirmed that the VvmiRNAs' cleavage on their target genes occurred at the 9th-11th sites, amongst which the 10th was their main cleavage site with the highest cleavage frequency (Fig. 5). This result indicated the specificity of cleavage sites of miRNAs, consistent with previous findings [32,33]. Although VvmiR8 and VvmiR31 had two cleavage sites on their corresponding targets (VvDCAF1 and VvCCoAOMT), the cleavage sites of the former were at the 9th and the 10th (the highest frequency 18/20), and those of the latter were at the 10th (the highest frequency 10/16) and the 11th. Similarly, VvmiR38-5p on VvGAI1 possessed three cleavage sites at the 9th-11th, amongst which the 10th on VvGDSL possessed only one cleavage site at the 10th with the same high frequency 18/22. Next, our developed PPM-RACE was used to further confirm the target genes of VvmiR8-5p, VvmiR31-3p, VvmiR38-5p, and VvmiR53-3p and their cleavage sites. The sequencing of the amplified 5′-end products identified the same cleavage sites as those of the 3′-end sequencing in the RLM-RACE experiment, but all their 5′-end cleavage product accumulation levels and their cleavage frequency detected in PPM-RACE were lower than their corresponding 3′-end cleavage frequency examined in RLM-RACE ( Fig. S1; Fig. 5). This phenomenon may be due to the fact that the 5′-end cleavage products were more easily degraded than the 3′end cleavage products, similar to previous reports [32,33]. The consistent results of the RLM-RACE and PPM-RACE experiments demonstrated that VvDCAF1, VvCCoAOMT, VvGAI1, and VvGDSL were the true target genes of VvmiR8-5p, VvmiR31-3p, VvmiR38-5p, and VvmiR53-3p, thereby verifying their cleavage-interaction mode in the stone-hardening stage of grape berries.

Expression modes of VvmiRNAs and their target genes during GA-induced grape seedless-berry development
To gain insight into seed development during grape stone-hardening stage and GA-induced parthenocarpy process, we compared berries phenotype and seed morphology of GA-treated and untreated control (CK) 'Wink' grape cultivar at 5, 20, 45, and 90 DAF were carried out (Figs. 6a-d). Amongst them, the berries at 45 DAF in untreated control group had full seeds and hardened seed coats (Fig. 1a). GA-treated grape berries showed more distinct increased in vertical diameter of berries grains at 20 and 45 DAF as compared with untreated CK groups, but no variation were observed in horizontal diameter (Fig. 6b). Additionally, GA-treated grapes showed 99.6% seedless rate (Fig. 6d), compared with untreated CK groups, GA-treated grape seeds were drastically inhibited growth, and there were almost no growth from fruit setting to maturation. These results confirmed the profound effects of GA-induced 'Wink' grape parthenocarpy, suggesting that complicated regulatory networks might exist during GA-mediated grape berries and seed development.
To determine the long-term roles of VvmiRNAs and their target genes validated above during GA-induced grape seedless-berry development, the expression levels of VvMIR8-5p, VvMIR31-3p, VvMIR38-5p, and VvMIR53-3p and their corresponding target genes VvDCAF1, VvGAI1, VvCCoAOMT, and VvGDSL were examined in berries at 5, 20, 45, and 90 DAF, respectively. Results (Fig. 6e) showed that except for VvMIR53- 3p and its target gene, the remaining three VvMIRNAs and their target genes exhibited significant expression differences in response to GA treatments relative to untreated control plants, whereas VvMIR53-3p and its target gene VvGDSL hardly differed in response to GA treatment than the control. Notably, the former three miRNAs VvmiR8-5p, VvmiR31-3p, and VvmiR38-5p expression were strongly up-regulated by GA treatment. Specifically, VvmiR8-5p and VvmiR31-3p displayed the highest expression level at the grape stone-hardening stage (45 DAF) in response to GA treatment relative to other stages (Fig. 6e), indicating that these two miRNAs may play significant roles by responding to GA during the stone-hardening stage of grape berries. By contrast, the expression of their target genes VvDCAF1, VvCCoAOMT, and VvGAI1 exhibited strong downregulation in response to GA treatment (Fig. 6e). VvmiRNAs and their target genes above also displayed the opposite expression trends during GA-induced grape seedless-berry development, confirming that these VvmiRNAs negatively modulated their target genes' expression during this process.
Interestingly, GA treatments significantly down-regulated the expression of VvGAI1, which is a key DELLA protein negative interaction factor in the GA signal pathway (Fig.  6e). Similarly, during grape stone-hardening stage GA treatment also significantly down-regulated the key genes VvCCoAOMT and VvDCFA1, which were involved in lignin biosynthesis and embryo development, (Fig. 6e). These results indicated that GA may repress grape stone hardening and embryo development by inducing the expression of VvmiR31-3p and VvmiR8-5p to negatively regulate the expression levels of VvCCoAOMT and VvDCAF1, as a key molecular mechanism involved in the modulation of GAinduced grape seedless-berry development. However, GA exhibited no effect on the expression levels of VvmiR53-3p and VvGDSL target genes compared with the untreated control, due to the fact that VvGDSL was a gene related to cell-wall development. Meanwhile, we revealed that VvGDSL exhibited the highest expression level at 20 DAF (berry-expansion stage), implying that it may participate in the modulation of the cell-wall expansion development of young berries.

Dynamic accumulation of cleavage products of the four VvmiRNAs' target genes during GA-induced grape-berry development
Monitoring the accumulation patterns of these four VvmiRNAs' cleavage products and target genes during GA-induced grape seedless-berry development could contribute to determining the variation of their cleavage roles. Here, the 3′-and 5′-end cleavage product accumulation levels were examined by RLM-RACE and PPM-RACE, respectively. Results showed similar dynamic Fig. 6 Characterization of grape berries at the seed coat-hardening stage and Gibberellin (GA 3 )-induced grape parthenocarpic berries. a-b Variation in morphology and growth curve of grape berries derived from GA 3 -treated and untreated control (CK) plants. c-d Variation in morphology and growth curve of seeds derived from GA 3 -treated and CK control plants. Vertical diameter, VD; horizontal diameter, HD. e Expression levels of VvmiRNAs responsive to GA 3 -treatment during grape berry development [5 days after flowering (5DAF), 20DAF, 45DAF, 90DAF]. The mean and SD values were obtained from three biological samples. ANOVA test was used to identify significant differences, Asterisks indicated statistically significant differences at (*P < 0.05; **P < 0.01) as determined by Student's t-test accumulation modes of both ends-cleavage products during different stages of grape-berry development (Fig. 7), confirming the dynamic variation in cleavage roles of these VvmiRNAs on their target genes during this process. Moreover, the accumulation modes of cleavage products resembled the expression modes of the corresponding VvmiRNAs, indicating that miRNAs may be the main factors affecting their interactions. Interestingly, we found that GA evidently promoted the cleavage roles of VvmiR8-5p, VvmiR31-3p, and VvmiR38-5p on their corresponding targets by obviously up-regulating accumulation levels of the corresponding target cleavage products, while VvmiR53-3p/VvGDSL pair with nearly no change under GA treatment. In particular, the cleavage products of VvmiR31 and VvmiR8 on their corresponding targets VvCCoAOMT and VvDCFA1 were boosted at the highest level at 45 DAF (grape stone-hardening stage), which may derive from the correlation of these two targets' potential functions with embryo and stone development [28,30]. Unlike this, those of VvmiR38-5p on VvGAI1 were enhanced drastically by GA at all stages of GA-induced seedless berries used in this work. These results suggested that GA might involved in manipulating grape seedless-berry development primarily by promoting the cleavage role of VvmiR38-5p on VvGAI1 at all stages here. It is inferred that GA might repress grape stone hardening and embryo development through negatively regulating the expression levels of the lignin biosynthesis gene VvCCoAOMT and embryo development related gene VvDCAF1 by inducing the expression of VvmiR31-3p and VvmiR8-5p, as a potential key molecular mechanism involved in the modulation of GA-induced grape seedless-berry development. Fig. 7 Accumulation levels of 3′−/5′-end cleavage products of VvmiRNAs on their target genes in GA-treated and untreated control (CK) plants at different stages of grape berry development by RLM-RACE and PPM-RACE. 3'DCAF1 denotes 3′-end cleavage products of target genes by VvmiRNAs, whereas 5'DCAF1 indicates 5′-end cleavage products of target genes by VvmiRNAs. Each experiment was repeated three times. ANOVA test was used to identify significant differences, Asterisks indicated statistically significant differences at (*P < 0.05; **P < 0.01) as determined by Student's t-test

Discussion
Grape stone-hardening stage is the critical stage for seed embryo development [8]. When grape inflorescences are treated with GA at 10 days before anthesis, seed embryo development and seed coat formation of grape berries are inhibited, leading to seedless-berry formation [12]. Thus, GA is a key phytohormone regulator for seed embryo development and stone-hardening in grape. Previous studies have revealed that some VvmiRNAs are differentially expressed during grape-berry development in response to GA application [8,12,21,34]. However, research on stone-hardening stage-specific VvmiRNAs and their regulatory networks during grape-berry development remains imperative. Accordingly, we identified and characterized stone-hardening stage-specific miR-NAs and their mRNA target genes in response to GA in the current work to gain in-depth insights into the molecular basis of seedless-berry development for the molecular breeding of novel seedless-grape varieties.

VvmiRNA-mediated TFs and methylation−/acylationrelated genes in grape during seed-hardening stage
In this work, the potential target genes of stonehardening stage-specific VvmiRNAs were predicted. Out of them, some target genes were the members of TF families, such as SPL (2/6/7/9/10/ , GATA (5/24), TCP (2/9), GLK 2, and WRKY20 (Fig. 8), indicating VvmiRNA could regulate berry development including seeds through regulating expression of these TFs at stone hardening state. Similar cases were reported by previous studies [35]. For example, ZmMYB138 and ZmMYB115 for ZmmiR159k regulated the endosperm development of corn seeds by affecting the transcriptional activities of Du1/Wx and Ae1/Bt2 genes [36]. MYB89 was expressed predominantly in developing seeds during maturation, which inhibited seed oil accumulation by directly repressing WRI1 [37]. AP2, WRKY TF SHB1, WRKY10 and LEC1 controlled the development of seeds in rice and Aribidopsis [38][39][40]. Additionally, VvHD-Zips and VvMADS39 were closely related to embryo abortion in grape berries [41,42]. BpMADS12 promotes lignin accumulation in B. platyphylla [43,44]. All these findings supported the view that miRNA-mediated TFs might participate in the regulatory network of seed development and seed coathardening formation (lignin) in grape.
In addition, another VvmiRNAs at grape stonehardening stage were predicated to target some methylation/acylation-related genes. As shown in Fig. 8, VvDNMT2, VvEHMT, VvMT13, VvCCoAOMT1, VvPAT, VvAcyl, and VvLSD were predicated as the potential target genes for VvmiR28, VvmiR31, VvmiR40, VvmiR56, VvmiR166, VvmiR2950, VvmiR394a/b/c, VvmiR3633a/b, VvmiR396c/d, and VvmiR399e, suggesting these VvmiRNAs above might regulate their potential target genes related to methylation/acylation, whereby involving in modulatiion of grape berry development process. Similarly, previous studies reported that DNA methylation plays important roles in modulation of fruit development [45,46]. Specific tissues have their own methylome patterns, the epigenome is not static, and extensive reprogramming occurs during tomato fruit development [47]. In apple and pear, DNA methylation of MYB transcription factor can regulate fruit pigment accumulation and fruit color [48,49]. Therefore, the current work provides the several potential genes for gaining in-depth insights into fruit development.
VvmiRNA-targeted genes involved in hormone signal pathways in grape berry at the stone-hardening stage GA and AUX play essential roles in grape berry and seed development [8,12,50]. Here, amongst our identified VvmiRNAs in grape berries at the stone-hardening stage, several VvmiRNA-mediated target genes were found to be involved in hormone metabolism and signal transduction according to functional annotations (Fig. 8). From the Fig. 8, VvGA3ox is the potential target of VvmiR3633, which is involved in the GA biosynthesis pathway, whereas VvGID (the receptor of GA signal), as the predicted target for VvmiR396, plays a key role in the GA signal-transduction pathway (Fig. 8). Similarly, VvDELLAs for VvmiR477 and VvGAMYB for VvmiR159 and VvmiR319 are the core elements of the GA signaltransduction pathway. Furthermore, GA treatment induced the expression levels of miR159s, leading to the reduction in the expression of GAMYB level, which delayed flowering, perturbed anther development, and promoted parthenocarpy and subsequently seedless-fruit formation in Arabidopsis, tomato, and grape plants [18,51]. Likewise, GA treatment down-regulated the expression levels of miR319 and miR166 during the modulation of Arabidopsis and tomato plant development [51,52].
Our previous study also revealed that VvmiR160s may mediate auxin signal transduction to regulate grape berry and seed development at the stone-hardening stage by mediating target genes ARFs [13]. Similarly, it has been reported that miR160 and miR167 have important regulatory roles in female and male reproduction and the parthenocarpy of Arabidopsis plants [51,52]. Moreover, VvAP2 and VvERF, ethylene signaling-related genes were identified as target genes for VvmiR172 and VvmiR3629 (Fig. 8). These findings suggested that VvmiRNA-mediated phytohormone signaling was an essential step during early seed development in grape. This result was in agreement with that of Curbara et al. [53], who reported that ARFs and TIR1, as auxin receptors, are negatively regulated by miR160, miR167, miR390, and miR393; and the ABA-insensitive gene ABI3 is repressed by miR516 during the early seed development of barley [53]. All these findings implied that VvmiRNAs may negatively regulate several hormone-signal-related genes during the modulation of grape berry and seed development at the stone-hardening stage.
VvmiRNA-mediated regulatory networks in berry and seed development at the grape stone-hardening stage The characterization of the potential target genes of VvmiRNAs at the grape stone-hardening stage is the key step for elucidating the miRNA-mediated regulatory networks associated with grape berry and seed development. Recently, some miRNAs and their target genes have been reported during fruit development in tomato [54], eggplant [55], pear [56], grape [23,57], citrus [58,59], banana [25], melon [60], apple [61], and strawberry [62,63]. In this work, we identified many key target genes related to berry development for VvmiRNAs Fig. 8 A mode chart of miRNA-mediated regulatory network related to grape berry seed and stone development. Green rectangle represents diverse VvmiRNAs which might regulate seed and seed coat development, whereas yellow rectangle indicates their corresponding target genes. The blue ellipse donates the metabolism pathways and specific tissues of VvmiRNAs and their target gene mediated in grape seed and seed coat development involved in sugar, acid, pigmentation and hormone metabolism, lignin synthesis, and methylation/acylation process in grape berries and seeds (Fig. 8). On the basis of the identification and characterization of VvmiRNAs and their target genes in grape berries at the stonehardening stage, a putative schematic mode of VvmiRNA-mediated berry and seed development was proposed. As shown in Fig. 8: i) some VvmiRNAs might be involved in the development of grape seeds and embryos by mediating their target genes, including VvmiR8, VvmiR27-3p, VvmiR28, VvmiR55, VvmiR56, VvmiR64, VvmiR156a, VvmiR166a-h, VvmiR171a/c/d, VvmiR319e, VvmiR3633a, VvmiR390, and VvmiR396c/d (corresponding yellow regions, same as below); ii) 19 other VvmiRNAs might participate in the modulation of lignin biosynthesis of grape stone by mediating their target genes in lignin biosynthesis, including VvmiR53-5p, VvmiR156b, VvmiR164a/c, VvmiR166a-h, VvmiR3633b, VvmiR396c/d, VvmiR397a, and VvmiR399b/c; iii) two VvmiRNAs of VvmiR19 and VvmiR44-3p might regulate seed and stone development by targeting the genes related to proline metabolism; iv) 10 VvmiRNAs might participate in the modulation of sugar metabolism, including VvmiR13-5p, VvmiR27-5p, VvmiR37-3p, VvmiR56, VvmiR63, VvmiR2950, VvmiR396c/d, and VvmiR477b; and v) six other VvmiRNAs might manipulate cellulose biosynthesis by mediating VvCSLE6 and VvCESA4, including VvmiR394a/b/c, VvmiR54, VvmiR40, and VvmiR1. Our putative schematic mode was supported by previous studies [11], which reported that miR156, miR164, miR1132, miR5077, and miR396b regulate sugar and acid metabolism in Lyciumbarbarum and pear fruits, whereas miR397a mediated-LACs manipulate the lignin synthesis of pear fruits, poplar, and Arabidopsis [7,[64][65][66]. Furthermore, miR156, miR166, miR167, miR168, miR393, miR172, and miR396 are preferentially/highly expressed during embryo development, whereas miR164 is primarily expressed in seeds [67]. These insights into the miRNA-mediated seed and stone regulatory networks in grapes could contribute to the understanding of the molecular regulatory mechanism during grape-berry development in the stone-hardening stage at the global transcriptome-wide level.
Regulatory modes of VvmiRNA-mediated seedless berries during GA-induced grape parthnocarpy GA is one of the key hormones inducing parthenocarpy to produce seedless fruit, and expression changes of genes in the GA signaling pathway could induce parthenocarpy fruit set and fruit development. Meanwhile, DELLA protein is the key repressor of GA signal transduction, and the reduction of DELLA protein activity could lead to released GA signal and the appearance of corresponding phenotypes, such as parthenocarpy [68].
Nowadays, exogenous GA is extensively used to induce grape parthnocarpy, and recent studies have shown that GA signaling induces the expression level of miR159 and miR160 to regulate grape parthnocarpy, producing seedless berries [12,13]. Here, we also found GA might repress the expression of the lignin biosynthesis enzyme gene VvCCoAOMT and the embryo developmental gene VvDCAF1 to produce grape seedless berries by upregulating VvmiR31-3p and VvmiR8-5p during GAinduced grape parthenocarpy process. With research development, more miRNAs are identified to be involved in the modulation of GA-mediated parthenocarpy, which would contribute to gain more insight into the regulatory mechanism of miRNA-mediated grape seedless berries induced by exogenous GA.

Conclusion
In the present study, a total of 161 conserved and 85 species-specific miRNAs/miRNAs* (precursor) were identified in grape berries at stone-hardening stage using Solexa sequencing. Out of them, 30 were stone hardening stage-specific VvmiRNAs. And, high SNP variations in VvmiRNA sequences resulted into the generation of new VvmiRNA family memberslike VvmiR168, VvmiR479, VvmiR3636 families and so on. Among stone-hardening stage VvmiRNAs,13 VvmiRNAs might be involved in the regulation of embryo development, 11 in lignin and cellulose biosynthesis, and 28 in the modulation of hormone signaling, sugar, and proline metabolism. The target genes for 4 novel VvmiRNAs were validated using RLM-RACE and PPM-RACE methods, and their cleavage primarily occurred at the 9th-11th sites from the 5′ ends of miRNAs at their binding regions. Furhtermore, among the 4 novel VvmiRNAs above, GA could up-regulate the expressions of VvmiR31-3p and VvmiR8-5p, thereby inhibiting those of their target genes VvCCoAOMT (the lignin biosynthesis enzyme gene) and VvDCAF1 (the embryo developmental gene) to produce grape seedless berries, as a potential key molecular mechanism involved in GA-induced grape seedless berry development. Finally, a schematic model of miRNA-mediated grape seed and stone-hardening development was proposed. Our results can serve as valuable references for the molecular breeding of seedless grape berries.

Plant materials and GA treatment
Five-year-old 'Wink' grape cultivar used as the experimental material was grown in experimental field at Jiangpu Farm of Nanjing Agricultural University, Nanjing City, Jiangsu Province, China. First, the grape berries at 45DAF stage were used for high throughput sequencing. Second, based on our preliminary trails, the total 18 inflorescence clusters with similar growth status from 6 grape plants were selected as materials, of which the 9 clusters were dipped into 50 mg/L GA 3 for 30s at 10 days before flowering to induce grape seedless berries. The other remaining 9 clusters were treated with water and used as control set. In the early morning (9 to 10 AM), 3-4 grains from the middle of each cluster of GA 3 -treated and water-treated control plants at different time points [5 days after flowering (5DAF), 20DAF, 45DAF and 90DAF] were collected and immediately frozen in liquid nitrogen and stored at − 80°C until use. Each type of samples consisted of three biological replicates.
Small RNA (sRNA) library construction and high throughput sequencing analysis Total RNA was extracted from 200 μg each sample from four grape berries at 5, 20, 45, 90DAF stages respectively using our modified CTAB method for small RNA highthroughput sequencing and qRT-PCR library construction [69]. Libraries were prepared with 1 μg total RNA for each sample. Total RNA samples were purified by electrophoretic separation on a 15% urea denaturing polyacrylamide gel electrophoresis (PAGE) gel and small RNA regions corresponding to the 18-30 nt bands in the marker lane (14-30 ssRNA Ladder Marker, TAKARA) were excised and recovered. Then the 18-30 nt small RNAs were ligated to adenylated 3′ adapters annealed to unique molecular identifiers (UMI), followed by the ligation of 5′ adapters. The adapter-ligated small RNAs were subsequently transcribed into cDNA by SuperScript II Reverse Transcriptase (Invitrogen, USA) and then several rounds of PCR amplification with PCR Primer Cocktail and PCR Mix were performed to enrich the cDNA fragments. The PCR product for cDNA at 45DAF were selected by agarose gel electrophoresis with target fragments 110~130 bp, and then purified by QIAquick Gel Extraction Kit (QIAGEN, Valencia, CA). The libraries were quality and quantitated in two methods: check the distribution of the fragments size using the Agilent 2100 bioanalyzer, and quantify the library using real-time quantitative PCR (qPCR) (TaqMan Probe). The final ligation PCR products were sequenced using the BGISEQ-500 platform (BGI-Shenzhen, China). On the other hand, four cDNA libraries at 5, 20,45,90 DAF for miRNAs were used as templates to detect miRNA expression levels (see qRT-PCR method section).
After obtaining the sequencing results, further trimming and filtering the adaptor and low-quality tag sequences, and the high-quality sRNA clean reads were mapped into the Rfam (https://rfam.xfam.org) to filter the rRNA, tRNA, snRNA, and snoRNA. The filtered reads were then compared against known plant miRNA database in the miRBase 21.0 (http://www.mirbase.org/) with BLASTn. After searching against Rfam database and miRBase, the remaining reads were further mapped to the grape reference genome (http://genomes.cribi. unipd.it/DATA/V2/).

Bioinformatics analysis and identification of VvmiRNA and VvmiRNA SNPs
The clean reads were screened from raw data by filtering out the corrupted adapter sequences, poly-A tails and sequences with ≤18 nt and ≥ 30 nt. The clean read sequences were mapped into the Rfam (https://rfam.xfam. org) to filter the rRNA, tRNA, snRNA and snoRNA etc. The filtered reads were then compared against known plant miRNA database existing in the miRBase 21.0 with BLASTn. Only matching (0-3 mismatches) sequences in their sequences' ends were considered as known VvmiR-NAs, while other sequences have one base variation with the known VvmiRNAs in the middle sites of their sequences and thus can be considered as miRNA SNV. On the other hand, the identification criteria of novel miRNA as follow: 1) except for the identified known VvmiRNAs and VvmiRNA SNVs, the remaining sequences were mapped to the grape reference genome (http://genomes.cribi.unipd.it/DATA/V2/) and mRNA sequence. The reads mapped to genome but mRNA were used to predict the potential miRNA precursor with mireap, and then these reads werre processed by miRCat (http://srna-tools.cmp.uea.ac.uk) [70] using default parameters to generate the secondary structures; 2) The negative free energy of folding structure were less than -20kj; 3) The both arms of stem-loop structures contained the bubbles with less than 6 mismatched bases; 4) The first base of miRNAs possessed the "U" preference; 5) The length of miRNA is usually in the range of 19 -24 nt. In addition, as to the depth coverage and frequency filters for reliable calling of SNVs on miR-NAs, here the depth coverage was required to be more than 4 (> 4), and the frequency filters was more than 0.05 (> 0.05).

Identification of precise sequences of VvmiRNAs by miR-RACE
The cDNA was amplified with miR-Racer 5′ primer, 3′ primer and their corresponding gene-specific primers to generate 5′-and 3′-miR-RACE products, respectively, for the identification of precise sequences of miRNAs [16,69]. The clone products of 5′-and 3′-miR-RACE were approximately 56 and 87 bp, respectively. 3′-miR-RACE was performed using common primer 1(CP1) (ATTCTAGAGGCCGAGGCGGCCGACATG) and miRNA specific primer 1(MGSP1), while 5′-miR-RACE was performed using CP2 (GGAGCACGAGGACACTGA CATGGACT) and MGSP2. The design procedure of MGSP1 and MGSP2 primers as follow: MGSP1 primer consists of 10 bp adaptor sequence (GGAGTAGAAA) add 17 bp sequence intercepted from 5'end of miRNA, while MGSP2 primer includes 10 bp Poly(T) and 17 bp complementary sequence cut off from 3'end of miRNA [16,69]. Here, the 17 nucleotides complimentary to the miRNA were sufficient for the accurate and efficient PCR amplification of the opposite ends. The primer specificity was validated by inspecting the specific band of PCR product. All specific primers are listed in Supplementary Table S7.

Validation of potential target genes for VvmiRNAs with RLM-RACE and PPM-RACE
The mRNA library was ligated with 5′-adapter or 3′-PolyA tail [29] and then reverse transcribed as cDNA. RLM-RACE and PPM-RACE [28] were performed with their corresponding cDNA and primers, respectively (Supplementary Table S6). The products of RLM-or PPM-RACE were mapped into the target genes for validation of potential target genes and identification of cleavage sites and frequency. RLM-RACE was carried out using the common primer 1 (GGAGCACGAGGACACTGACA TGGACT) and the specific primers (P1); and PPM-RACE was performed with the common primer 2 (ATTCTAGA GGCCGAGGCGGCCGACATG) and the specific primers (P2). P1 is the reverse primer at the downstream of the predicted cleavage site in target gene, while P2 is the forward primer at the upstream of corresponding cleavage site in target gene. The primer specificity was validated by inspecting the specific band of PCR product. P1 and P2 were listed in Supplementary Table S8.

qRT-PCR analysis of VvmiRNAs and their target genes
For qRT-PCR expression analysis of VvmiRNAs and their target genes, the cDNA libraries for mRNA and miRNA from control and GA treatment at diverse stages were constructed using our developed methods [13], and then these cDNAs were used as templates to detect their corresponding expression levels using their primers with three replicates (Table 3), and the expression levels were normalized using 5.8S rRNA. The relative expression levels of the miR-NAs and their targets were calculated by using the 2 −ΔΔCT method. miRNA qRT-PCR was amplified using common primer (qP1: ATTCTAGAGGCCGAGGCGGCCGAC ATG) and specific primer (qP2: miRNA sequence). The primer specificity was validated by inspecting the specific band of PCR product. The U6 gene was used as the reference gene for the normalization of all miRNAs' relative expression values, and the actin gene was used as the referenced one for that of all target genes' relative expression values. All primers were listed in Supplementary Table S9.  VvmiR53-3p GGCAGCAGCAUACUACUUUG GGCAGCAGCAUACUACUUUG Yes the experimental cost and publication fee for this work. However, the funding agencies had no participation in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Availability of data and materials
The RNA-seq data have been deposited into the NCBI GEO under the accession number GSE182618 (https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE182618). All data generated or analyzed during this study are included in this published article and its additional files.