The plant materials consisted of flowers from the purple and white varieties of M. malabathricum that had been grown under the same conditions in the peony base of the Guangzhou Institute of Landscape Gardening (Guangdong, China) in July 2016. Petal samples were collected at three flower color developmental stages for flavonoid content analysis and RNA extraction (Figure 1). Flowers were sampled from separate plants for each stage. The plants were at different stages, and the two varieties were randomized locations.

Total RNA was isolated from each sample using a TRIzol® Reagent RNA Isolation Kit (Invitrogen, Grand Island, NY) following the manufacturer’s protocol. Three biological replicates of RNA were extracted for each sample. RNA degradation was detected on 1.2% agarose gels, and RNA purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). Approximately 5 μg of total RNA (pooled in equal amounts from three biological replicates) for each sample was used for library construction. mRNA was enriched by oligo(dT) beads, and sequencing libraries were created using an Illumina TruSeq™ RNA Sample Prep Kit (Illumina, SD, USA) following the standard high-throughput protocol. Sequencing was conducted on an Illumina HiSeq™ 4000 platform at Gene Denovo Biotechnology Co. (Guangzhou, China), and paired-end reads were generated.

Raw reads in the FASTQ format were first processed through in-house Perl scripts. Clean reads were obtained by removing reads containing adapters, reads containing poly-N, and low-quality reads from raw data. At the same time, the Q20, Q30, GC content, and sequence duplication level of the clean data were calculated. The remaining clean reads were assembled without a reference by Trinity for the two M. malabathricum varieties with a min_kmer_cov set to 2 by default, and all other parameters were set to the default (Grabherr et al., 2011). The obtained sequences were defined as unigenes. All of the unigenes were subjected to the Nr, SwissProt, COG, KEGG, and Pfam public databases to annotate gene functions. To identify the homologous genes of two M. malabathricum varieties by similarity searching with an E < 10⁻⁵ using BLASTn, the GATK3 software (McKenna et al., 2010) was used to perform SNP calling. Raw vcf files were filtered with the GATK standard filter method with default parameters.

High-quality clean reads were compared to de novo assemblies using Bowtie2, and then the RPKM (reads per kb per million read) values were calculated and normalized (Mortazavi et al., 2008). To identify DEGs across samples without biological replicates, for each sequenced library, the read counts were adjusted by edgeR package through one scaling normalized factor (Robinson et al., 2010). Differential expression analysis of any two samples was performed using the edgeR package. We identified genes with a log2 fold change ≥2 and a corrected p <0.01 in comparison as significant DE genes. We used the KOBAS software with a hypergeometric test and the Benjamini-Hochberg FDR correction to test the statistical enrichment of DE genes in KEGG pathways (Xie et al., 2011).

Flavonoid content was determined according to commonly used spectrophotometric methods based on the formation of a complex flavonoid, aluminum, which has an absorption maximum at 510 nm (Pękal and Pyrzynska, 2014), and using quercetin as the internal standard. Briefly, the flower petal materials were dried to a specific weight, chopped very finely, and sieved through a plastic sieve (<0.5 mm). The dried material (1.0 g) was crushed by ultrasonication and extracted for 1 h with 2.5 ml of 60% (v/v) aqueous ethanol at 60°C. Then, the precipitate solution was separated by centrifugation, and pellets were discarded. One milliliter of the diluted supernatant sample was separately mixed with 1 ml of 2% aluminum chloride methanol solution (Djeridane et al., 2006). After incubation at room temperature for 15 min, the absorbance of the reaction mixture was measured at 510 nm with a Milton Roy 601 UV-vis spectrophotometer, and the flavonoid content was expressed in milligrams per gram of dry weight. All analyses were carried out with at least four replicates, and ANOVA was performed on the means of the flavonoid content. SPSS Statistics v. R23.0 (IBM) was used for data analysis with ANOVA (PROC ANOVA). Treatment means were separated at the 0.05 significance level by Tukey’s test.

To validate the transcriptome, RNA samples from three biological replicates of each stage were adjusted to 1 μg/μl with nuclease-free water. Two micrograms of total RNA was reverse transcribed in a 20-μl reaction volume using the Prime Script™ RT reagent Kit with gDNA Eraser (Takara, Dalian, China). Nine genes were chosen, and quantitative PCR was performed with an ABI Step-One plus Real-Time PCR system (Applied Biosystems, Canada). Each reaction was carried out in triplicate, and the primers used for qPCR were designed based on the transcriptome data except for actin and are listed in Supplementary Table S1. Relative expression levels of target genes were determined by the 2^−△△Ct method (Livak and Schmittgen, 2001). The correlation between the expression profiles of the selected genes measured by qPCR and RNA-seq was determined using the R package.

To profile the flower transcriptome of two M. malabathricum varieties with contrasting purple and white petals, we sequenced six RNA-seq libraries from three successive developmental stages of flowers, including young bud (stage 1), buds about to open (stage 2), and buds fully opened (stage 3) from the purple and white varieties of M. malabathricum (Figure 1). The early petal color of the purple variety is unpigmented, slightly pigmented, and completely turned to purple during stages 1–3. However, the petal color of the white variety in the three stages is white.

After adapter removal and quality trimming, we obtained a total of 195,675,054 reads for the purple variety and 180,277,772 reads for the white variety, for a total of 29 and 26 Gb of clean nucleotides, respectively. The percentages of clean reads among the raw tags (Q20) ranged from 97.41 to 97.65% (Table 1). We conducted de novo assembly of the transcriptomes of the two M. malabathricum varieties by Trinity (Grabherr et al., 2011) and obtained 52,498 and 49,380 unigenes for the purple and white varieties, respectively. The average unigene length and other characteristics of the assembly in the two M. malabathricum varieties were similar (Table 1 and Supplementary Figure S1). The sequencing coverage (estimated as the mean number of reads per unigene) was assessed as 3,169 for the purple variety and 3,052 for the white variety. BUSCO analysis (Simão et al., 2015) revealed a completeness score of approximately 70% for each assembly (Supplementary Figure S2). These results demonstrate that the assembled unigenes could be a useful source of relatively complete coding sequences for M. malabathricum.

In total, 34,856 (67%) and 34,715 (70%) unigenes were successfully annotated with a significance threshold (E ≤ 1E−5) in at least one of the Nr, SwissProt, COG, and KEGG public databases for the purple and white varieties, respectively, and the results are summarized in Table 2. Interestingly, for M. malabathricum, the top hit species was Eucalyptus grandis, or flooded gum, followed by Theobroma cacao and then Cephalotus follicularis (Supplementary Figure S3A), suggesting that it is more closely related to flooded gum. Annotated unigenes were categorized into 25 functional groups according to the COG classification (Supplementary Figure S3B). The distribution of COG categories in both varieties is largely similar; only some categories showed slight differences between them. For example, the purple variety showed a considerably higher number of genes than the white variety, including general function prediction, signal transduction mechanisms, and transcription. In turn, two pathways were particularly abundant in the white variety, “energy production and conversion” and “translation, and ribosomal structure, and biogenesis.” These annotated transcriptomes can facilitate the discovery of lineage-specific biological functions of M. malabathricum genes. Additionally, the transcripts were annotated against all plant transcription factors (TFs), resulting in 1,597 and 1,603 annotated transcripts for the purple and white varieties, respectively. Comparative analysis of the TF gene family compositions revealed significant differences between the two varieties (Supplementary Figure S3C, Supplementary Table S2). This exceptionally large and diverse inventory of TFs may play an important role in floral ontogeny in shrubs.

As mentioned earlier, the raw data characteristics and assembly, as well as the annotation characteristics, were largely similar for the purple and white varieties. We further investigated these transcriptomes for homology based on the corresponding BLASTn hits that aligned to one another (E < 1E−5). In total, 41,026 unigenes covered more than 70% of their best BLAST hit and showed that there is a strong correlation between the assembly length of predicted unigenes that are orthologous in the purple variety and the white variety (Figures 2A,B). To make a comparison, SNPs in the assembled M. malabathricum transcriptome were also examined. As a result, the distribution of each SNP type was similar in the two varieties (Figures 2C,D). Our results also showed similarity in the frequency distribution of RNA editing in three stages of flowering for either variety, indicating that this is a general pattern of RNA editing across the floral developmental stages. However, there is a clear distinction in the frequency distribution for RNA editing between the two varieties, and the editing efficiency was significantly higher in the purple variety compared to the white variety (Figures 2E,F), which may be directly linked to the expressed genes of the contrasting genotypes. Taken together, comparative analysis of two transcriptomes demonstrated their overall similarity but also revealed a presumably lineage-specific transcriptomic signature.

To estimate the expression levels, high-quality reads from each sample were mapped to the corresponding M. malabathricum variety transcriptome assembly. The overall transcriptional activity of all the genes was normalized, and abundance was estimated based on RPKM values. In general, the distributions of the expression levels of most of the genes were similar for the two M. malabathricum varieties in the corresponding stages. Higher expression levels were observed in flower development stage 1 followed by stage 2, and the lowest expression was observed in stage 3 (Figure 3A). To observe the orthologous gene expression patterns, a heatmap that represents the hierarchical clustering of homologous transcripts between different developmental stages revealed that samples of the two M. malabathricum varieties at the same stage were grouped together and that stage 1 was strongly differentiated from stage 2 and stage 3 (Figure 3B). Similarly, despite the small relative distinction between varieties at the same stage, specific gene expression levels for individual plant varieties were clearly presented in this study. As expected, the correlation of the same stage between the purple variety and the white variety was high (r from 0.886 to 0.918) compared with the correlation of the different stages from the purple or white variety (Figure 3C). These data revealed that many shared transcriptional signatures across floral developmental stages were activated in both M. malabathricum varieties.

Pairwise comparisons among the preceding stages, with respect to each M. malabathricum variety, were performed to gain insight into the gene expression regulation of flower development. A false discovery rate (FDR) cutoff of 0.01 and log2 fold change of ±2 were set as the thresholds for putative differentially expressed genes (DEGs). The number of DEGs at the different developmental stages is shown in Supplementary Table S3. The enrichment analysis with KEGG terms (FDR < 0.05) in both DEG sets is largely similar, with only a small number of categories being unique for either variety. The four major metabolic pathway categories included (1) secondary metabolite synthesis, such as phenylpropanoid, flavonoid, cutin, wax, and diterpenoid biosynthetic pathways; (2) plant hormone signaling and production; (3) energy homeostasis and nutrient assimilation pathways, including the energy substrate biosynthesis and interconversions (such as pentose, glucuronate, starch and sucrose, fatty acids, and steroid biosynthesis), amino acid, nucleotide, and glutathione metabolism and photosynthesis; and (4) cellular proliferation, which appeared significantly overrepresented in flower development in the two varieties of M. malabathricum (Figure 4). It is also worth mentioning that most of the pathways were found to be significantly enriched in the last half of floral development, from stages 2 to 3 (stage 3 vs. stage 2). In contrast, in the initial pigmented period of stages 1 and 2 (stage 2 vs. stage 1), more KEGG enrichment categories were only marginally significant or not significant (FDR > 0.05), except for a limited number of pathways, such as phenylpropanoid biosynthesis, glutathione metabolism, and DNA replication.

Although many of the enriched KEGG terms were common, several of them were unique in different sets of genes. For example, phenylpropanoid biosynthesis was the most represented in a linage-specific expression manner in the purple variety compared to the white variety over flower developmental stages 1–3. Furthermore, three overrepresented pathways related to flavonoid biosynthesis, diterpenoid biosynthesis, and plant hormone signal transduction were also identified as the specific category for the initial flower pigmentation developmental process (P1 vs. W1). These results suggest that there is a potentially unique constitutive metabolism in initial flower development in each variety, possibly representing differences in lineage-specific molecular expression-derived morphological characteristics.

To assess the accuracy of the pathway enrichment from the DEG analysis obtained using RNA-seq data, we quantified flavonoid production in all three flower development stages in both M. malabathricum varieties to validate whether flavonoid biosynthesis was overexpressed (Figure 5). The results of the variance analysis showed that the accumulation of flavonoids drastically decreased with the flower developmental stages, particularly in the purple variety. The flavonoid content was significantly higher in white petals in stage 1 relative to that in purple petals (p < 0.01). Furthermore, we performed qPCR analysis for six selected genes associated with flavonoid synthesis, and the qPCR results were consistent with the RNA-seq data (r² = 0.83). The temporal expression patterns were highly similar between the M. malabathricum purple and white varieties, with the exception of the magnitude of differences; for example, the magnitude of variation of F3’H and ANS was markedly increased in the purple variety. The final biosynthetic steps for flavonoids that were catalyzed by FLS showed strong downregulation in the developing phase, which had an overall effect of decreasing the pool of FLS transcripts that, in turn, could limit the metabolic flux of flavonoids. The enzymes of ANS that catalyze the late steps of anthocyanin biosynthesis had significantly higher transcript induction for the purple variety development at both stage 2 and stage 3 relative to that observed for the white variety; low expression of the ANS gene resulted in little accumulation of anthocyanins in the white variety, which supported the hypothesis regarding its importance in the regulation of plant colors (Shimizu et al., 2011). These analyses suggest that there is a temporal regulation of flavonoid biosynthesis with flower development in both varieties, which agrees with the overrepresentation of the flavonoid biosynthesis pathway in the transcriptome analysis. Thus, flavonoid biosynthesis is important to the regulation of flower color development in M. malabathricum.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.