Integrated sRNAome and RNA-Seq analysis reveals miRNA effects on betalain biosynthesis in pitaya

Background MicroRNAs (miRNAs) and their regulatory functions in anthocyanin, carotenoid, and chlorophyll accumulation have been extensively characterized in many plant species. However, the miRNA regulatory mechanism in betalain biosynthesis remains mostly unknown. Results In this study, 126 conserved miRNAs and 41 novel miRNAs were first isolated from Hylocereus monacanthus, among which 95 conserved miRNAs belonged to 53 miRNA families. Thirty-four candidate miRNAs related to betalain biosynthesis were differentially expressed. The expression patterns of those differential expressed miRNAs were analyzed in various pitaya tissues by RT-qPCR. A significantly negative correlation was detected between the expression levels of half those miRNAs and corresponding target genes. Target genes of miRNAs i.e. Hmo-miR157b-HmSPL6-like, Hmo-miR160a-Hpcyt P450-like3, Hmo-miR6020-HmCYP71A8-like, Hmo-novel-2-HmCYP83B1-like, Hmo-novel-15-HmTPST-like, Hmo-miR828a-HmTT2-like, Hmo-miR858-HmMYB12-like, Hmo-miR858-HmMYBC1-like and Hmo-miR858-HmMYB2-like were verified by 5′RACE and transient expression system in tobacco. Conclusions Hmo-miR157b, Hmo-miR160a, Hmo-miR6020 Hmo-novel-2, Hmo-novel-15, Hmo-miR828a and Hmo-miR858 play important roles in pitaya fruit coloration and betalain accumulation. Our findings provide new insights into the roles of miRNAs and their target genes of regulatory functions involved in betalain biosynthesis of pitaya.


Background
Mature microRNAs (miRNAs) are a type of endogenous non-coding small RNAs with 20-24 nucleotide (nt) length. miRNAs regulate their target genes by mRNA at the posttranscriptional level via the RNA-induced silencing complex (RISC) by binding with the Argonaute (AGO) protein to cleavage target mRNA or repress translation of target mRNA [1]. miRNAs play vital roles in plant growth and development, (a)biotic stress response, post-transcriptional regulation, and pigment regulation [2][3][4].
Pitaya, also known as pitahaya or dragon fruit, is a perennial climbing fruit crop belonging to the genus Hylocereus (Cactaceae) under the order Caryophyllales. As a member of the Cactaceae, pitaya exhibits a range of specific adaptations to arid lands in term of succulent stems with spines instead of leaves, the crassulacean acid metabolism (CAM) pathway [26][27][28]. It is an excellent plant material for basic and applied biological research. The potential economic impact of pitaya lies in its diverse uses not only as agricultural produce and processed foods but also in industrial and medicinal products. Pitaya is a fast-return fruit crop with production in the second year after planting and full production in 3-4 years. Pitaya fruit is mature in 28-50 days (28-35 days in summer and 35-50 days in autumn) after flowering and has 7-12 separate fruiting cycles per year due to climatic or nutritional limitations in South China. Therefore, pitaya has become a favorite fruit of many farmers and home gardeners in Southeast Asia, China, the United States, Israel, Australia, Cyprus and the Canary Islands.
The color of pitaya is attributed to the presence of betalains [29][30][31][32]. Betalains are red and yellow alkaloid pigments that are found in all families of the Caryophyllales with the exception of Molluginaceae and Caryophyllaceae which produce anthocyanins [33]. Betalains in pitaya fruit are not only good for human health but also can help consumers distinguish cultivars [29,34,35]. Betalains also play vital roles in the protection against drought, UV radiation, high saline soils, and diseases [35][36][37][38][39]. Pitaya is the only at large-scale commercially grown fruit containing abundant betalains for the consumer. Previous studies are mainly focused on characterizations of key genes and TFs involved in betalain biosynthesis. Key genes such as tyrosinase (TYR), cytochrome P450 (Cyt P450), 4,5-dihydroxy-phenylalanine (DOPA)-dioxygenase (DOD) and glucosyltransferases (GTs) [40] and TFs such as WRKY and MYB involved in betalain biosynthesis have been investigated in detail [41][42][43]. However, the roles of miRNAs in betalain biosynthesis has not been reported yet. In this study, candidate miRNAs and their target genes related to betalain biosynthesis were identified based on small RNA and transcriptome databases of pitaya pulp at different developmental stages. The aim of the present study is to explore the roles of miRNA in pitaya betalain biosynthesis, which may contribute to a better understanding of betalain biosynthesis in Hylocereus.

Sequencing of sRNAs and the transcriptome
Six sRNA libraries were generated to identify miRNAs using pulps from 'Guanhuahong' pitaya on the 19 th day after flowering (DAF) (white pulp stage, Hp19d_1, and Hp19d_2), 25 th DAF (pulp coloration stages, Hp25d_1, and Hp25d_2) and 29 th DAF (mature stage, Hp29d_1, and Hp29d_2). The Illumina sequencing data of sRNAs from the 19 th , 25 th and 29 th 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 were constructed to identify the target genes of miRNAs using pulps from 'Guanhuahong' pitaya on the 19 th DAF (white pulp stage, Hp19d), 25 th DAF (pulp coloration stages, Hp25d) and 29 th DAF (mature stage, Hp29d). All screened reads were de novo assembled into 68,505 transcripts with an N50 of 1700 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 miR-NAs 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]. Fortyone novel candidate miRNAs with a clear precursor including stem-loop secondary structure were identified (Table S2). The length of novel miRNAs was between 20 and 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).

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 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 (|log 2 fold change| ≥ 1) and p values< 0.05 as the thresholds (Fig. 5). In the pairwise comparisons between any two stages, 3988 genes were found to be significantly differentially expressed. The highest amount of DEGs was obtained between the Hp19d and Hp29d libraries, including 1486 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 (1961 DEGs including 1274 down-regulated and 687 up-regulated). Among those DEGs, 124 genes were significantly differentially expressed in all three fruit development stages of 'Guanhuahong' 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 Tables S2, S5, S6, and   (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 ubiquitinprotein ligase. Based on transcriptome annotation, six targets involved in betalain biosynthesis were selected [23,[40][41][42]. 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 Fig. 7, the The comparison of DEGs between any two stages of the pitaya pulp 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 ( Figure 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, HmMYB12like, 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 Hmonovel-2 on HmCYP83B1-like and Hmo-novel-15 on HmTPST-like were both occurred at the 10 th 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 Hmonovel-15 can guide the cleavage of the mRNA of HmCYP83B1-like and HmTPST-like, respectively. Cleavage occurred mostly at the 10 th nucleotide from the 5′-end in the binding sites. However, HmMYBC1-like occurred at the 9 th 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.

Discussion
Great progress has been made in pitaya betalain studies in term of physical and chemical properties [47,48], purification and identification [49][50][51][52][53][54], antioxidant and radical scavenging capacity [55][56][57] as well as metabolic and transcriptional analyses [45,46]. However, no information is available about miRNAs involved in betalain biosynthesis of pitaya. Therefore, identification of pitaya miRNAs associated with their target genes will help our understanding of molecular regulatory mechanisms of pitaya betalain biosynthesis. In this study, transcriptome and sRNAome were performed to explore the role of miRNAs in pitaya coloration mechanisms at different developmental stages. The percentage of 24 nt sRNAs (an average of 58.03%) was much higher than that of 21 nt sRNAs (an average of 12.16%) in pitaya pulp ( Figure  S1). The result was consistent with the length distribution of sRNAs in apple and litchi [15,20,58], but was inconsistent with the findings in strawberry, orange, Brassica juncea, and apple [21,[59][60][61] suggesting that different plant species have different length sRNA distribution. A total of 95 known miRNAs belonging to 53 miRNA families and 41 new miRNAs were identified from 'Guanhuahong' pitaya (H. monacanthus) (Table  S2). miRNAs showed different expression levels during fruit developmental stages of pitaya, indicating that miR-NAs play essential roles in pitaya fruit development. Those results suggested that miRNAs are involved in pitaya growth and development through different expressions in various pitaya tissues.  [13,20,21,23]. Bgy-mir156 regulates the target gene phytoene synthase (PSY) and affects the accumulation of carotenoids in carrot (Daucus carota) [12]. Betalains are secondary metabolites but cannot co-exist naturally in one plant at the same time [62]. Pitaya is a high-value, functional fruit containing abundant betalains. In our study, the highest expression level of Hmo-miR157b was detected on the 23 rd DAF (color conversion) in 'Guanhuahong' pitaya, suggesting that Hmo-miR157b may play significant roles in betalain accumulation. Results from 5′RACE showed that HmSPL6-like was targeted by Hmo-miR157b (Fig. 7). Those results indicated that Hmo-miR157b might positively regulate betalain biosynthesis by SPL TFs, while SPL TFs negatively regulates betalain accumulation in pitaya.
Hpcyt P450-like3 is involved in betalain biosynthesis in H. monacanthus [45]. In this study, Hpcyt P450-like3 was targeted by Hmo-miR160a, suggesting that Hmo-miR160a was involved in betalain biosynthesis (Fig. 9B5). Novel miRNAs are involved in accumulations of chlorophyll, carotenoid, and anthocyanin. Ttu-novel-48 could regulate chlorophyll accumulation in leaves of durum wheat [6]. Csi-novel-03 regulates carotenoid pathways by AP2 TFs in Citrus [8]. In strawberry receptacle fruit  Figure S6 ripening, ABA-induced Fa_novel23 resulted in the rapid accumulation of fruit anthocyanin [25]. In our study, two novel miRNAs: Hmo-novel-2 and Hmo-novel-15 were obtained. Hmo-novel-2 influenced betalain accumulation via the regulation of HmCYP83B1-like (Fig. 9A6-B6 and Fig. 10A11-C11 and A12-C12). Hmo-novel-15 was identified as a regulator of betalain biosynthesis regulating the expression of HmTPST-like (Figs. 9B7 and 10A13-C13 and A14-C14). Hmo-miR6020-HmCYP71A8-like was not verified in 5′RACE but in the tobacco transient expression system (Fig. 10A1-C1 and A2-C2) suggesting that the translation inhibited the regulation modes of miRNAs on their targets rather than the degradation of mRNAs. The rest of the differentially expressed miRNAs and their targets in HTS were not confirmed in 5′RACE and tobacco transient expression system, indicating that these genes might not be targets of these miRNAs. Further work is necessary to elucidate their roles in betalain biosynthesis of pitaya.

Conclusions
In this study, sRNAome and RNA-Seq were first used to identify differentially expressed miRNAs and their target genes involved in betalain biosynthesis. Comprehensive sRNAome analyses uncovered 126 conserved miRNAs and 41 novel miRNAs were obtained from 'Guanhuahong' pitaya (H. monacanthus), among which 95 conserved miR-NAs belonged to 53 miRNA families. 26.79 Gb raw RNA-Seq data were generated and de novo assembled into 68, 505 transcripts, in which 39,737 were annotated. miRNAs and their target genes involved in betalain accumulation were compared at different developmental stages of pitaya fruit. Seven target genes were verified by 5′RACE and a tobacco transient expression system. Those Hmo-miRNAs negatively regulated expression of their target mRNAs through guiding corresponding target mRNA cleavage or inhibiting the translation. Hmo-miR157b-HmSPL6-like, Hmo-miR160a-Hpcyt P450-like3, Hmo-miR6020-HmCYP 71A8-like, Hmo-novel-2-HmCYP83B1-like, Hmo-novel-15-HmTPST-like, Hmo-miR828a-HmTT2-like, Hmo-miR858-HmMYB12-like, Hmo-miR858-HmMYBC1-like and Hmo-miR858-HmMYB2-like are possibly involved in betalain biosynthesis in pitaya. The present study provides new information that miRNAs are actively involved in betalain accumulation of pitaya fruit by regulating the upstream TFs, which may contribute to a further understanding of miRNAs in betalain biosynthesis of pitaya.

Plant materials
Two pitaya cultivars, i.e., 'Guanhuahong' (red peel with red pulp, H. monacanthus) and 'Guanhuabai' (red peel with white pulp, H. undatus) and Nicotiana benthamiana were used as plant materials. 'Guanhuahong' and 'Guanhuabai' pitayas, authenticated by Professor Guibing Hu and Yonghua Qin (College of Horticulture, South China Agricultural University), were selected from 860 seedlings of 'Hongshuijing' (H. monacanthus) [63][64][65]. Nicotiana benthamiana was grown in a greenhouse with a condition of 16 h/8 h day/night at 25°C and was used for interactions between miRNAs and their target gene assays in vivo. The South China Agricultural University provided all plant materials used in this study, and no specific permissions were required for the collection of those samples for research purposes following institutional, national and international guidelines. Fruits of 'Guanhuahong' and 'Guanhuabai' pitayas from the same orchard of Dalingshan Forest Park were separated into peels and pulps on the 13 th , 16 th , 19 th , 23 rd , 25 th , 27 th and 29 th DAF ( Figure S7) for expression analyses of crucial miRNAs and their targets. Pulps from 'Guanhuahong' pitaya on the 19 th DAF (white pulp stage, Hp19d, Hp19d_1, and Hp19d_2), 25 th DAF (pulp coloration stages, Hp25d, Hp25d_1, and Hp25d_2) and 29 th DAF (mature stage, Hp29d, Hp29d_1, and Hp29d_2) were used for RNA-Seq and sRNAome. All samples were frozen immediately in liquid nitrogen and stored at − 80°C until use.  [20,66]. Six small RNA libraries (Hp19d_1, Hp19d_2, Hp25d_1, Hp25d_2, Hp29d_1 and Hp29d_2) from two biological replicates and three RNA-Seq libraries (Hp19d, Hp25d and Hp29d) were constructed (https://dataview.ncbi.nlm.nih.gov/object/PRJNA588519? reviewer=ko91rr55muepqo1d25plnp0kv4). All HTS was performed by LC-BIO (Hangzhou, China).
Analyses of transcriptome raw data was identical to that of sRNAome raw data. Clean transcriptome data were assembled into non-redundant unigenes using Trinity (http://trinityrnaseq.github.io/). Unigenes were tentatively identified based on the best hits against known sequences in the database.

Prediction and functional annotation of target gene
Target genes of differentially-expressed miRNAs were predicted by Target Finder software. GO, KOG and KEGG were used to analyze the functions of target genes. E-value ≤10 − 5 was considered as significant enrichment.

Sequence alignment analyses
Multiple sequences were aligned using DNAMAN software (version 8). The miRNA target sites were identified using sequence alignments and manual analyses.

Analyses of miRNAs and target genes by RT-qPCR
Stem-loop RT-qPCR was used to confirm expression of miRNAs since it is a highly sensitive method for detection of miRNAs [68]. cDNAs were produced from 1.0 μg of total RNA samples using the MMLV-reverse transcriptase (Invitrogen) with miRNA specific stem-loop and oligo(dT) primers, respectively. The specific primers and PCR reactions were performed according to our previous method [20]. RT-qPCR was conducted in ABI 7500 real-time PCR System (Applied Biosystems, CA, USA) using the SYBR qPCR Mix (Vazyme). Twenty microliters reaction mixture contained 2.0 μL of diluted cDNAs (~15 ng/μL), 10.0 μL 2 × SYBR qPCR Mix (Vazyme), 0.5 μL of each primer (10 μM) and 7.0 μL ddH 2 O. All experiments were performed in triplicate. U6 and actin gene were used as reference genes [69,70]. The sequences of miRNAs and target genes were shown in Tables S2 and S10, respectively. All primers used for RT-qPCR analyses were listed in Tables S4, S5 and S6. The expression levels of miRNAs and target genes were calculated by 2 -△△C T method [71].

5′RACE analyses
To verify the miRNA-mediated cleavage events, RNA ligase-mediated 5′RACE (RLM-RACE) was performed using the SMARTer® RACE 5′/3′Kit User Manual (012615) (TaKaRa, Dalian, China) according to the manufacturer's manual. 1.0 μg total RNA from the pulps of the 29 th DAF 'Guanhuahong' and 'Guanhuabai' pitayas was ligated with 5′RNA adapters, respectively. The ligated mRNA was reversely transcribed by oligo (dT) primer. 5′ end products were obtained using 5′ adaptor primers and 3′ gene-specific primers. PCR products were inserted into pMD18-T vector (TaKaRa). Specific primers used for nested PCR were shown in Table S7.

Transient expression analysis
A transient expression system was used to confirm the interaction between miRNAs and their target genes in vivo. Construction of expression vectors were constructed following the procedure of Liu et al. [69]. Over-expression vectors of miRNAs related to the betalain biosynthesis and a control miRNA were constructed respectively. Co-expression of miRNAs and their targets in N. benthamiana leaves using Agrobacterium tumefaciens GV3101 infiltration. Transient expression in N. benthamiana was performed as described by Sparkes et al. (2006) [72]. Three days after infiltration, leaves were observed with a laser scanning confocal microscope (Zeiss LSM800) according to the following parameters: laser wavelength of 488 nm (laser intensity of 0.2%), master gain values ranging from 550 to 700 V, ESID gain value of 3 and digital gain value of 1.0. The transient expression assays were repeated at least three times.