Uncovering anthocyanin biosynthesis related microRNAs and their target genes by small RNA and degradome sequencing in tuberous roots of sweetpotato

Background Compared with white-fleshed sweetpotato (WFSP), purple-fleshed sweetpotato (PFSP) is a desirable resource for functional food development because of the abundant anthocyanin accumulation in its tuberous roots. Some studies have shown that the expression regulation mediated by miRNA plays an important role in anthocyanin biosynthesis in plants. However, few miRNAs and their corresponding functions related to anthocyanin biosynthesis in tuberous roots of sweetpotato have been known. Results In this study, small RNA (sRNA) and degradome libraries from the tuberous roots of WFSP (Xushu-18) and PFSP (Xuzishu-3) were constructed, respectively. Totally, 191 known and 33 novel miRNAs were identified by sRNA sequencing, and 180 target genes cleaved by 115 known ib-miRNAs and 5 novel ib-miRNAs were identified by degradome sequencing. Of these, 121 miRNAs were differently expressed between Xushu-18 and Xuzishu-3. Integrated analysis of sRNA, degradome sequencing, GO, KEGG and qRT-PCR revealed that 26 differentially expressed miRNAs and 36 corresponding targets were potentially involved in the anthocyanin biosynthesis. Of which, an inverse correlation between the expression of ib-miR156 and its target ibSPL in WFSP and PFSP was revealed by both qRT-PCR and sRNA sequencing. Subsequently, ib-miR156 was over-expressed in Arabidopsis. Interestingly, the ib-miR156 over-expressing plants showed suppressed abundance of SPL and a purplish phenotype. Concomitantly, upregulated expression of four anthocyanin pathway genes was detected in transgenic Arabidopsis plants. Finally, a putative ib-miRNA-target model involved in anthocyanin biosynthesis in sweetpotato was proposed. Conclusions The results represented a comprehensive expression profiling of miRNAs related to anthocyanin accumulation in sweetpotato and provided important clues for understanding the regulatory network of anthocyanin biosynthesis mediated by miRNA in tuberous crops. Electronic supplementary material The online version of this article (10.1186/s12870-019-1790-2) contains supplementary material, which is available to authorized users.


Background
Sweetpotato (Ipomoea batatas L.), a hexaploid (2n = 6x = 90) dicotyledonous plant of Convolvulaceae family, is an important crop around the world due to its high yield, wide adaptability and rich nutrition [1]. The flesh of the tuberous roots has multiple colors, such as white, yellow, orange and purple [2]. The purple-fleshed sweetpotato (PFSP) is not only as nutritious as the white-fleshed sweetpotato (WFSP) but also enriches high content of anthocyanin. Anthocyanin is water-soluble pigment belonging to the flavonoid group. Numerous studies have showed that anthocyanin has a strong antioxidant activity and therefore has therapeutic effects on a variety of diseases like obesity and cancers [3]. More importantly, the anthocyanin of PFSP has high thermal stability and light stability with a high level of acylation than that of strawberry, raspberry, and apple [4][5][6]. Therefore, the PFSP has been recently proposed as a potential pharmaceutical crop for developing drugs, such as antineoplastic, antiinflammatory, and antioxidant agents [6,7].
In this study, sRNA and degradome sequencing were used to identify miRNAs and their corresponding targets that potentially involved in anthocyanin biosynthesis in sweetpotato. Small RNA and degradome libraries from the tuberous roots of WFSP (cultivar Xushu-18, XS-18) and PFSP (cultivar Xuzishu-3, XZS-3) were constructed and sequenced. Integrated analysis of sRNA, degradome sequencing, GO, KEGG and qRT-PCR revealed a comprehensive account of the ib-miRNA populations, corresponding targets, expression abundance, as well as miRNAs potentially involved in the anthocyanin metabolism. Expression levels of eight differently expressed miR-NAs and their targets were validated by qRT-PCR, which were consistent with the sequencing data. In addition, we demonstrated that the over-expression of ib-miR156 in Arabidopsis strongly suppressed the abundance of ibSPL and resulted in a purplish phenotype. Finally, a possible ib-miRNA-target regulatory model associated with anthocyanin biosynthesis in tuberous roots of sweetpotato was illustrated. Our findings provided a comprehensive expression profiling of ib-miRNAs, and suggested that miR-NAs were involved the regulation of anthocyanin biosynthesis in sweet potato.

Result
Small RNA populations in the tuberous roots of sweetpotato To identify miRNAs involved in the process of anthocyanin biosynthesis, the tuberous roots of WFSP (XS-18) and PFSP (XZS-3) were used to construct sRNA libraries and sequenced on Illumina HiSeq2000 platform, respectively. A sum of 27,705,914 (XS-18) and 28,947,914 (XZS-3) raw reads were produced. After removing of the low quality reads and contaminated adapter sequences, 26,825,634 (96.82%) and 27,331,707 (94.42%) clean reads were obtained from XS-18 and XZS-3, respectively. A total of 15,764,489 (XZS-3, 91.61%) and 13,904,486 (XS-18, 73.32%) sRNAs were mapped to the corresponding assembled unigenes of sweetpotato (Additional file 1). The mapped reads were further annotated against Pfam database, and subsequently divided into rRNAs, tRNAs, snRNAs, snoRNAs, ta-siRNA (TAS) and others. The endogenous sRNAs were identified as known miRNAs and novel miRNAs (Additional file 2). The length distribution patterns of the sRNAs were similar in the two libraries. They were ranged from 18 to 30 nt, of which 24 nt sRNAs were the most abundant size with 18. 63 and 15.19% in XS-18 and XZS-3 libraries, respectively (Fig. 1).

Identification of known and novel ib-miRNAs in sweetpotato
The known miRNAs were found to be important in plant growth, development and many other biological processes [38]. To identify known miRNAs in sweetpotato, the sRNA sequences were aligned to the mature miRNA sequences deposited in the miRBase 22.1. A total of 191 known miRNAs belonging to 43 miRNA families were identified in the two libraries, with 185 and 145 known miRNAs in XS-18 and XZS-3, respectively (Additional file 3). Among the 43 identified families, the member number of each miRNA family varied significantly. The ib-miR396 family contained the largest number with 18 members, followed by ib-miR166 and ib-miR159 families, with 17 and 15 members, respectively. By contrast, 9 families had only one member, such as ib-miR157, ib-miR394, ib-miR5083 and ib-miR5658 (Fig. 2).
The abundance of the miRNAs were estimated as transcripts per million (TPM), which highly varied among the 43 known miRNA families (Additional file 4). Some miRNAs highly expressed with more than 1000 TPM, such as ib-miR159a-1 and ib-miR319a-1, whereas some miRNAs expressed at a lower level with less than 2 TPM, such as ib-miR160-1, ib-miR1030a-j and ib-miR1128. In addition, different members in the same miRNA family exhibited different expression levels. For instance, in miR159 family, ib-miR159a-1 showed the highest abundance with 5978 and 1106 TPM in the XS-18 and XZS-3 libraries, respectively. However, ib-miR159f only presented 4 and 1 TPM in the XS-18 and XZS-3 libraries, respectively (Additional file 4).
The length of known miRNAs ranged from 18 nt to 24 nt, in which 21 nt miRNAs were the most abundant size (Additional file 3). To understand the base preference of sweetpotato miRNAs, we analyzed the base distribution for each position of known miRNAs. In both libraries, nucleotide sequences analysis revealed that uridine (U) was the most common nucleotide at the 5′ end. The cleavage site of the target gene usually matched with the 10 th or 11 th nucleotide of miRNA. The majority nucleotides of known ib-miRNAs was adenine (A) at position 10 th and thymine (T) or A at position 11 th (Fig. 3a).
The remaining sRNA sequences were mapped to sweetpotato transcriptome database, and the hairpin structures were used to identify novel miRNAs in sweetpotato. Novel miRNAs were predicted based on the published criteria [39]. Totally, 33 and 33 novel miRNAs were predicted in the XS-18 and XZS-3 libraries, respectively (Additional file 5). The length distribution of the novel miRNAs was between 20 nt to 25 nt. Similar to known miRNAs, the nucleotide bias analysis showed that U appeared mainly at the 5′ end of novel miRNAs (Fig. 3b). The length of the novel miRNA precursors ranged from 57 nt to 293 nt, with an average length of 140 nt (Additional file 5), which was in accordance with   star (miR*) strands lays strong support for the authenticity of novel miRNA [39]. Similar to the known miR-NAs, the abundance of the identified novel ib-miRNAs in sweetpotato varied considerably. MiR* strands were found to be present in 17 of the 33 novel ib-miRNAs and expressed at lower levels than their corresponding miR strands (Additional file 5), which agreed with the report that miR* strands were mostly degraded [40].

Target validation of ib-miRNAs by degradome sequencing
To explore the function of sweetpotato miRNAs, computational program psRNATarget (http://plantgrn.noble. org/v1_psRNATarget/) was performed to predict their target genes. A total of 3278 transcripts were predicted to be targets of the 224 miRNAs (Additional file 7). Then, high-throughput degradome sequencing was performed to validate the miRNA targets. A total of 1,886,192 unique reads of 14.51 million raw reads were obtained, and 766,328 unique reads matched with the reference sequences assembled from sweetpotato transcriptome database (Additional file 8). Totally, 180 target genes cleaved by 115 known ib-miRNAs and 5 novel ib-miRNAs were identified by CleaveLand4 analysis (http://sites.psu.edu/axtell/software/cleaveland4/) (Additional file 9). The cleaved-target transcripts were classified into five categories, named category 0, 1, 2, 3 and 4, according to the relative abundance of the tags at the target sites. The miRNAs and corresponding targets in the five categories were shown in Additional file 9. Among the 180 identified targets, 55, 8, 34, 1 and 23 targets were found in categories 0, 1, 2, 3 and 4, respectively.

Validation of ib-miRNAs and their target genes by qRT-PCR
Eight differentially expressed miRNAs and their corresponding targets between XS-18 and XZS-3 were randomly selected for qRT-PCR to verify the results obtained from sRNA and degradome sequencing. In general, the expression patterns of the eight selected miRNAs showed by qRT-PCR were consistent with that showed by high-throughput sequencing, indicating that the sRNA sequencing data were reliable (Fig. 6). Up-regulated expression levels of ib-miR858b, ib-miR2111a-5p and ib-miR156a-5p were detected in XZS-3 compared to XS-18, whereas the abundance of their corresponding target genes MYB, ERF and SPL were down-regulated, respectively ( Fig. 6a, b, c). The expression amounts of ib-miR160e-5p and ib-miR166m were decreased, while their corresponding target genes ARF and HD-ZIP had an increased expression in XZS-3 compared to XS-18 ( Fig.  6d, e). The inversely correlated expression patterns between miRNA and their corresponding targets confirmed degradome sequencing data. However, the expression of ib-miR396g-5p, novel_miR-1 and novel_-miR-5 showed the same trend with that of their corresponding transcripts, Cluster-18,233.16052 (WD40), Cluster-18,233.73593 (MYB) and Cluster-18,233.59778 (ARF), respectively, suggesting that the transcripts may not be the cleaved targets of these miRNAs (Fig. 6f, g, h). The expressions of another five miRNAs listed in Table 1 were also verified by qRT-PCR (Additional file 11).

GO and KEGG analysis of targets regulated by differentially expressed ib-miRNAs
A total of 3278 transcripts were predicted to be targets of the 224 miRNAs by computational program psRNA-Target. Then, 180 targets identified by degradome sequencing were subjected to gene ontology (GO) analysis, and 124 of the 180 targets were found in GO database, which were involved in biological processes (79), molecular function (83), and cellular components (71) (Additional file 12). For biological processes category, more than 40% of the target genes were involved in regulation of transcription, DNA-templated annotations (GO:0006355), followed by DNA-templated transcription, initiation (GO:0006352) and DNA replication (GO:0006260) (Fig. 7a). In the molecular function category, the three most dominant terms were DNA binding (GO:0003677), sequence-specific DNA binding transcription factor activity (GO:0003700), and protein binding (GO:0005515) (Fig. 7b). Within the cell component category, the dominant two enrichments were the nucleus (GO:0005634) and the transcription factor complex (GO:0005667) (Fig. 7c).

Identification of ib-miRNAs and their targets related to anthocyanin biosynthesis in sweetpotato
Identification of differentially expressed miRNAs in tuberous roots of WFSP and PFSP could help to better understand the biological function of miRNAs involved in anthocyanin biosynthesis. Combined with sRNA, degradome sequencing, GO and KEGG analysis, 26 of the 121 miRNAs differentially expressed miRNAs were identified to be potentially involved in anthocyanin biosynthesis in sweetpotato (Table 1, Additional file 14).

Over-expression of ib-pri-MIR156 induced anthocyanin accumulation in Arabidopsis
The orthologous gene of miR156 in Arabidopsis has been demonstrated to regulate anthocyanin biosynthesis [14]. In this study, the inverse correlation between the accumulation of ib-miR156 and ibSPL, the targets of ib-miR156, in WFSP and PFSP suggested that ib-miR156 may also participate in anthocyanin biosynthesis. To test this hypothesis, the recombinant vector, pc2300-pO-T2-ib-pri-MIR156, was generated and transformed into Arabidopsis. Interestingly, a purplish phenotype was observed in young seedlings of the ib-pri-MIR156 over-expressing Arabidopsis under normal growth conditions (Fig. 10). Analysis by qRT-PCR showed that the abundance of miR156 was greatly increased, whereas the expression of SPL decreased significantly in the transgenic Arabidopsis (Fig. 10e). Besides, upregulated expressions were detected for four structural genes related to anthocyanin biosynthesis, including CHS, CHI, DFR and ANS in the transgenic Arabidopsis (Fig. 10e). The results indicated that ib-pri-MIR156 positively regulated anthocyanin accumulation by repressing the expression of SPL.

Identification of ib-miRNAs by high-throughput sequencing
MiRNAs, as the key post-transcriptional regulators, participate in various biological processes in plant. Recently, increasing evidences have showed that plant miRNAs are also involved in secondary metabolism [14,38,42]. Anthocyanin is one of the important secondary metabolite products. A number of miRNAs have been reported to play important roles in regulating anthocyanin biosynthesis in plant, such as miR156 [14], miR393 [43], miR828 [18,25,26], miR163 [44], miR165/166 [17] and miR858 [45]. However, the identification of miRNAs related to anthocyanin biosynthesis has not yet been reported in the tuberous root of sweetpotato.
In this study, a total of 191 known ib-miRNAs were identified from the WFSP and PFSP libraries by analysis of sRNA sequencing data. The number was much higher than previous reports on sweetpotato. In the first report on sweetpotato miRNAs, only 8 potential miRNA candidates were computationally identified [46]. Then, 24 known miRNAs were identified in Solexa Sequencing of Ningzishu 1 [36]. In the two studies, the number of miRNA identified is considerably limited because EST sequences were used as the references. The known miR-NAs number of this study is comparable with the recent report on miRNAs responded to chilling stress of sweetpotato, in which 190 known miRNAs were found [47]. The research is based on Kazusa sweetpotato GARDEN database, which scaffold the genome sequences of a wild diploid ancestor of sweetpotato, Ipomoea trifida, although the total lengths of the assembled sequences was only 712 Mb and the size is far less than estimated genome sequence [48].
Most of the known miRNA families detected by Xie et al were also present in our dataset, except miR1508, miR5253, miR5298 and miR2911 [47]. On the other hand, we identified miR6478, miR6300, miR6173 and miR5508 which were not reported in sweetpotato previously. Since the expression levels of miRNAs existed temporal and spatial specificity, this result is not surprising [49]. Based on the hairpin structures of pre-miRNAs (Additional file 6), 33 novel candidate miRNAs were identified in this study. The number of novel miRNAs is clearly more than the research with the EST sequences as reference, but it is far less than the report of miRNAs responded to chilling stress of sweetpotato. Although the genomic sequence of Taizhong 6 sweetpotato cultivar was released last year, the incompleteness and inaccuracy of assembly limited its application as a good reference genome of hexaploid sweetpotato [50,51]. In order to eliminate the difference caused by genetic background, species and development stages, the transcriptome sequences developed with the same sweetpotato cultivars with sRNA deep sequencing were used as the reference to predict the novel miRNAs in this study. Degradome sequencing has been successfully applied to identify miRNA targets in many plant species [52,53]. In this study, a total 180 targets for 115 ib-miRNAs were obtained by the degradome sequencing approach. The number of targets supported by degradome was far less than that predicted by using computational psRNATarget. As expected, a majority target genes were transcription factors in plant species, and some ib-miRNA had more than one targets. Such as, ib-miR156a-5p and ib-miR5658 both had 10 targets and ib-miR159a had 6 targets. In addition, some ib-miRNAs belonging to the same family shared the same transcript. According to the annotation of targets, these target genes were found to participate in wide biological processes.

Ib-miRNAs involved in anthocyanin biosynthesis in sweetpotato
Auxin is essential for plant development. High levels of auxin can repress the expression of the MBW complex and thereby regulate the biosynthesis of anthocyanin [54,55]. ARF, WRKY and HD-ZIP genes play important roles in auxin-mediated signaling, which regulates anthocyanin biosynthesis through the Aux/IAA-ARF pathway in apple [56]. In plants, ARF genes were regulated by miR160 and miR390 [24,57]. The WRKY mainly involved in auxin-mediated signaling and had considerable effect on flavonoid and anthocyanin biosynthesis [58,59]. HD-ZIP TF plays critical roles in shoot apical meristem and organ polarity in plant [23]. Our previous study showed that blockage of miR165/166 caused the up-regulation of HD-ZIP TFs and increased IAA content accompanied by enhanced anthocyanin content in Arabidopsis [17]. In this study, the analysis of sRNA and degradome sequencing demonstrated that WRKY and ARF were targeted by ib-miR172e-3p and ib-miR160e-5p, respectively; while HD-ZIPs were targeted by ib-miR164c and ib-miR166m (Table 1). These results suggested that ib-miR172e-3p, ib-miR160e-5p, ib-miR164c and ib-miR166m might regulate anthocyanin biosynthesis through the auxin signaling.
CAD is one of the key enzyme genes of lignin biosynthesis [60], which has been found to be targeted by ib-miR156, ib-miR157 and ib-miR396 (Fig. 9). MiR156 and miR159 have been reported as potential regulators of secondary wall biosynthesis in Acacia mangium [61]. In this study, four members of ib-miR159 family were found to have a common novel target I2'H ( Fig. 9), which can promote the accumulation of medicarpin [62]. MiR6300, only detected in few leguminous plants like Catharanthus roseus and Camellia sinensis, was found to target UFGT [63,64]. Eight miRNAs were found to regulate the enzyme genes involved in the phenylpropanoid pathway in sweetpotato in this study. Based on the interactions of the miRNAs and enzyme (See figure on previous page.) Fig. 6 Comparison of the expression levels of the miRNAs and their corresponding target genes determined by qRT-PCR and deep sequencing. a-h The Y-axis on the left and right represented the relative expression levels determined by qRT-PCR and deep sequencing, respectively. Ib5S rRNA and ibActin were used as an internal control for miRNAs and targets, respectively. The expression levels of the miRNAs and targets in XS-18 were set as 1.0. Relative expression was calculated using the 2 -ΔΔCt method. The bars represented the mean ± SD values of three biological experiments. * and ** indicated a statistically significant difference between XS-18 and XZS-3 at P < 0.05 and 0.01 according to DMRT analysis, respectively genes in phenylpropanoid biosynthesis, we proposed a pathway panel for phenylpropanoid regulation in sweetpotato (Fig. 9).
The MYB TFs have been identified to be one of the major regulators in the pathway of anthocyanin biosynthesis [65][66][67]. For example, AtPAP1 (AtMYB75) and its orthologs have been shown to effectively induce anthocyanin production in various plant species [68][69][70].
Over-expression of MYB1 induced anthocyanin accumulation in sweetpotato and apple [65,71]. MYB has been demonstrated to be regulated by many miRNA members. For instance, miR828 was involved in the anthocyanin production, trichome and cotton fiber development by regulating MYB in plant [16,18,72,73]. MiR858 was found to play a positive role in the accumulation of anthocyanin by cleaving MYBL2 in Arabidopsis [9]. MiR858 has also been identified to target up to 66 MYB members in apple [74,75]. Similarly, our recent studies showed that Sl-miR858 regulated two SlMYB transcripts and functioned in anthocyanin accumulation in tomato [45]. In this study, MYB genes were also targeted by ib-miR858, ib-miR319, ib-miR159 and ib-miR156. Compared with WFSP, ib-miR858b was significantly up-regulated whereas its corresponding target gene ibMYB down-regulated in PFSP (Table 1; Fig. 6). The WD40 protein was reported to be essential and irreplaceable in the MYB-bHLH complex for anthocyanin biosynthesis in sweetpotato [76]. The transcript of ibWD40 was not  cleaved by ib-miR396g-5p in degradome sequencing, although it was a potential target of ib-miR396g-5p predicted by psRNATarget (Additional file 7). Ib-miR396g-5p and ibWD40 both expressed higher in PFSP than in WFSP showed by qRT-PCR (Fig. 6). One possible reason is that the miRNAs may impose a translational regulation on their targets, which could not be detectable by degradome sequencing. Transcription factor bHLH is essential for the activity of the R2R3-MYB partner and can directly bind to the promoters of DFR and UFGT to activate or inhibit anthocyanin biosynthesis [77,78]. However, none of miRNA was identified to target bHLH in this study. The enzyme gene LcUFGT plays a major role in anthocyanin accumulation [79]. In sweetpotato, we found ib-miR6300 directly targeted ibUFGT, suggesting that b-miR6300 may participate in anthocyanin biosynthesis by regulating ibUFGT (Additional file 9). However, the regulatory mechanism needs further validation. SPLs are transcription factors widely existing in plants, which play an important role in plant growth and development, primary and secondary metabolism, as well as other biological processes. In Arabidopsis, SPL was demonstrated to be regulated by miR156, and acted as negative regulators of anthocyanin accumulation by destabilization of the MBW complex in Arabidopsis [14]. In this study, ib-miR156a-5p was identified to cleave ibSPL genes (Additional file 9). The results of sRNA sequencing and qRT-PCR showed that miR156a-5p significantly up-regulated whereas its target ibSPL expressed downregulated in PFSP compared with WFSP (Fig. 6). By over-expressing ib-miR156 in Arabidopsis, a purple phenotype with high anthocyanin accumulation in the main stem of transgenic Arabidopsis was observed. The results of qRT-PCR showed that the over-expression of ib-miR156 in Arabidopsis strongly suppressed the abundance of SPL (Fig. 10), which consistent with the results of Gou et al. [14]. In addition, up-regulated transcripts for four structural genes of anthocyanin pathway were detected in ib-miR156 over-expressing transgenic plant (Fig. 10), suggesting that ib-miR156 may modulate anthocyanin biosynthesis through regulation of the structure genes in the phenylprepanoid pathway. These results indicated that the regulatory function of ib-miR156 is similar to that of at-miR156.
Putative ib-miRNA-target model involved in anthocyanin biosynthesis in sweetpotato Anthocyanin biosynthesis was modulated by regulating genes, including MYB, WDR, bHLH, SPL, ARF and MADS-box. In addition, sucrose synthase, ABC transporter and sugar/inositol transporter were also potentially participated in anthocyanin biosynthesis. In this study, according to the annotation of targets, we found that 26 differentially expressed miRNAs and 36 their corresponding regulator genes were more likely to participate in anthocyanin biosynthesis in sweetpotato (Table 1, Additional file 14). Based on the present data and previous reports, a possible ib-miRNA-target model related to anthocyanin biosynthesis was proposed. As shown in Fig. 11, miRNA families of ib-miR156, ib-miR159, ib-miR396 and ib-miR858 potentially targeted SPL, MYB and WDR to shape a regulatory MBW complex; miRNA families of ib-miR164, ib-miR165/166 and ib-miR172e-3p.2 possibly targeted WRKY and HD-ZIP to involve in auxin signaling; ib-miR160, miRNA families of ib-miR172 and ib-miR2111 might regulate ARF and ERF to activate auxin-mediated signaling and sucrose signaling. All of these biological processes would influence the expression of structural genes, and subsequently form a complex regulatory network to modulate anthocyanin biosynthesis.

Conclusions
In summary, this is the first report on systematic identification, expression analysis and potential roles of miR-NAs and their targets in regulating anthocyanin biosynthesis in tuberous roots of sweetpotato. A total of 26 differentially expressed miRNAs and 36 corresponding targets were more likely to be related to anthocyanin biosynthesis by bioinformatic analysis and expression validation. By Agrobacterium mediated genetic transformation, ib-miR156 over-expressing transgenic Arabidopsis with purplish phenotype were obtained. The expression amount of ib-miR156 was up-regulated in the transgenic lines, while its target SPL showed down-regulated expression. Four anthocyanin-specific enzyme genes, CHS, CHI, DFR and ANS, expressed significantly higher in transgenic Arabidopsis than in the wild type plants (WT), suggesting that ib-miR156 could positively mediate anthocyanin biosynthesis by modulating related structural genes. Base on above results, a putative ib-miRNA-target model associated with anthocyanin biosynthesis in sweetpotato was proposed. Our findings provided comprehensive information for anthocyanin-specific miRNAs and their targets, as well as a starting point for mechanism investigation of miRNAs in anthocyanin biosynthesis in sweetpotato.

Plant materials
Sweetpotato (Ipomoea batatas L.) cultivars of white-fleshed Xushu-18 (XS-18) and purple-fleshed Xuzishu-3 (XZS-3) were used for sRNA and degradome sequencing. The cultivars were provided by the Sweetpotato Research Institute in Jiangsu province of China and grown in the experimental field of the Shanxi Agricultural University (Taigu, China) under the normal conditions. After root formation at the stage of 90 days, white and purple tuberous roots with similar size and shape were selected. After gently washed, the roots were immediately frozen by liquid nitrogen and stored at − 80°C freezer. We mixed three tuberous roots samples together as one biological replicate to eliminate effects of individual genetic variance, respectively. In total, three independent biological replicates were collected for each sample.

Small RNA library construction and sequencing
Total RNAs were extracted from sweetpotato roots using Trizol reagent (Invitrogen, USA) according to the manufacturer's instructions. Total RNA from each genotype was isolated in three biological replicates and the Identification and expression analysis of sweetpotato miRNAs using deep sequencing The process for known and novel miRNAs in sweetpotato tuberous roots was performed following the previous reported approach [39]. Firstly, clean reads were obtained from raw data by removing reads containing ploy-N, larger with 5'adapter contaminants, without 3'adapter or the insert tag, containing ploy A or T or G or C and low quality reads. Sequences smaller than 18 nucleotide (nt) and larger than 30 nt were also removed. Then the sRNA tags were mapped to reference sequence by Bowtie [80] without mismatches allowed. The transcriptome sequences using the same sweetpotato cultivars (XS-18 and XZS-3) with sRNA deep sequencing were developed in our lab (not published). The assembled sequences from the transcriptome database were used as reference sequence for predicting miRNA precursors. The mapped sequences were subjected to BLASTn search against Repeat Masker (http://www.repeatmasker. org/) and Rfam database (http://rfam.xfam.org/) to remove tags originating from protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, or those types of data from the specified species itself. Usually, we followed the following priority rule: known miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > gene > NAT-siRNA > gene > novel miRNA > ta-siRNA. The remaining reads were aligned with the miRNA sequences deposited in the miRBase 22.1 database (http://www.mirbase.org/) [81]. The matched sequences with no mismatches were considered to be known miRNAs. The unaligned reads were then subjected to software miREvo and miRdeep2 to predict novel miRNA candidates through exploring the secondary structure, the Dicer cleavage site and the minimum free energy [82]. The stem loop hairpin structures of pre-miRNAs were also predicted using RNAfold software. MiRNA expression levels were estimated by TPM (transcript per million) through the following Normalization formula: Normalized expression = mapped read count/total reads*1000000 [83].
To reveal the differentially expressed miRNAs related to anthocyanin biosynthesis between XS-18 and XZS-3, the miRNAs expression was analyzed using the DEGseq R package, respectively. P-value was adjusted using qvalue [84]. Qvalue < 0.01 and |log2 (fold change)| > 1 was set as the threshold for significantly differential expression by default.

Degradome library construction, sequencing and data analysis
To identify potential target mRNAs for sweetpotato miR-NAs, a degramdome library was constructed using the mixed roots of XS-18 and XZS-3. Firstly, by using the Oligotex kit (Qiagen, Germany), 200 μg of total RNA was used for extracting poly (A) RNA, which was ligated to a 5′ adapter with an EcoP15 I recognition site in its 3′ end. After ligated, the first-strand cDNA was generated and amplify by PCR. The PCR product was purified and digested with PAGE-gel and EcoP15 I, respectively. Then, the EcoP15 I cleaved fragments were ligated to a 3′ double-strand DNA adapter and followed by PAGE-gel purification to obtain the ligated products. After PCR amplification, PAGE-gel was used for the third time to purify the final products. Finally, the purified cDNA library was ready for deep sequencing on Illumina HiSeq2000 sequencing system (LC-BIO Sciences, China).
After sequencing, the adapter sequences, low-quality reads and N-containing fragments were filtered from the raw reads. The remaining sequences mapped to the sweetpotato transcriptome database of XS-18 and XZS-3 were used to identify potentially cleaved targets by CleaveLand4 pipeline as previously described [85]. The degradome event mediated by sRNA was categorized as 0, 1, 2, 3, and 4 based on the abundance of its cleaved tag reads across the target transcript [86]. Category 0: Abundance equal to the maximum of the target transcript abundance and the cleaved tag reads had only one maximum value. Category 1: Abundance equal to the maximum of the target transcript abundance, and the cleaved tag reads had more than one maximum value. Category 2 represented the cleaved tag abundance was less than the maximum but higher than the mean of the transcript abundance. Category 3: the cleaved tag abundance was less than or equal to the mean of the transcript abundance. Category 4 represented the depth of the position equal to 1. T-plots were built to analyze the miRNA targets and RNA degradation patterns according to the distribution and abundances along these transcripts. The potential targets of miRNAs were analyzed by PAREsnip software with P-value < 0.05 [87].

GO enrichment and KEGG pathway analysis
GO terms of the miRNAs targets were annotated according to their biological role, molecular function and cellular component by using the online GO analysis tool (http://www.geneontology.org/page/go-enrichment-analysis) [88]. The statistical enrichment of the candidate target genes in KEGG pathways was set by KOBAS software [89,90].

Construction of ib-pri-MIR156 over-expressing vector and genetic transformation of Arabidopsis
Ib-pri-MIR156 was isolated from PFSP, and inserted into a pOT2 vector between the 35S2 promoter and the 35S terminator. Then, the pOT2 was mobilized into a modified pCAMBIA2300 binary vector to construct over-expressing vector pc2300-pOT2-ib-pri-MIR156. The vector was introduced into Agrobacterium GV3101 and followed by Agrobacterium mediated transformation to obtain transgenic Arabidopsis [91]. All the primers used for plasmid construction were listed in Additional file 15.
Analysis of the expressions of ib-miRNAs, corresponding targets and anthocyanin-specific enzyme genes by qRT-PCR The expressions of miRNAs and their corresponding targets were evaluated by qRT-PCR. Total RNA was extracted from XS-18 and XZS-3 by using Trizol (Invitrogen, USA). For miRNA expression analysis, the specific stem loop RT, forward and universal reverse primers were designed according to the reported method with some modifications [92]. One μg of DNase treated total RNA was used for cDNA synthesis following the procedures of PrimeScript™ 1st Strand cDNA synthesis kit (TaKaRa, China). All qRT-PCR analysis was carried out using SYBR Premix Ex Taq™ (Tli RNaseH Plus) (TaKaRa, China) and performed on a Bio-Rad CFX96 Real-Time PCR Systems (Bio-Rad, USA). For sweetpotato, the ib5S and ibActin gene were used as internal reference for accurate normalization in each reaction for the miRNA, and target genes, respectively [93]. For transgenic Arabidopsis over-expressing miR156, the atU6 and atActin were used. The relative expression levels were calculated using the 2 -ΔΔCt method [94]. Each sample from three biological replicates was performed in triplicate. Results were presented as means ± SD. The method of Dunnett's two-tailed t-test was used to conduct the statistical analysis of RT-PCR results, and the statistical significant differences were shown at p ≤ 0.05 (marked *) and p ≤ 0.01 (marked **). The sequences for all the primers were listed in Additional file 15.

Availability of data and materials
The datasets generated and/or analyzed during the current study are included in Additional Tables.