Sequencing of sRNAs and the transcriptome
Six sRNA libraries were generated to identify miRNAs involved in the fruit development of the pitaya. The Illumina sequencing data of sRNAs from the 19th, 25th and 29th day after flowering (DAF) showed that 24 nt sRNAs are the most abundant, followed by 21 nt sRNAs (Figure S1). A total of 82,318,241 reads were obtained from the sRNA datasets. After removal of the adaptor, insert, polyA and short RNAs of <18 nt in length, 62,729,725 (76.20%) valid reads were obtained, including the rRNA, tRNA, snRNA, snoRNA and some other Rfam RNA (Table S1).
Three transcriptome libraries (Hp19d, Hp25d and Hp29d) of pitaya pulps were constructed to identify the target genes of miRNAs. All screened reads were de novo assembled into 68,505 transcripts with an N50 of 1,700 bp. And then, these transcripts were assembled into 39,737 unique sequences with an average length of 879 bp. The length distributions of those transcripts and unigenes were shown in Figure S2. The distribution of assembled transcripts and genes with different GC contents in the transcriptome datasets were summarized in Figure S3. The majority of transcripts and genes were in the range of 35%–50% in GC contents.
Identification of miRNAs
A total of 126 conserved miRNAs were identified based on BLAST searches and sequence analyses in comparison of sRNA sequences with known mature plant miRNAs in miRBase (Table S2). Among them, 95 known miRNAs belonging to 53 families were obtained. The number of miRNA members for each family varied from 1 (for MIR535) to 7 (for MIR482) (Figure S4). Novel miRNAs were predicted using the sRNA valid reads based on structure and expression criteria [44]. 41 novel candidate miRNAs with a clear precursor including stem-loop secondary structure were identified (Table S2). The length of novel miRNAs was between 20-25 nt. The most abundant sequences (49%) were 21 nt-length, and followed by 24 nt-length (29%), which was consistent with typical length distribution of mature miRNAs. The higher expression levels of the miRNA, the more copies of the miRNAs were sequenced. The abundance of known miRNAs was higher than that of putative novel miRNAs, except Hmo-novel-2, Hmo-novel-6, Hmo-novel-21 and Hmo-novel-22, which possessed more than 100 normalized reads (Table S2).
Analyses of Differentially Expressed miRNAs
High-throughput sequencing (HTS) was performed to explore the expression changes of miRNAs involved in betalain biosynthesis of ‘Guanhuahong’ pitaya. Only miRNAs with expression values at p<0.05 were considered to be significantly regulated. In Hp29d/Hp25d/Hp19d, 33 known miRNAs and five novel miRNAs were found to be differentially expressed (Table 1).
Changes in expression levels of miRNAs during fruit development of pitaya possibly reflected their potential functional differences. In different pulp coloration stages of ‘Guanhuahong’ pitaya, the expression levels of Hmo-miR396b, Hmo-miR535, Hmo-novel-7 and Hmo-novel-12 at the red pulp stages (25th and 29th DAP) were significantly higher than that of the white pulp stage (19th DAP) (Table 1). Expression levels of Hmo-miR156, Hmo-miR157b, Hmo-miR164a, Hmo-miR164b and Hmo-miR399a increased gradually from 19th DAF to 29th DAF in ‘Guanhuahong’ pitaya. Those results were in consistent with betalain accumulation pattern [45, 46]. However, the expression levels of Hmo-miR171c, Hmo-miR408, Hmo-miR393, Hmo-miR398c, Hmo-miR8175, Hmo-miR398a, Hmo-miR168a, Hmo-miR529b, Hmo-miR398b, Hmo-miR172a, Hmo-miR159a and Hmo-miR397b decreased gradually from 19th DAF to 29th DAF during fruit maturation of ‘Guanhuahong’ pitaya. Twelve differentially expressed miRNAs i.e. Hmo-miR390a, Hmo-miR6149a, Hmo-miR159b, Hmo-miR530, Hmo-miR6300, Hmo-miR390b, Hmo-miR6020, Hmo-miR394, Hmo-miR5072, Hmo-novel-15, Hmo-miR171d and Hmo-miR482f showed a lower expression level at the red pulp stages (25th and 29th DAP) compared to a relatively higher expression level at white pulp stage (19th DAP). Those results were in contrast to betalain accumulation pattern suggesting that those miRNAs negatively regulated their target genes in betalain accumulation.
Bioinformatics of transcriptome analyses
All 39,737 unique sequences were annotated according to BLAST (cut-off E-value ≤10−5) searches of Nr, Swiss-Prot Protein, Pfam, Gene Ontology (GO), euKaryotic Ortholog Groups (KOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Table S3).
The GO, KEGG and KOG databases were used to classify the functions of the predicted unigenes. 12,559 unigenes were classified into three main categories: ‘biological process’, ‘cellular component’, and ‘molecular function’ by GO database (Fig. 1). As for the ‘molecular function’ category, the largest number of unigenes was gathered in ‘ATP binding’, while the major groups for the ‘cellular component’ category were ‘integral to membrane’ (2,878 unigenes, 23%), ‘nucleus’ (2,211 unigenes, 18%) and ‘plasma membrane’ (1,646 unigenes, 13%). In the category of ‘biological process’, ‘regulation of transcription, DNA-dependent’ (729 unigenes, 5.8%) and ‘transcription, DNA-dependent’ (706 unigenes, 5.6%) were the two most abundant subcategories.
7,822 unigenes were mapped onto KEGG pathways. The highest number of unigenes was carbohydrate metabolism (796 unigenes), followed by energy metabolism (596 unigenes), amino acid metabolism (595 unigenes) and translation (510 unigenes) (Fig. 2). Differential expression genes (DEGs) in the Hp25d-VS-Hp19d, Hp29d-VS-Hp19d, and Hp29d-VS-Hp25d were analyzed. In the Hp25d-VS-Hp19d, 24, 20, 16 and 13 DEGs were involved in phenylpropanoid, stilbenoid, diarylheptanoid and gingerol biosynthesis, DNA replication and flavonoid biosynthesis pathway (p<0.0001), respectively (Fig. 3A). In Hp29d-VS-Hp19d, 140 genes were found to be significantly differentially expressed (p<0.0001), including biosynthesis of phenylpropanoid, flavonoid, stilbenoid, diarylheptanoid and gingerol, metabolism of phenylalanine, α-linolenic acid, starch and sucrose, xenobiotics by Cyt P450, as well as drug metabolism-Cyt P450 (Fig. 3B). The expression levels of DEGs among different groups in one pathway showed different scales of changes. 13 DEGs were detected from the flavonoid biosynthesis pathway in the Hp25d-VS-Hp19d and Hp29d-VS-Hp19d, respectively. However, no DEGs related to flavonoid biosynthesis were found in Hp29d-VS-Hp25d (Fig. 3C). Those results suggested that flavonoid biosynthesis is a crucial pathway in pitaya betalain biosynthesis.
The KOG database is used to study the classification and evolutionary rates of orthologous proteins. As shown in Fig. 4, group R (general function prediction only), group O (posttranslational modification, protein turnover, chaperones) and group T (signal transduction mechanisms) are the three most abundant groups in pitaya dataset, suggesting that a large number of transcriptional and posttranslational regulation of gene expression and function are involved in pitaya fruit development.
DEGs from three development stages (Hp25d-VS-Hp19d, Hp29d-VS-Hp19d, and Hp29d-VS-Hp25d) were evaluated by pairwise comparisons using the expression fold (|log2fold change|≥1) and p values<0.05 as the thresholds (Figure 5). In the pairwise comparisons between any two stages, 3,988 genes were found to be significantly differentially expressed. The highest amount of DEGs was obtained between the Hp19d and Hp29d libraries, including 1,486 down-regulated and 733 up-regulated (Fig. 5A). The lowest number of DEGs (1835) was detected between the Hp25d and Hp29d libraries (987 down-regulated and 848 up-regulated), followed by the Hp19d and Hp25d libraries (1,961 DEGs inluding 1,274 down-regulated and 687 up-regulated). Among those DEGs, 124 genes were significantly differentially expressed in all three fruit development stages of pitaya (Fig. 5B).
Verification of the accuracy of the sRNAome and RNA-Seq data
To verify the reliability of the sRNAome and RNA-Seq results, the expression of fifteen miRNAs and fifteen genes related to betalain biosynthesis were analyzed by RT-qPCR. The IDs, reads per kilobase of exon model per million mapped read (RPKM) value, and primers for the 30 transcripts were shown in Table S2, S5, S6, and S8, respectively. The overall correlation coefficients of 0.717** and 0.719** were obtained by linear regression [(Q-PCR value) = a (sRNAome or RNA-Seq value) + b] analysis (Fig. 6), respectively, indicating that the results of sRNAome and transcriptome were consistent with those of RT-qPCR. Those results suggested that sRNAome and RNA-Seq data can be used for subsequent experiments.
Prediction of targets for differentially expressed miRNAs
To investigate the functions of miRNAs in betalain biosynthesis of pitaya, it is crucial to predict their target genes. Target genes of differentially-expressed miRNAs were predicted by Target Finder software. Transcripts of transcriptome were identified as possible target genes for the majority of miRNAs. A total of 124 target mRNAs, 22 target mRNAs for 25 conserved miRNAs and 5 novel miRNAs were obtained, respectively (Table S9). More than one target gene was predicted for most miRNAs. miRNAs have the potential to regulate targets belonging to certain gene families with different biological functions. For example, the predicted target genes of Hmo-miR529b were found to be involved in squamosa promoter-binding-like protein 6 (SPL6), histone-lysine N-methyltransferase and ubiquitin-protein ligase. Based on transcriptome annotation, six targets involved in betalain biosynthesis were selected. Among them, HmCYP83B1-like, HmCYP71A8-like, HmTPST-like, HmWDTC1-like and HmSPL16-like were respectively targeted by Hmo-novel-2, Hmo-miR6020, Hmo-novel-15, Hmo-miR164a and Hmo-miR156 while HmSPL6-like was co-targeted by Hmo-miR156, Hmo-miR157b and Hmo-miR529b.
Tissue-specific analyses of differentially expressed miRNAs
miRNA preferentially expressed in specific tissues can provide clues to its physiological functions. Thirty differentially expressed miRNAs were analyzed their functions by RT-qPCR using ten different tissues from ‘Guanhuabai’ and ‘Guanhuahong’ pitayas. As shown in Figure 7, the 30 differentially expression miRNAs showed different expression levels in those pitaya tissues. Hmo-miR172a, Hmo-miR394, Hmo-miR530, Hmo-miR6020, Hmo-miR397b, Hmo-miR160b and Hmo-miR398b displayed similar expression patterns. Hmo-miR156 preferentially expressed in pitaya roots and fruits. Hmo-miR398a strongly expressed in pitaya stamens compared to weak expression in the other tissues. Higher expression levels of Hmo-novel-7 and Hmo-miR171d were detected in petals but moderately or weakly expressed in the other tissues. Hmo-miR164a and Hmo-miR164b displayed constitutive expression and the highest expression level was detected in pitaya fruits. Hmo-miR396b and Hmo-miR5072 had higher expression in petals of ‘Guanhuabai’ pitaya compared to higher expression of Hmo-miR171c and Hmo-novel-2 in petals of ‘Guanhuahong’ pitaya. Hmo-novel-12, Hmo-novel-21, Hmo-novel-15, Hmo-miR393, Hmo-miR390b, Hmo-miR159a and Hmo-miR408 were preferentially expressed in petals and receptacles of ‘Guanhuahong’ pitaya. Hmo-miR157b and Hmo-miR6300 were highly expressed in ‘Guanhuahong’ and ‘Guanhuabai’ pitayas, respectively. Results from RT-qPCR indicated that expression of the 30 miRNAs from pitaya had tissue- and/or cultivar- characteristics.
Validation of miRNAs and target genes related to betalain biosynthesis
miRNAs involved in pitaya betalain biosynthesis were screened based on annotated transcripts and previous reports related to pigment synthesis. Fourteen target genes predicted from 11 miRNAs were assayed by RT-qPCR. The target genes from the ‘Guanhuahong’ and ‘Guanhuabai’ pitayas had the same sequences (Fig. S5). Sixteen targets showed decreased or increased expression trends along with the increased or decreased expression of the miRNA, suggesting that they might be actively cleaved by miRNAs (Fig. 8). For example, Hpcyt P450-like3 (Fig. 8A9 and B9) and Hpcyt P450-like2 (Fig. 8A10 and B10) targeted by Hmo-miR160a, HmMYB12-like (Fig. 8A5), HmMYBC1-like (Fig. 8A6) and HmMYB2-like (Fig. 8A7) targeted by Hmo-miR858, and HmCYP71A8-like (Fig. 8A2) targeted by Hmo-miR6020 showed increasing expression trend at first and decreasing thereafter while Hmo-miR160a, Hmo-miR858 and Hmo-miR6020 showed a reverse trend at all pulp coloration stages of pitaya. HmSPL6-like (Fig. 8A15) targeted by Hmo-miR156 showed decreasing while Hmo-miR156 showed increasing at all pulp coloration stages of pitaya. However, some targeted genes and their miRNAs such as HmTPST-like (Fig. 8B11) targeted by Hmo-novel-15 in ‘Guanhuahong’ pitaya and HmSPL6-like (Fig. 8A13) targeted by Hmo-miR529b in ‘Guanhuabai’ pitaya had similar expression patterns at all pulp stages of pitaya.
5′RACE analyses were further verified the fourteen candidate targets. The HmSPL6-like, HmTT2-like, HmMYB12-like, HmMYBC1-like, Hpcyt P450-like3, HmCYP83B1-like and HmTPST-like were confirmed to be cleaved by their corresponding miRNAs (Fig. 9). The cleavage sites of Hmo-novel-2 on HmCYP83B1-like and Hmo-novel-15 on HmTPST-like were both occurred at the 10th nucleotide from the 5′-end of miRNAs in the binding region. The cleavage frequency of Hmo-novel-2 on HmCYP83B1-like and Hmo-novel-15 on HmTPST-like was up to 10/10 both in ‘Guanhuabai’ and ‘Guanhuahong’ pitayas, respectively. Those results confirmed that Hmo-novel-2 and Hmo-novel-15 can guide the cleavage of the mRNA of HmCYP83B1-like and HmTPST-like, respectively. Cleavage occurred mostly at the 10th nucleotide from the 5′-end in the binding sites. However, HmMYBC1-like occurred at the 9th nucleotide in ‘Guanhuahong’ pitaya, which may be due to the wide range of miRNA cutting mRNA caused by siRNA interference. Besides, the same gene from ‘Guanhuahong’ pitaya is different from ‘Guanhuabai’ pitaya (Fig. 9), which may be responsible for the different pulp colors of the two pitaya cultivars.
The relationship between miRNAs and their target genes
A total of 5 miRNAs and 7 corresponding target genes i.e. Hmo-miR160a-Hpcyt P450-like2, Hmo-miR6020-HmCYP71A8-like, Hmo-novel-2-HmCYP83B1-like, Hmo-novel-15-HmTPST-like, Hmo-miR858-HmMYB12-like, Hmo-miR858-HmMYBC1-like and Hmo-miR858-HmMYB2-like were verified their interactions in a tobacco transient expression system. The interactions could affect miRNA processing resulting in higher precursor accumulation and reduced mature miRNA in pitaya, and further regulate pitaya betalain biosynthesis (Fig. 10).