Skip to main content

Integrated microRNA and transcriptome profiling reveal key miRNA-mRNA interaction pairs associated with seed development in Tartary buckwheat (Fagopyrum tataricum)

Abstract

Background

Tartary buckwheat seed development is an extremely complex process involving many gene regulatory pathways. MicroRNAs (miRNAs) have been identified as the important negative regulators of gene expression and performed crucial regulatory roles in various plant biological processes. However, whether miRNAs participate in Tartary buckwheat seed development remains unexplored.

Results

In this study, we first identified 26 miRNA biosynthesis genes in the Tartary buckwheat genome and described their phylogeny and expression profiling. Then we performed small RNA (sRNA) sequencing for Tartary buckwheat seeds at three developmental stages to identify the miRNAs associated with seed development. In total, 230 miRNAs, including 101 conserved and 129 novel miRNAs, were first identified in Tartary buckwheat, and 3268 target genes were successfully predicted. Among these miRNAs, 76 exhibited differential expression during seed development, and 1534 target genes which correspond to 74 differentially expressed miRNAs (DEMs) were identified. Based on integrated analysis of DEMs and their targets expression, 65 miRNA-mRNA interaction pairs (25 DEMs corresponding to 65 target genes) were identified that exhibited significantly opposite expression during Tartary buckwheat seed development, and 6 of the miRNA-mRNA pairs were further verified by quantitative real-time polymerase chain reaction (qRT-PCR) and ligase-mediated rapid amplification of 5′ cDNA ends (5′-RLM-RACE). Functional annotation of the 65 target mRNAs showed that 56 miRNA-mRNA interaction pairs major involved in cell differentiation and proliferation, cell elongation, hormones response, organogenesis, embryo and endosperm development, seed size, mineral elements transport, and flavonoid biosynthesis, which indicated that they are the key miRNA-mRNA pairs for Tartary buckwheat seed development.

Conclusions

Our findings provided insights for the first time into miRNA-mediated regulatory pathways in Tartary buckwheat seed development and suggested that miRNAs play important role in Tartary buckwheat seed development. These findings will be help to study the roles and regulatory mechanism of miRNAs in Tartary buckwheat seed development.

Peer Review reports

Background

MicroRNAs (miRNAs), a class of endogenous non-coding small RNAs (sRNAs) that are 20–24 nt in length, play crucial regulatory functions in animals and plants by repressing their target genes expression at the transcription or post-transcription level [1,2,3,4]. In plants, the biosynthesis of miRNAs is a multistep pathway that is involved in many genes and enzymes [5]. Plant miRNA genes are initially transcribed into primary miRNAs (pri-miRNAs) by RNA polymerase II in the nucleus [5, 6]. Then, the pri-miRNAs are processed into miRNA precursors (pre-miRNAs) with stem-loop structures by the DICER-LIKE1 (DCL1), HYPONASTIC LEAVES1 (HYL1) and SERRATE (SE) protein complex [5,6,7]. Next, these pri-miRNAs are further cleaved into miRNA::miRNA* duplexes under the action of DCL1, HYL1 and SE protein complex, and then the 3΄ end of the duplexes are methylated by the methyltransferase HUA ENHANCER1 (HEN1) [5]. After methylation, the duplexes are exported to the cytoplasm by the HASTY protein [5, 8]. Finally, the miRNA duplex is bound by ARGONAUTE1 (AGO1) to form the RNA-induced silencing complex (RISC) to carry out its function by either cleaving target mRNAs or repressing the translation process [5, 9]. At present, miRNAs have been identified in numerous plants, and more and more evidence indicates that they play crucial roles in plant growth and development, secondary metabolism, biotic and abiotic stress tolerance, and signal transduction [4, 10,11,12,13,14,15,16].

Seed is the reproductive and the primary nutrient storage organ in many crop plants, and its developmental success or failure directly determines the final crop yield and seed quality, as well as whether genetic information can be successfully transmitted to the next generation [4]. Crop seed development is terribly a complex biological process that involves many gene regulatory pathways [17, 18]. An increasing body of evidence shows that miRNAs participate in the regulation of seed development in crop plants. To date, thousands miRNAs have been identified in the development of seed in multiple crop plants including rice [19, 20], maize [21, 22], wheat [23, 24], barley [25, 26], soybean [27, 28], peanut [18], Brassica napus [29, 30], narrow-leafed lupin [6], and common buckwheat [4] by using a high-throughput sequencing approach. These studies found that the expression of various miRNAs is extremely dynamic during seed development and some miRNAs are specifically expressed in seed, implying that miRNAs have very vital regulatory roles in seed development. In fact, a few miRNAs have been functionally demonstrated to play crucial regulatory roles in the seed development of model crop rice and a few other crops by negatively regulating their target genes expression. For example, miR159, miR160, miR397, miR398, and miR408 positively regulate rice grain size [31,32,33,34,35], while miR1432, miR156, miR167, miR396c, miR396e, miR396f, and miR530 have opposite roles in regulating the grain size [35,36,37,38,39,40,41]. Notably, the conserved role of miR408 and miR160 in regulating seed size was also found in Arabidopsis and tobacco, and cotton, respectively [42, 43]. In addition to the role in seed size, miRNAs were also demonstrated to regulate nutrient accumulation in developing seed [4]. For instance, miR160 positively regulate starch accumulation in rice seed [4, 32]. EgmiR5179 and csa-miR167A regulate the oil and linolenic acid biosynthesis in oil palm and Camelina sativa seeds, respectively [44, 45].

Tartary buckwheat (Fagopyrum tataricum) is an annual medicinal and edible crop, belonging to the eudicot family Polygonaceae [46]. It is widely cultivated in Asia and Eastern Europe, especially in the mountainous areas of Southwest China [46]. Tartary buckwheat seed is a good source of nutrients including starch, protein, dietary fiber, fatty acid (linoleic acid), and various minerals [46]. Importantly, Tartary buckwheat seed also contains rich flavonoids especially rutin, which have been proved to be effective in preventing liver injury and especially inflammatory liver injury [43, 46]. Therefore, it is of great significance to understand the molecular mechanism of Tartary buckwheat seed development, which will be helpful in high-yield and quality breeding of Tartary buckwheat [46]. To date, several transcriptome analyses have been reported that examine the molecular mechanism of Tartary buckwheat seed development [4, 43, 46,47,48]. However, to our knowledge, no one has studied miRNAs in Tartary buckwheat, and miRNAs whether and how to regulate the development of Tartary buckwheat seed is largely unclear. In this study, we first examined the conserved evolution of miRNA biosynthesis in Tartary buckwheat compared to other plants through homology identification and phylogenetic analysis of miRNA biosynthesis genes. Then, we identified the known and novel miRNAs in the developing seed of Tartary buckwheat and predicted their target genes. Finally, we performed integration analysis between miRNA and mRNA expression to insight into the miRNA-mediated molecular mechanisms of Tartary buckwheat seed development and identified key miRNA-mRNA interaction pairs for Tartary buckwheat seed development. Our results provide valuable information for enhancing the understanding of the miRNA-mediated regulatory mechanism of Tartary buckwheat seed development and aid in the Tartary buckwheat seed improvement.

Results

Identification, phylogeny, and expression profiles analysis of the miRNA biosynthesis genes in Tartary buckwheat

RDR, DCL, HYL, SE, HEN, HST and AGO genes have been demonstrated to play essential roles in plant miRNA biosynthesis [5]. As the first and most important step in studying the miRNA in Tartary buckwheat, we identified the orthologs of these genes in Tartary buckwheat. A total of 8 RDR, 4 DCL, 1 HYL, 2 SE, 1 HEN, 2 HST, and 8 AGO genes were identified in the Tartary buckwheat genome (Additional file 1: Table S1), respectively. Phylogenetic analysis revealed that plant RDR proteins can be divided into four clades, as previously defined by Qian et al. [49] (Fig. 1). Among the 8 FtRDR proteins, 3, 1, 3, and 1 belong to clades 1, 2, 3, and 4, respectively. In addition, the FtRDR proteins were closely homologous to A. thaliana RDR proteins. The DCL proteins can also be classed into four groups, and 4 FtDCL proteins showed a one-to-one counterpart with A. thaliana DCL proteins (Fig. 1). Similarly, the AGO proteins can be separated into four clades, but clade 4 only contained the grass AGO proteins as defined by Zhang et al. [50] (Fig. 1). For Tartary buckwheat AGO proteins, clade 1 was the biggest clade, which contained 6 FtAGO proteins. In contrast, both clades 1 and 2 only contained 1 FtAGO protein (Fig.1).

Fig. 1
figure1

Phylogenetic analysis of Tartary buckwheat RDR, DCL and AGO proteins. The corresponding proteins of Arabidopsis, maize, and rice were also used to construct the phylogenetic tree. The phylogenetic tree was constructed by using MUSCLE and MEGA 7.0 with the Maximum Likelihood method

RNA-seq data were further used to investigate the expression profiles of the 26 identified genes in roots, stems, leaves, flowers, and during seed development. As shown in Fig. 2, most RDR genes (6) had lower expression in all tissues. In contrast, 4 genes (FtAGO1, FtAGO2, FtAGO6, and FtSE1) had constitutive high expression in all tissues. The remaining 16 genes were expressed at a moderate level. Among them, 4 genes (FtHYL1, FtHST2, FtAGO4, and FtAGO8) exhibited specifically high expression in seeds. Furthermore, 6 genes (FtDCL1, FtDCL3, FtDCL4, FtAGO2, FtAGO3, and FtRDR4) showed significant differential expression (|log2(fold change)| > 1 and FDR value < 0.05) during seed development (Fig. 2). The 6 differentially expressed genes (DEGs) showed three expression patterns in the developing Tartary buckwheat seed (1) FtAGO2 and FtAGO3 were down-regulated only at the initial maturity stage (S3), (2) FtDCL1 and FtRDR4 showed a sustained decrease during seed development, and (3) FtDCL3 and FtDCL4 were up-regulated at the peak filling stage (S2) and down-regulated at the initial maturity stage (S3) (Fig. 2).

Fig. 2
figure2

Expression profiles of Tartary buckwheat miRNA biosynthesis genes in different organs and during seed development. The different colors represent the expression abundance of the genes

Sequencing of Tartary buckwheat seed sRNAs

To investigate the effects of miRNAs on seed development in Tartary buckwheat, we constructed and sequenced six sRNA libraries from seeds at three different developmental stages (Fig. 3). A total of 120.63 million raw reads were generated in six sRNA libraries (Additional file 1: Table S2). After filtration, 13.89 million, 17.09 million, 17.19 million, 14.54 million, 13.64 million, and 12.00 million clean reads were obtained for the six libraries, respectively (Additional file 1: Table S2). The length distribution displayed most the small RNA reads were 21–24 nt (Additional file 2: Figure S1). Among them, the 24 nt small RNA was the most abundant type and showed differential accumulation during Tartary buckwheat seed development (Additional file 2: Figure S1). To further obtain sRNA reads containing miRNAs, non-coding RNAs (including rRNAs, tRNAs, snRNAs, and snoRNAs) and repeat sequences were removed by comparing the clean reads to the Silv, GtRNAdb, Rfam, and Repbase databases. After discarding the non-coding RNAs and repeat sequences, a total of 10.45 million, 14.16 million, 14.60 million, 13.05 million, 6.80 million, and 5.63 million unannotated sRNA reads containing miRNAs were obtained for the six libraries, respectively (Additional file 1: Table S2). Among these unannotated clean reads, a total of 61.94, 66.76, 46.98, 39.20, 37.22, and 39.42% reads were mapped to the Tartary buckwheat reference genome, respectively (Additional file 1: Table S2).

Fig. 3
figure3

Seed phenotypes of Tartary buckwheat at the initial filling stage (S1), peak filling stage (S2), and initial maturity stage (S3)

Identification of known and novel miRNAs

To identify known miRNAs in the developing Tartary buckwheat seed, the mapped sRNA reads were subjected to Blastn search against miRBase v21.0. In total, 101 conserved miRNAs, belonging to 25 known miRNA families, were identified (Additional file 1: Table S3). Of these families, 17 families contained more than two members, three families (MiR390, MiR399, and MiR858) contained two members, and the remaining five families (MiR159, MiR394, MiR397, MiR828, and MiR845) only contained one member. The length of these known miRNAs ranged from 19 to 22 nt, and the numbers of miRNAs with length of 19, 20, 21, and 22 nt were 8, 29, 56, and 8, respectively (Additional file 1: Table S3). Among these identified conserved miRNAs, 29 miRNAs were highly expressed, which have more than 1000 read counts at least in one library (Additional file 1: Table S3). The rest of miRNAs were expressed at a moderate level (36) or low level (36) which have < 10 read counts for each library, respectively (Additional file 1: Table S3). Notably, the highest expression miRNA was fta-miR159, and read counts ranged from 14,806 to 95,396 for each library.

To further identify novel miRNAs in the developing Tartary buckwheat seed, the remaining sRNA sequences that were not mapped to the Tartary buckwheat reference genome were analyzed using the miRdeep2 program. As a result, a total of 129 novel miRNAs were identified (Additional file 1: Table S4). The length of these novel miRNAs ranged from 18 nt to 24 nt, and the minimum free energy (AMFE) distribution ranged from − 116.2 kcal moL− 1 to − 26.2 kcal moL− 1 (Additional file 1: Table S4). Among these novel miRNAs, 66 miRNAs were assigned to 35 known miRNA families in the miRNA database, while the remaining 63 miRNAs had no similarity to any known family. Like the above identified conserved miRNAs, more than two-thirds of the novel miRNAs were expressed at a moderate or high level, and the highest expression was fta_novel_miR71 (Additional file 1: Table S4).

miRNAs target prediction and functional analysis

To understand the potential function of these identified miRNAs, the putative target genes of the 230 miRNAs were predicted using the TargetFinder software. A total of 3268 potential target genes were successfully predicted for among 213 miRNAs, including 2052 target genes from 101 known miRNAs and 1372 target genes from 112 novel miRNAs (Additional file 1: Table S5). These predicted target genes mainly encoded transcription factors (TFs, 285), protein kinases (254), phosphatase, E3 ubiquitin-protein ligase, proteins in hormone signal transduction and other cellular processes, and enzymes in various metabolisms (Additional file 1: Table S6). The main target TFs of these miRNAs are displayed in Additional file 1: Table S7, and the top three TFs were MYB (65), AP2/ERF (18), and NAC (18), respectively. Consistent with the results of previous studies, some conserved miRNAs targeted the known TFs. For instance, miR156, miR160, miR164, miR171, miR172, miR319, miR396, and miR858 targeted SPL, ARF, NAC, GRAS, AP2-ERF, TCP, GRF, and MYB TFs, respectively (Additional file 1: Table S7). In addition, based on a homologous search of the putative target genes, the predicted target genes of 61 miRNAs, including 39 conserved and 22 novel miRNAs, were found to be the homologous genes of 47 known seed or organ size genes (Additional file 1: Table S8). Among these, the 39 conserved miRNAs were from 14 known miRNA families including MiR156/157 (10), MiR159 (1), MiR162 (1), MiR164 (1), MiR166 (2), MiR167 (1), MiR168 (4), MiR169 (1), MiR172 (5), MiR319 (2), MiR390 (2), MiR395 (2), MiR396 (4), and MiR530 (3). Notably, among these 61 miRNAs, most miRNAs were first found to target the known seed or organ size genes, except the members of the MiR156/157 and MiR396 families which are well known to target the seed size SPL and GRF TFs, respectively. Furthermore, based on a homologous search, 8 miRNAs (4 known and 4 novel) were identified to target the structural or regulatory genes of flavonoid biosynthesis (Additional file 1: Table S9). Among these, fta_novel_miR26 targeted phenylalanine ammonia-lyase (PAL), fta_novel_miR74 targeted 4-coumarate-CoA ligase (4CL), fta-miR394a-5p targeted flavonol synthase (FLS), fta_novel_miR1 targeted the homologous genes of AtMYB123 that regulate anthocyanindins biosynthesis, fta_novel_miR58 targeted the homologous genes of AtbHLH42/TT8 that regulate anthocyanindins biosynthesis, and fta-miR828-5p targeted the homologous genes of AtMYB75/90/113/114 that regulate anthocyanindins biosynthesis. In addition, fta-miR858a-5p and fta-miR858b-3p targeted 13 and 14 MYB TFs involved in the regulation of flavonol biosynthesis and anthocyanindins biosynthesis, respectively.

DEMs in the development Tartary buckwheat seed and their targets

To identify Tartary buckwheat seed development-associated miRNAs and understand their potential regulatory mechanisms, the DEMs were identified by comparing the TPM expression value. A total of 76 miRNAs, including 39 conserved and 37 novel miRNAs, displayed significant differential expression (Fig. 4). Among these, 45, 43, and 22 miRNAs were found in the comparisons between S1 and S2, S1 and S3, and S2 and S3, respectively (Fig. 4a). When compared with S1, 1 miRNA (fta_novel_miR110) was differentially expressed in both S2 and S3 stages, and 22 and 17 miRNAs were specifically differentially expressed in the S2 and S3 stage, respectively (Fig. 4a). When compared with S2, 4 miRNAs were specifically differential expression in S3 stage (Fig. 4a). In addition, between stages S1 and S2, 19 and 26 miRNAs were up- or down-regulated, respectively (Fig. 4b). In the S1 vs. S3 comparison, 23 up-regulated and 20 down-regulated miRNAs were found. In the comparison of S2 vs. S3, 11 and 11 miRNAs were up- or down-regulated, respectively (Fig. 4b).

Fig. 4
figure4

DEMs during Tartary buckwheat seed development. a Venn diagram represented the overlap of DEMs among the comparisons. b Number of DEMs in the comparisons of S1 vs. S2, S1 vs. S3, and S2 vs. S3. Green and yellow bars represented the number of up- or down-regulated genes, respectively. c Expression heat map of DEMs. The different colors represent the expression abundance of the miRNAs

The expression heat map of all DEMs as shown in Fig. 4c. Among these DEMs, 26, 18, and 9 miRNAs displayed specifically high expression at the initial filling stage, peak filling stage, and initial maturity stage, respectively. Twelve miRNAs exhibited high expression at both the initial filling stage and peak filling stage. Notably, the expression of 9 novel miRNAs (fta_novel_miR17–25) was sustained growth during seed development (Fig. 4c). To verify the results of miRNA sequencing, quantitative stem-loop RT-PCR was performed to determine the expression level of 9 DEMs during Tartary buckwheat seed development. As shown in Fig. 5, the results were broadly consistent with those obtained in miRNA sequencing.

Fig. 5
figure5

Quantitative real-time PCR verification of 9 DEMs at different developmental stages in Tartary buckwheat seeds

To better understand the functions of these DEMs, their predicted target genes were investigated and further subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. A total of 1534 target genes were identified for among 74 DEMs (Additional file 1: Table S5). Of these target genes, 850, 765, and 735 genes were targeted by 44, 41, and 22 DEMs in the S1 vs. S2, S1 vs. S3, and S2 vs. S3 comparisons, respectively. GO analysis showed that 514, 448, and 419 targets were assigned to biological processes, cellular component, and molecular function categories for the S1 vs. S2, S1 vs. S3, and S2 vs. S3 comparisons, respectively (Additional file 2: Figure S2). There were 18, 14, and 14 GO terms in biological processes, cellular component, and molecular function categories for each comparison, respectively. The “metabolic process”, “cell part”, and “catalytic activity” were the greatest abundance terms for these three categories (Additional file 2: Figure S2). KEGG analysis displayed that a total of 332, 309, and 272 target genes were assigned to 61, 54, and 55 KEGG pathways for the S1 vs. S2, S1 vs. S3, and S2 vs. S3 comparisons (Additional file 1: Table S10), respectively. Most of these pathways were involved in the metabolism process, of which biosynthesis of “amino acids biosynthesis” (ko01230), “starch and sucrose metabolism” (ko00500) “phenylpropanoid biosynthesis” (ko00940), “purine metabolism” (ko00230), and “amino sugar and nucleotide sugar metabolism” (ko00520) were the major pathways. In addition, a high number target genes were assigned to “plant hormone signal transduction” (ko04075). The top 50 pathways for each comparison are shown in Additional file 2: Figure S3.

Identification of the key miRNA-mRNA pairs related to Tartary buckwheat seed development

To identify potential miRNA-mRNA pairs related to Tartary buckwheat seed development, we performed expression correlation analyses between DEMs and target mRNAs that were differentially expressed during Tartary buckwheat seed development. Based on our previous transcriptome data [46], which have the same samples with small RNA sequencing, 439 out of 1534 target genes of 74 DEMs exhibited significant differential expression during Tartary buckwheat seed development. Correlation analyses between the 439 DEGs and their corresponding DEMs identified that 117 miRNA-mRNA pairs were negatively correlated (Additional file 1: Table S11). Among them, 65 miRNA-mRNA pairs, which consisted of 25 miRNAs and 65 target genes, showed significantly expression negative correlation (R ≥ 0.5, P < 0.05) (Table 1 and Fig. 6). Based on a homologous annotation of these 65 target genes, 56 miRNA-mRNA pairs achieved homologous functional annotation (Table 1). These 56 miRNA-mRNA pairs were found to be significantly involved in cell differentiation and proliferation, cell elongation, hormone response and balance, organogenesis, embryo and endosperm development, transport of mineral elements (Pi, Fe and Mn), and fatty acid and flavonoid biosynthesis (Table 1). In addition, 2 calcium binding proteins, 1 receptor-like kinase, 1 CBL-interacting protein kinase, and 1 SUMO protease were identified in these 65 target genes, although their homologues in other plants were not reported to participate in seed development (Table 1). Notably, among these identified miRNA-mRNA pairs, 5 miRNA-mRNA pairs, including fta-miR156c-5p-FtPinG0000496000.01, fta-miR167c-5p-FtPinG0002560000.01, fta-miR396a-5p-FtPinG0000680700.01, fta-miR530c-5p-FtPinG000259–4600.01, and fta_novel_miR17-FtPinG0007576600.01, might be involved in seed size regulation of Tartary buckwheat because the target mRNAs were homologous to the known genes related to rice seed size (Table 1). Furthermore, fta-miR858a-5p exhibited conspicuous negative correlation with several MYB TFs which are the well-known regulators of plant flavonoid biosynthesis (Table 1).

Table 1 Identified key miRNA-mRNA interaction pairs for tartary buckwheat seed development
Fig. 6
figure6

miRNA-mRNA interaction pairs showed significant negative correlation (R ≥ 0.5, P < 0.05) of expression during Tartary buckwheat seed development. Left: Heat map of DEMs. Right: Heat map of the differentially expressed target mRNAs of DEMs. The different colors represent the expression abundance of the miRNAs and their corresponding target mRNAs

qRT-PCR and 5′-RLM-RACE validation of the key miRNA-mRNA pairs related to Tartary buckwheat seed development

To further confirm these identified key miRNA-mRNA pairs related to Tartary buckwheat seed development, qRT-PCR and 5′-RLM-RACE were performed to analyze fta-miR156c-5p-FtPinG0000496000.01, fta-miR164a-5p-FtPinG0000240100.01, fta-miR166c-3p-FtPinG0008845600.01, fta-miR167c-5p-FtPinG0002560000.01, fta-miR396a-5p-FtPinG0000680700.01, and fta_novel_miR118-FtPinG0003617400.01 pairs, which were involved in seed size, cell differentiation and proliferation, and auxin response, respectively. qRT-PCR results indicated that the expression of these miRNAs was obviously negatively correlated with their corresponding target genes expression (Fig. 7). 5′-RLM-RACE results showed that these miRNAs could cleave their corresponding target genes (Fig. 8).

Fig. 7
figure7

Quantitative real-time PCR validation of six identified key miRNA-mRNA interaction pairs for Tartary buckwheat seed development. Green and blue bars represented the miRNAs and their corresponding target mRNAs, respectively. The error bar represents the error values of three biological replicates

Fig. 8
figure8

5′ RLM-RACE verification of six identified key miRNA-mRNA interaction pairs for Tartary buckwheat seed development. The numbers above sequences indicate the detected cleavage site of independent clones

Discussion

MiRNAs, a class negative regulators of gene expression, play key roles in numerous plant developmental processes including seed development [10, 51]. However, it whether and how to regulate the development of Tartary buckwheat seed remains unexplored. Here, to better understand the miRNAs-mediated molecular mechanisms of Tartary buckwheat seed development, we systematically identified the miRNA biosynthesis genes in the Tartary buckwheat genome and miRNA during Tartary buckwheat seed development. More importantly, we performed integrated analysis of miRNA and target mRNA expression profiles during Tartary buckwheat seed development, and identified the key miRNA-mRNA pairs for Tartary buckwheat seed development.

miRNA biosynthesis orthologs within the Tartary buckwheat genome

Numerous studies have shown that several different genes, including RDR, DCL, HYL, SE, HEN, HST, and AGO, play crucial roles in plant miRNA biosynthesis [5]. In this study, we identified 8 RDR, 4 DCL, 1 HYL, 2 SE, 1 HEN, 2 HST, and 8 AGO genes in the Tartary buckwheat genome. The phylogenetic analyses showed that 8 RDR, 4 DCL, and 8 AGO could be divided into 4, 4, and 3 subfamilies, respectively. The results were consistent with the previous studies in other plants [6, 49, 50, 52, 53], and indicated that the evolution of these gene families was conserved in these plants. Expression analysis revealed that 1 AGO-like FtAGO1 had the highest expression in seeds in these identified genes, 1 AGO-like FtAGO4 exhibited specific high expression in seeds, and 2 AGO-like FtAGO2 and FtAGO3 exhibited significant differential expression during seed development. Notably, OsAGO17, which belongs to the same subfamily as these four FtAGO genes, had been functionally demonstrated to positively regulate the grain size and grain weight in rice [54, 55]. This suggested that these four FtAGO genes might also play a similar role with OsAGO17 in Tartary buckwheat seed development. In addition, we also found that FtHYL1, FtHST2, and FtAGO8 displayed especially high expression in seeds, and FtDCL1, FtDCL3, FtDCL4, and FtRDR4 were obviously differential expression during seed development. These indicated that these miRNA biosynthesis genes might also have important regulatory roles in Tartary buckwheat seed development.

Characteristics of sRNAs during Tartary buckwheat seed development

By sequencing six sRNA libraries from three differently developmental stages seeds, we obtained abundant sRNAs with a length of 18–30 nt. Among them, the 24 nt sRNAs were the most abundant in developmental Tartary buckwheat seed, which was similar to previous observations in the developmental seed of many plants [6, 56]. Furthermore, the 24 nt sRNAs exhibited differential accumulation during Tartary buckwheat seed development. These suggested that the 24 nt sRNAs might play crucial roles in Tartary buckwheat and other plant seed development. It has been reported that many of the 24-nt sRNAs were heterochromatic siRNAs (hetsiRNAs), which mediate transcriptional gene silencing through DNA methylation (RdDM) [56]. Generally, the 24-nt sRNAs need from cytoplasm imported into the nucleus before they were methylated, and the process was mediated by AGO4 [57]. Notably, our study observed that FtAGO8, which was the AGO4 ortholog, exhibited special high expression in the development Tartary buckwheat seed. These observations implied that FtAGO8 might has a role similar to AGO4 in 24-nt sRNAs methylation, and these 24-nt sRNAs might mediate transcriptional gene silencing through the RdDM pathway during Tartary buckwheat seed development.

Known and novel miRNAs and their target genes in the developmental Tartary buckwheat seed

In the Tartary buckwheat genome, 278 miRNAs have been predicted through a genome-wide bioinformatics analysis [58]. In our study, a total of 230 miRNAs, including 101 known and 129 novel miRNAs, were identified during Tartary buckwheat seed development. Among the miRNAs, more than two-thirds were expressed at high and moderate levels. These indicated that most miRNAs in the Tartary buckwheat genome were involved in the seed development, and further suggested that Tartary buckwheat seed development was an extremely complex biological process. Notably, the fta-miR159 exhibited the highest expression in the developmental Tartary buckwheat seed. In strawberry, Fa-miR159b and Fa-miR159b were more highly expressed in developing fruit and played a key regulatory role in fruit development [4, 59]. In rice, osa-miR159 was also more highly expressed in developing seed and positively regulated grain length and width [4, 20, 33, 35]. Therefore, our results suggested that fta-miR159 might also possess a crucial regulation role in the Tartary buckwheat seed development.

Using TargetFinder software, we identified 3268 potential target genes from 213 miRNAs. Among these target genes, the largest number were encoded TFs. Notably, among these miRNAs with target gene encoded TFs, some miRNAs such as miR156, miR160, miR396, and miR858 targeted the SBP, ARF, GRF, and MYB TFs, respectively. Interestingly, these miRNAs have been experimentally verified to target these corresponding TFs in other plants [32, 36, 37, 60]. This suggested that these miRNA targets were highly conserved in different plants and further confirmed the high reliability of target identification in our study. Based on a homologous search of the predicted target genes, 61 miRNAs have target genes were homologous genes of 47 known seed or organ size genes [61]. Notably, among these miRNAs, most were first identified targeting the known seed or organ size genes, except for miR156 and miR396 targeting the known seed size SPL and GRF TFs, respectively [22, 32, 36]. These showed that these miRNAs might play an important regulatory role in Tartary buckwheat seed size and most of these known seed or organ size regulatory genes might exist a post-transcription regulatory mechanism in regulating seed size. In plants, miRNAs have been reported to regulate flavonoid biosynthesis through target regulating the structural or regulatory genes of flavonoid biosynthesis [62]. In our study, 3 miRNAs (fta_novel_miR26, fta_novel_miR74, and fta-miR394a-5p) were identified to target the flavonoid biosynthesis structural genes PAL, 4CL, and FLS, respectively. Furthermore, 5 miRNAs (fta_novel_miR1, fta_novel_miR58, fta-miR858a-5p, fta-miR858a-5p, and fta-miR858b-3p) were found to target the orthologs of known flavonoid biosynthesis regulatory genes in Arabidopsis thaliana [63,64,65,66]. These implied that these miRNAs could regulate flavonoid biosynthesis in Tartary buckwheat by regulating the expression of the structural or regulatory genes of flavonoid biosynthesis.

DEMs during Tartary buckwheat seed development

In many plants, a larger number of DEMs have been identified in the developing seed [4, 6, 18,19,20,21,22,23,24,25,26,27,28,29,30]. In this study, we identified 76 DEMs during Tartary buckwheat seed development. Expression profiles analysis revealed that these DEMs possessed stage-specific highly expressed patterns during Tartary buckwheat seed development, suggesting that these DEMs might perform their regulatory roles in a specific stage during Tartary buckwheat seed development and miRNAs with the same expression pattern might have a similar regulatory role in the developing Tartary buckwheat seed. Notably, some DEMs such as miR156, miR160, miR166, miR167, miR168, miR395, miR396, miR397, miR398, miR399, and miR408 were also found to be differentially expressed during seed development in many seed plants [4, 18, 19, 21,22,23,24,25,26,27,28,29,30]. In rice, miR156, miR160, miR167, miR396, miR397, miR389, and miR408 have been functionally verified to regulate the rice seed development [31,32,33,34,35,36,37,38, 40]. These suggested that these known miRNAs were the key and conserved miRNAs in different plant seed development regulation and could also be used as candidate miRNAs for Tartary buckwheat seed development. In addition, 37 novel miRNAs and several conserved miRNAs, including miR530, miR828, and miR858, were found to be specifically differential expression in developing Tartary buckwheat seed when compared with other seed plants. This indicated that these miRNAs might have a specific regulatory function in Tartary buckwheat seed development. Furthermore, all these information indicated that conserved and diverse miRNA-mediated regulatory mechanisms in seed development might exist in different seed plants.

Integration analysis of miRNA and target mRNA

Integrated miRNA and its target mRNA expression analysis could be helpful in understanding the function of miRNAs and identifying the functional miRNA-mRNA pairs related to seed development [67]. To our knowledge, there was no integration analysis of miRNA and target mRNA performed in previous studies of plant seed development. In this study, we identified 65 significantly negatively correlated miRNA-mRNA pairs (R ≥ 0.5, P < 0.05), which consisted of 25 miRNAs and 65 corresponding target genes. Based on homologous queries of these miRNAs corresponding target genes, 56 miRNA-mRNA pairs obtained functional annotations. These miRNA-mRNA pairs are involved in different aspects, most significantly cell differentiation and proliferation, cell elongation, hormones response and balance, organogenesis, embryo and endosperm development, mineral elements transport and accumulation, and flavonoid biosynthesis. These suggested that these identified miRNA-mRNA pairs were the key miRNA-mRNA pairs related to Tartary buckwheat seed development.

In rice, OsSPL16 [68, 69], OsNAC024 [70], OsSUS [71], OsGIF1 [38, 72], and OsAGO17 [54] genes have been functionally demonstrated to regulate the seed size. In this study, we found the target mRNAs from 5 miRNA-mRNA pairs (fta-miR156c-5p-FtPinG0000496000.01, fta-miR167c-5p-FtPinG0002560000.01, fta-miR396a-5p-FtPinG0000680700.01, fta-miR530c-5p-FtPinG0002594600.01 and fta_novel_miR17-FtPinG0007576600.01) were homologous with these 5 rice seed size genes, respectively, which indicated that these genes might have a conserved regulatory role in the seed size of Tartary buckwheat seed and the corresponding miRNAs were the key miRNAs for controlling Tartary buckwheat seed size. Notably, to date, these seed size genes had not been identified as the target genes of any miRNAs in previous studies. Therefore, our results provided the first-hand information on the post-transcription regulation of these seed size genes. It is well known that the Ca2+ signal transduction pathways play important roles in many developmental processes including seed development in plants [73, 74]. In this study, we found that three of the 56 identified miRNA-mRNA pairs were involved in Ca2+ signal transduction because the target mRNAs were homologous to the calcium binding protein and CBL-interacting protein kinases, respectively. This suggested that these miRNAs could regulate Tartary buckwheat seed development though Ca2+ signal transduction pathways. In a previous study, we found that the total seed flavonoid content was dynamically accumulated during Tartary buckwheat seed development, and identified one SG7 subgroup R2R3-MYB TF gene (FtPinG0009153900.01) was the key regulatory gene of flavonoid biosynthesis in the developing Tartary buckwheat seed [46]. Interestingly, in this study, we found fta-miR858a-5p-FtPinG0009153900.01 pair displayed significantly expression negative correlation. In addition, fta-miR858a-5p also exhibited conspicuous negative correlation with the orthologs of Arabidopsis thaliana MYB123, which was the well-known regulator of proanthocyanidin accumulation in developing seed [64]. These indicated that the flavonoid biosynthesis of Tartary buckwheat seed is involved in miRNA-mediated post-transcription regulation and fta-miR858a-5p is the key regulator.

To further verify the reliability of these identified miRNA-mRNA pairs related to Tartary buckwheat seed development, we carried out qRT-PCR and 5′-RLM-RACE analyses for among 6 miRNA-mRNA pairs that are involved in seed size, cell differentiation and proliferation, and auxin response. These results showed that the expression of these miRNAs and their corresponding target genes were obviously negatively related and could cleave their corresponding target genes. The above results indicated that these identified miRNA-mRNA pairs could be considered as the candidate miRNA-mRNA pairs for regulating the development of Tartary buckwheat seed, and also suggested that combining small RNA and transcriptome analysis is an effective method for identifying key miRNAs that are involved in plant seed development.

Conclusion

The present work is the first attempt to integrate miRNA and mRNA expression data to identify key regulatory miRNA-mRNA pairs in developing Tartary buckwheat seed. A total of 230 miRNAs, including 101 known and 129 novel miRNAs, were first identified in Tartary buckwheat. Among these miRNAs, 76 showed significant differential expression during Tartary buckwheat seed development, and 1543 of their target genes were identified. Additionally, 25 miRNAs corresponding to 56 target genes were identified as key candidate miRNA-mRNA pairs for Tartary buckwheat seed development. The integrated analyses of miRNAs and target mRNAs in this study not only provide the first new insights into miRNA-mediated regulation in Tartary buckwheat seed development, but also provide a basis for further research into the functions of these candidate miRNAs and their targets in Tartary buckwheat seed development and for Tartary buckwheat seed improvement.

Methods

Plant materials and sample preparation

Tartary buckwheat cultivar Jinqiao No. 2, which was obtained from the Research Center of Buckwheat Industry Technology of Guizhou Normal University (Guiyang, Guizhou, China), was used in this study. The seed samples used for sRNA sequencing were same as our previous transcriptome analysis [46]. In brief, seeds were collected at 12 (S1: initial filling stage), 15 (S2: peak filling stage), and 21 (S3: initial maturity stage) days after pollination (DAP) [46].

Identification, phylogeny, and expression profiles analysis of the miRNA biosynthesis genes in the Tartary buckwheat genome

The Tartary buckwheat genome information was downloaded from http://www.mbkbase.org/Pinku1 [58]. All protein sequences from annotated genes were used to build a local protein blast database in BioEdit software (Version 7.0.9.0). The protein sequences of RDR16, DCL14, AGO110, HEN1, HYL1, SE, and HST genes from Arabidopsis thaliana were used as queries to carry out local blastp with e-values <e− 5 to identify orthologs in the Tartary buckwheat genome. The conserved domain of all obtained non-redundant sequences were further examined in the Pfam (http://pfam.sanger.ac.uk) [75], and the final orthologs were ensured. The RDR, DCL, and AGO protein sequences from Tartary buckwheat, Arabidopsis thaliana, rice, and maize were used to construct the phylogenetic trees using the Maximum Likelihood method in MEGA 7.0.21 software. Expression analysis of these identified miRNA biosynthesis genes in Tartary buckwheat was performed using RNA-seq data from 4 organs (root, stem, leaf, and flower) [58] and three developmental-stage seeds [46]. The expression heatmaps of these genes were generated by using TBtools software [76].

RNA extraction, sRNA library construction, and sequencing

Total RNA was extracted from three samples as previously described [46] and used for sRNA and qRT-PCR assay. RNA quality and concentration were determined via 1.2% agarose gel electrophoresis and NanoDrop 2000c spectrophotometer (NanoDrop, Wilmington, DE, USA), respectively.

For construction and sequencing of the sRNA libraries, two replicates were performed for each sample. Firstly, sRNAs were isolated from total RNA by polyacrylamide gel electrophoresis (PAGE). Next, the isolated sRNAs were added to a 5′ RNA adaptor and a 3′ RNA adaptor by using T4 RNA ligase (TaKaRa, Dalian, China). Then, sRNAs with added 5′ and 3′ RNA adaptors were reverse transcribed into single-stranded cDNA using RT-PCR. Follows, the single-stranded cDNA was further synthesized into double-stranded cDNA by PCR amplification using adapter primers. Finally, the PCR product was purified and subjected to high-throughput sequencing by using the Illumina SE50 system at Biomarker Technologies Co., Ltd. (Beijing, China).

Bioinformatics analysis of sequencing data

Raw reads were processed to generate clean reads by removing the low-quality reads, reads containing poly-N or poly-A sequences, reads lacking the 3′ adaptor sequence, reads with length < 18 nt or > 30 nt, and adapter sequences. The clean reads were then mapped against the Tartary buckwheat genome (http://www.mbkbase.org/Pinku1/) [58]. The Silva database, GtRNAdb database, Rfam database and Repbase database were used to filter the rRNA, tRNA, snRNA, and snoRNA to produce the unannotated reads containing miRNA [4]. Unannotated reads were blasted against miRBase to search for known miRNAs [77]. The potential novel miRNAs were predicted by using the miRdeep2 program [78]. The miFam.dat (http://www.miRbase.org/ftp.shtml) was used to investigate the miRNA family class of these identified miRNAs. The miRNA expression abundances were calculated and normalized by transcript per million (TPM), and differential expression analysis was carried out by the DEGSeq R package with FDR (false discovery rate) value < 0.05 and |log2(fold change)| > 1 as the threshold for significant difference. The potential target genes of miRNA were predicted by using TargetFinder software [79]. The predicted target genes were subjected to GO (http://geneontology.org/) [80] and KEGG (www.kegg.jp/kegg/kegg1. html) [81] pathway analysis to predict and classify possible functions.

Identification of the key miRNA-mRNA pairs for Tartary buckwheat seed development

Expression profiles of all the DEMs target genes were analyzed by using transcriptomic data from our previous report [46]. Then, the significantly different expression target genes were identified and functional analysis was performed by GO, KEGG and homologous annotation. Finally, expression correlation analysis between DEMs and differently expressed target genes was performed. The miRNAs and target mRNA pairs with R ≥ 0.5 and P value < 0.05 were identified, and the key miRNA-mRNA pairs for Tartary buckwheat seed development were obtained by homologous annotation of target mRNAs.

qRT-PCR validation of DEMs and its differently expressed target genes

Total RNA was isolated from three different developmental stages of Tartary buckwheat seeds by using the EASYspin Plus Plant RNA Kit (Aidlab, Beijing, China) and digested by DNase I (TaKaRa, Dalian, China) to remove the genomic DNA. Then, 2 μg total RNA was reverse transcribed into single-stranded cDNA using the miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China) following the manufacturer’s instructions. The qRT-PCR was performed on a CFX96 Real-time System (BIO-RAD, USA) using the miRcute Plus miRNA qPCR Kit (SYBR Green) (Tiangen, Beijing, China). Three independent biological replicates were used, and the FtU6 was used as the reference gene. The relative expression level of each miRNA was calculated by using the 2−ΔΔCt method. All primers used in this experiment are listed in Additional file 1: Table S12.

For differently expressed target genes, reverse transcription and qRT-PCR reactions were carried out by using PrimeScript™ RT Master Mix (Perfect Real Time) (TaKaRa, Beijing, China) and TB Green® Premix Ex Taq™ (Tli RNaseH Plus) (TaKaRa, Beijing, China), respectively. FtActin was used as the reference gene. All samples were performed three independent biological replicates. The relative expression level of each gene was calculated by using the 2−ΔΔCt method. The primers used in this experiment are listed in Additional file 1: Table S12.

Validation of the miRNA-directed cleavage of their predicted targets

The 5′-RLM-RACE method was used to investigate the miRNA-directed cleavage of their predicted target mRNA in vitro. Total RNA was ligated to an RNA adaptor and was reverse transcribed using the FirstChoice® RLM-RACE Kit (Ambion, USA) according to the manufacturer’s instructions. The 5′-RLM-RACE was performed as described previously by DeBoer et al. [6]. The primers used in this experiment are listed in Additional file 1: Table S13.

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article and its additional files. The dataset and materials presented in the investigation is available by request from the corresponding author.

Abbreviations

DEMs:

Differentially expressed miRNAs

qRT-PCR:

Quantitative real-time polymerase chain reaction

5′-RLM-RACE:

Ligase-mediated rapid amplification of 5′ cDNA ends

miRNAs:

MicroRNAs

sRNAs:

Small RNAs

DCL1:

Dicer-Like 1

HYL1:

Hyponastic Leaves 1

SE:

Serrate

HEN1:

Hua enhancer 1

AGO1:

Argonaute 1

RISC:

RNA-induced silencing complex

DEGs:

Differentially expressed genes

TFs:

Transcription factors

hetsiRNAs:

Heterochromatic siRNAs

AMFE:

Minimum free energy

PAL:

Phenylalanine ammonia-lyase

4CL:

4-coumarate-CoA ligase

FLS:

Flavonol synthase

GO:

Gene ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

DAP:

Days after pollination

References

  1. 1.

    Filipowicz W, Bhattacharyya SN, Sonenberg N. Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight? Nat Rev Genet. 2008;9(2):102–14.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  2. 2.

    Voinnet O. Origin, biogenesis, and activity of plant microRNAs. Cell. 2009;136(4):669–87.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  3. 3.

    He X, Shenkute AG, Wang W, Xu S. Characterization of conserved and novel microRNAs in Lilium lancifolium Thunb. by high-throughput sequencing. Sci Rep. 2018;8(1):2880.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  4. 4.

    Li H, Sun X, Ma C, Meng H, Chen Q. Identification and characterization of microRNAs during common buckwheat (Fagopyrum esculentum) seed development. Int J Agric Biol. 2020;24:1115–24.

    CAS  Google Scholar 

  5. 5.

    Yu Y, Jia T, Chen X. The ‘how’ and ‘where’ of plant microRNAs. New Phytol. 2017;216(4):1002–17.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    DeBoer K, Melser S, Sperschneider J, Kamphuis LG, Garg G, Gao LL, Frick K, Singh KB. Identification and profiling of narrow-leafed lupin (Lupinus angustifolius) microRNAs during seed development. BMC Genomics. 2019;20(1):135.

    PubMed  PubMed Central  Article  Google Scholar 

  7. 7.

    Shuai P, Liang D, Zhang Z, Yin W, Xia X. Identification of drought-responsive and novel Populus trichocarpa microRNAs by high-throughput sequencing and their targets using degradome analysis. BMC Genomics. 2013;14:233.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  8. 8.

    Park MY, Wu G, Gonzalez-Sulser A, Vaucheret H, Poethig RS. Nuclear processing and export of microRNAs in Arabidopsis. Proc Natl Acad Sci U S A. 2005;102(10):3691–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. 9.

    Li B, Qin Y, Duan H, Yin W, Xia X. Genome-wide characterization of new and drought stress responsive microRNAs in Populus euphratica. J Exp Bot. 2011;62(11):3765–79.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  10. 10.

    D'Ario M, Griffiths-Jones S, Kim M. Small RNAs: big impact on plant development. Trends Plant Sci. 2017;22(12):1056–68.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Sheng J, Zheng X, Wang J, Liu S, Pang M, Zhou F, et al. MicroRNAs and their targets by high-throughput sequencing and differentially expressed analysis in five Miscanthus species. Ind Crop Prod. 2018;118:196–205.

    CAS  Article  Google Scholar 

  12. 12.

    Li JB, Ding J, Yu X, Li H, Ruan CJ. Identification and expression analysis of critical microRNA-transcription factor regulatory modules related to seed development and oil accumulation in developing Hippophae rhamnoides seeds. Ind Crop Prod. 2019;137:33–42.

    CAS  Article  Google Scholar 

  13. 13.

    Chen CH, Zhong YD, Yu FX, Xu M. Deep sequencing identifies miRNAs and their target genes involved in the biosynthesis of terpenoids in Cinnamomum camphora. Ind Crop Prod. 2020;145:111853.

    CAS  Article  Google Scholar 

  14. 14.

    Yang Y, Guo J, Cheng J, Jiang Z, Xu N, An X, et al. Identification of UV-B radiation responsive microRNAs and their target genes in chrysanthemum (Chrysanthemum morifolium Ramat) using high-throughput sequencing. Ind Crop Prod. 2020;151:112484.

    CAS  Article  Google Scholar 

  15. 15.

    Ye YJ, Wang JW, Ni ZX, Meng X, Feng YH, Yang ZQ, et al. Small RNA and degradome sequencing reveal roles of miRNAs in strobilus development in masson pine (Pinus massoniana). Ind Crop Prod. 2020;154:112724.

    CAS  Article  Google Scholar 

  16. 16.

    Li C, Zhang BH. MicroRNAs in control of plant development. J Cell Physiol. 2016;231:303–13.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  17. 17.

    Savadi S. Molecular regulation of seed development and strategies for engineering seed size in crop plants. Plant Growth Regul. 2018;84:401–22.

    CAS  Article  Google Scholar 

  18. 18.

    Ma X, Zhang X, Zhao K, Li F, Li K, Ning L, et al. Small RNA and degradome deep sequencing reveals the roles of microRNAs in seed expansion in peanut (Arachis hypogaea L.). Front Plant Sci. 2018;9:349.

    PubMed  PubMed Central  Article  Google Scholar 

  19. 19.

    Zhu QH, Spriggs A, Matthew L, Fan L, Kennedy G, Gubler F, et al. A diverse set of microRNAs and microRNA-like small RNAs in developing rice grains. Genome Res. 2008;18(9):1456–65.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Peng T, Sun H, Qiao M, Zhao Y, Du Y, Zhang J, et al. Differentially expressed microRNA cohorts in seed development may contribute to poor grain filling of inferior spikelets in rice. BMC Plant Biol. 2014;14:196.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  21. 21.

    Jin X, Fu Z, Lv P, Peng Q, Ding D, Li W, et al. Identification and characterization of microRNAs during maize grain filling. PLoS One. 2015;10(5):e0125800.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  22. 22.

    Li D, Liu Z, Gao L, Wang L, Gao M, Jiao Z, et al. Genome-wide identification and characterization of microRNAs in developing grains of Zea mays L. PLoS One. 2016;11(4):e0153168.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  23. 23.

    Meng F, Liu H, Wang K, Liu L, Wang S, Zhao Y, et al. Development-associated microRNAs in grains of wheat (Triticum aestivum L.). BMC Plant Biol. 2013;13:140.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  24. 24.

    Hou G, Du C, Gao H, Liu S, Sun W, Lu H, et al. Identification of microRNAs in developing wheat grain that are potentially involved in regulating grain characteristics and the response to nitrogen levels. BMC Plant Biol. 2020;20(1):87.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  25. 25.

    Curaba J, Spriggs A, Taylor J, Li Z. Helliwell C miRNA regulation in the early development of barley seed. BMC Plant Biol. 2012;12:120.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Bai B, Shi B, Hou N, Cao Y, Meng Y, Bian H, et al. microRNAs participate in gene expression regulation and phytohormone cross-talk in barley embryo during seed development and germination. BMC Plant Biol. 2017;17(1):150.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  27. 27.

    Song QX, Liu YF, Hu XY, Zhang WK, Ma B, Chen SY, et al. Identification of miRNAs and their target genes in developing soybean seeds by deep sequencing. BMC Plant Biol. 2011;11:5.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. 28.

    Yu L, Guo R, Jiang Y, Ye X, Yang Z, Meng Y, et al. Genome-wide identification and characterization of novel microRNAs in seed development of soybean. Biosci Biotechnol Biochem. 2019;83(2):233–42.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  29. 29.

    Wang J, Jian H, Wang T, Wei L, Li J, Li C, et al. Identification of microRNAs actively involved in fatty acid biosynthesis in developing Brassica napus seeds using high-throughput sequencing. Front Plant Sci. 2016;7:1570.

    PubMed  PubMed Central  Google Scholar 

  30. 30.

    Wei W, Li G, Jiang X, Wang Y, Ma Z, Niu Z, et al. Small RNA and degradome profiling involved in seed development and oil synthesis of Brassica napus. PLoS One. 2018;13(10):e0204998.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  31. 31.

    Zhang YC, Yu Y, Wang CY, Li ZY, Liu Q, Xu J, et al. Overexpression of microRNA OsmiR397 improves rice yield by increasing grain size and promoting panicle branching. Nat Biotechnol. 2013;31(9):848–52.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  32. 32.

    Huang J, Li Z, Zhao D. Deregulation of the OsmiR160 target gene OsARF18 causes growth and developmental defects with an alteration of auxin signaling in rice. Sci Rep. 2016;6:29938.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  33. 33.

    Zhang JP, Yu Y, Feng YZ, Zhou YF, Zhang F, Yang YW, et al. MiR408 regulates grain yield and photosynthesis via a phytocyanin protein. Plant Physiol. 2017;175(3):1175–85.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. 34.

    Zhang H, Zhang J, Yan J, Gou F, Mao Y, Tang G, et al. Short tandem target mimic rice lines uncover functions of miRNAs in regulating important agronomic traits. Proc Natl Acad Sci U S A. 2017;114(20):5277–82.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  35. 35.

    Peng T, Qiao M, Liu H, Teotia S, Zhang Z, Zhao Y, et al. A resource for inactivation of microRNAs using short tandem target mimic technology in model and crop plants. Mol Plant. 2018;11(11):1400–17.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  36. 36.

    Jiao Y, Wang Y, Xue D, Wang J, Yan M, Liu G, et al. Regulation of OsSPL14 by OsmiR156 defines ideal plant architecture in rice. Nat Genet. 2010;42(6):541–4.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  37. 37.

    Duan P, Ni S, Wang J, Zhang B, Xu R, Wang Y, et al. Regulation of OsGRF4 by OsmiR396 controls grain size and yield in rice. Nat Plants. 2015;2:15203.

    PubMed  Article  CAS  Google Scholar 

  38. 38.

    Li S, Gao F, Xie K, Zeng X, Cao Y, Zeng J, et al. The OsmiR396c-OsGRF4-OsGIF1 regulatory module determines grain size and yield in rice. Plant Biotechnol J. 2016;14(11):2134–46.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. 39.

    Zhao YF, Peng T, Sun HZ, Teotia S, Wen HL, Du YX, et al. miR1432-OsACOT (acyl-CoA thioesterase) module determines grain yield via enhancing grain filling rate in rice. Plant Biotechnol J. 2019;17(4):712–23.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  40. 40.

    Miao C, Wang D, He R, Liu S, Zhu JK. Mutations in MIR396e and MIR396f increase grain size and modulate shoot architecture in rice. Plant Biotechnol J. 2020;18(2):491–501.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  41. 41.

    Sun W, Xu XH, Li Y, Xie L, He Y, Li W, et al. OsmiR530 acts downstream of OsPIL15 to regulate grain yield in rice. New Phytol. 2020;226(3):823–37.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  42. 42.

    Pan J, Huang D, Guo Z, Kuang Z, Zhang H, Xie X, et al. Overexpression of microRNA408 enhances photosynthesis, growth, and seed yield in diverse plants. J Integr Plant Biol. 2018;60(4):323–40.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  43. 43.

    Liu M, Ma Z, Zheng T, Sun W, Zhang Y, Jin W, et al. Insights into the correlation between physiological changes in and seed development of Tartary buckwheat (Fagopyrum tataricum Gaertn.). BMC Genomics. 2018;19(1):648.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  44. 44.

    Gao LC, Wang YF, Zhu Z, Chen H, Li DD. Egmir5179 from the mesocarp of oil palm (Elaeis guineensis jacq.) regulates oil accumulation by targeting nad transporter 1. Ind Crop Prod. 2019;137:126–36.

    CAS  Article  Google Scholar 

  45. 45.

    Na G, Mu X, Grabowski P, Schmutz J, Lu C. Enhancing microRNA167A expression in seed decreases the α-linolenic acid content and increases seed size in Camelina sativa. Plant J. 2019;98(2):346–58.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  46. 46.

    Li H, Lv Q, Ma C, Qu J, Cai F, Deng J, et al. Metabolite profiling and transcriptome analyses provide insights into the flavonoid biosynthesis in the developing seed of Tartary buckwheat (Fagopyrum tataricum). J Agric Food Chem. 2019;67(40):11262–76.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  47. 47.

    Gao J, Wang T, Liu M, Liu J, Zhang Z. Transcriptome analysis of filling stage seeds among three buckwheat species with emphasis on rutin accumulation. PLoS One. 2017;12(12):e0189672.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  48. 48.

    Huang J, Deng J, Shi T, Chen Q, Liang C, Meng Z, et al. Global transcriptome analysis and identification of genes involved in nutrients accumulation during seed development of rice tartary buckwheat (Fagopyrum Tararicum). Sci Rep. 2017;7(1):11792.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  49. 49.

    Qian Y, Cheng Y, Cheng X, Jiang H, Zhu S, Cheng B. Identification and characterization of Dicer-like, Argonaute and RNA-dependent RNA polymerase gene families in maize. Plant Cell Rep. 2011;30(7):1347–63.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  50. 50.

    Zhang H, Xia R, Meyers BC, Walbot V. Evolution, functions, and mysteries of plant ARGONAUTE proteins. Curr Opin Plant Biol. 2015;27:84–90.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. 51.

    Ma J, Zhao P, Liu S, Yang Q, Guo H. The control of developmental phase transitions by microRNAs and their targets in seed plants. Int J Mol Sci. 2020;21(6):1971.

    CAS  PubMed Central  Article  Google Scholar 

  52. 52.

    Kapoor M, Arora R, Lama T, Nijhawan A, Khurana JP, Tyagi AK, et al. Genome-wide identification, organization and phylogenetic analysis of Dicer-like, Argonaute and RNA-dependent RNA Polymerase gene families and their expression analysis during reproductive development and stress in rice. BMC Genomics. 2008;9:451.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  53. 53.

    Yadav CB, Muthamilarasan M, Pandey G, Prasad M. Identification, characterization and expression profiling of Dicer-like, Argonaute and RNA-dependent RNA polymerase gene families in Foxtail Millet. Plant Mol Biol Report. 2015;33(1):43–55.

    CAS  Article  Google Scholar 

  54. 54.

    Zhong J, He W, Peng Z, Zhang H, Li F, Yao J. A putative AGO protein, OsAGO17, positively regulates grain size and grain weight through OsmiR397b in rice. Plant Biotechnol J. 2020;18(4):916–28.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  55. 55.

    Pachamuthu K, Swetha C, Basu D, Das S, Singh I, Sundar VH, et al. Rice-specific Argonaute 17 controls reproductive growth and yield-associated phenotypes. Plant Mol Biol. 2021;105(1–2):99–114.

    PubMed  Article  PubMed Central  Google Scholar 

  56. 56.

    Rodrigues AS, Miguel CM. The pivotal role of small non-coding RNAs in the regulation of seed development. Plant Cell Rep. 2017;36(5):653–67.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  57. 57.

    Xie M, Yu B. siRNA-directed DNA methylation in plants. Curr Genomics. 2015;16(1):23–31.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  58. 58.

    Zhang L, Li X, Ma B, Gao Q, Du H, Han Y, et al. The Tartary buckwheat genome provides insights into rutin biosynthesis and abiotic stress tolerance. Mol Plant. 2017;10(9):1224–37.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  59. 59.

    Csukasi F, Donaire L, Casañal A, Martínez-Priego L, Botella MA, Medina-Escobar N, et al. Two strawberry miR159 family members display developmental-specific expression patterns in the fruit receptacle and cooperatively regulate Fa-GAMYB. New Phytol. 2012;195(1):47–57.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  60. 60.

    Sharma D, Tiwari M, Pandey A, Bhatia C, Sharma A, Trivedi PK. MicroRNA858 is a potential regulator of phenylpropanoid pathway and plant development. Plant Physiol. 2016;171(2):944–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  61. 61.

    Li N, Xu R, Li Y. Molecular networks of seed size control in plants. Annu Rev Plant Biol. 2019;70:435–63.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  62. 62.

    Gupta OP, Karkute SG, Banerjee S, Meena NL, Dahuja A. Contemporary understanding of miRNA-based regulation of secondary metabolites biosynthesis in plants. Front Plant Sci. 2017;8:374.

    PubMed  PubMed Central  Google Scholar 

  63. 63.

    Nesi N, Debeaujon I, Jond C, Pelletier G, Caboche M, Lepiniec L. The TT8 gene encodes a basic helix-loop-helix domain protein required for expression of DFR and BAN genes in Arabidopsis siliques. Plant Cell. 2000;12(10):1863–78.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  64. 64.

    Nesi N, Jond C, Debeaujon I, Caboche M, Lepiniec L. The Arabidopsis TT2 gene encodes an R2R3 MYB domain protein that acts as a key determinant for proanthocyanidin accumulation in developing seed. Plant Cell. 2001;13(9):2099–114.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  65. 65.

    Mehrtens F, Kranz H, Bednarek P, Weisshaar B. The Arabidopsis transcription factor MYB12 is a flavonol-specific regulator of phenylpropanoid biosynthesis. Plant Physiol. 2005;138(2):1083–96.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  66. 66.

    Teng S, Keurentjes J, Bentsink L, Koornneef M, Smeekens S. Sucrose-specific induction of anthocyanin biosynthesis in Arabidopsis requires the MYB75/PAP1 gene. Plant Physiol. 2005;139(4):1840–52.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  67. 67.

    Tan J, Wu Y, Guo J, Li H, Zhu L, Chen R, et al. A combined microRNA and transcriptome analyses illuminates the resistance response of rice against brown planthopper. BMC Genomics. 2020;21(1):144.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  68. 68.

    Wang S, Wu K, Yuan Q, Liu X, Liu Z, Lin X, et al. Control of grain size, shape and quality by OsSPL16 in rice. Nat Genet. 2012;44(8):950–4.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  69. 69.

    Wang S, Li S, Liu Q, Wu K, Zhang J, Wang S, et al. The OsSPL16-GW7 regulatory module determines grain shape and simultaneously improves rice yield and grain quality. Nat Genet. 2015;47(8):949–54.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  70. 70.

    Dwivedi N, Maji S, Waseem M, Thakur P, Kumar V, Parida SK, et al. The mediator subunit OsMED15a is a transcriptional co-regulator of seed size/weight-modulating genes in rice. Biochim Biophys Acta Gene Regul Mech. 2019;1862(10):194432.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  71. 71.

    Fan C, Wang G, Wang Y, Zhang R, Wang Y, Feng S, et al. Sucrose synthase enhances hull size and grain weight by regulating cell division and starch accumulation in transgenic rice. Int J Mol Sci. 2019;20(20):4971.

    CAS  PubMed Central  Article  Google Scholar 

  72. 72.

    He Z, Zeng J, Ren Y, Chen D, Li W, Gao F, et al. OsGIF1 positively regulates the sizes of stems, leaves, and grains in rice. Front Plant Sci. 2017;8:1730.

    PubMed  PubMed Central  Article  Google Scholar 

  73. 73.

    Jain M, Pathak BP, Harmon AC, Tillman BL, Gallo M. Calcium dependent protein kinase (CDPK) expression during fruit development in cultivated peanut (Arachis hypogaea) under Ca2+-sufficient and -deficient growth regimens. J Plant Physiol. 2011;168(18):2272–7.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  74. 74.

    Li H, Lv Q, Deng J, Huang J, Cai F, Liang C, et al. Transcriptome analysis reveals key seed-development genes in common buckwheat (Fagopyrum esculentum). Int J Mol Sci. 2019;20(17):4303.

    CAS  PubMed Central  Article  Google Scholar 

  75. 75.

    Finn RD, Mistry J, Schuster-Böckler B, Griffiths-Jones S, Hollich V, Lassmann T, et al. Pfam: clans, web tools and services. Nucleic Acids Res. 2006;34(Database issue):D247–51.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  76. 76.

    Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, et al. TBtools: An integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  77. 77.

    Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47(D1):D155–62.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  78. 78.

    Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  79. 79.

    Allen E, Xie Z, Gustafson AM, Carrington JC. microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell. 2005;121(2):207–21.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  80. 80.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. Gene Ontol Cons Nat Genet. 2000;25(1):25–9.

    CAS  Google Scholar 

  81. 81.

    Okuda S, Yamada T, Hamajima M, Itoh M, Katayama T, Bork P, Goto S, Kanehisa M. KEGG atlas mapping for global analysis of metabolic pathways. Nucleic Acids Res. 2008;36:W423–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was funded by the National Key R&D Program of China (2019YFD1001300, 2019YFD1001302), the National Natural Science Foundation of China (31701494, 31760430, 31860408), the Science and Technology Foundation of Guizhou Province (QianKeHeJiChu [2019]1235), the National Natural Science Foundation of China-Guizhou Provincial Government Joint Program (U1812401), Guizhou provincial department of education youth science and technology talent growth project (Qianjiaohe KY Zi [2018]128), the Training Program of Guizhou Normal University (QianKeHePingTaiRenCai [2017]5726), and the Initial Fund for Doctor Research in Guizhou Normal University (11904/0517051). We thank all the foundation of economic support. The funding organizations provided the financial support to the research projects, but were not involved in the design of the study, data collection, analysis of the data, or the writing of the manuscript.

Author information

Affiliations

Authors

Contributions

HYL and QFC designed the research. HYL carried out the experiments and wrote the manuscript. HYL, HLM, XQS, JD, and LWZ performed data analysis. QYL, TXS and QFC reviewed and revised the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Hongyou Li or Qingfu Chen.

Ethics declarations

Ethics approval, guidelines and consent to participate

The plant material (Jinqiao No. 2) was a Tartary buckwheat cultivar in China, which was collected by our lab (Research Center of Buckwheat Industry Technology of Guizhou Normal University, Guiyang, Guizhou, China). All authors stated that this study comply with the Chinese legislation and the Convention on the Trade in Endangered Species of Wild Fauna and Flora.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Characteristic features of all the identified DCL, AGO, HYL1, SE, HST, and RDRs in tartary buckwheat genome. Table S2. Length distribution of sRNA sequences identified in developing tartary buckwheat seeds. Table S3. Information of identified conserved miRNAs. Table S4. Information of identified novel miRNAs. Table S5. Predicted target genes of miRNAs in tartary buckwheat. Table S6. Annotation of miRNAs target genes in tartary buckwheat. Table S7. Information of TFs targeting by miRNAs. Table S8. miRNAs have target genes are the orthologs of the known seed or organ size. Table S9. miRNAs have target genes are the orthologs of the known structural or regulatory genes of flavonoid biosynthesis. Table S10. KEGG pathways of the target genes of DEMs. Table S11. miRNA-mRNA interaction pairs show expression negative correlation during tartary buckwheat seed development. Table S12. Primers of sequences for qRT-PCR analysis. Table S13. Primers of sequences for RLM-5′RACE analysis.

Additional file 2: Figure S1.

Read length distribution of sRNAs. Figure S2. GO analysis of the target genes of DEMs. Figure S3. KEGG analysis of the target genes of DEMs.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, H., Meng, H., Sun, X. et al. Integrated microRNA and transcriptome profiling reveal key miRNA-mRNA interaction pairs associated with seed development in Tartary buckwheat (Fagopyrum tataricum). BMC Plant Biol 21, 132 (2021). https://doi.org/10.1186/s12870-021-02914-w

Download citation

Keywords

  • Tartary buckwheat
  • Integrated analysis
  • miRNAs
  • Target gene
  • Seed development
  • Seed size
  • Flavonoids
  • DCL
  • AGO
  • RDR