Integrated metabolome and transcriptome analysis of the anthocyanin biosynthetic pathway in relation to color mutation in miniature roses

Background Roses are famous ornamental plants worldwide. Floral coloration is one of the most prominent traits in roses and is mainly regulated through the anthocyanin biosynthetic pathway. In this study, we investigated the key genes and metabolites of the anthocyanin biosynthetic pathway involved in color mutation in miniature roses. A comparative metabolome and transcriptome analysis was carried out on the Neptune King rose and its color mutant, Queen rose, at the blooming stage. Neptune King rose has light pink colored petals while Queen rose has deep pink colored petals. Result A total of 190 flavonoid-related metabolites and 38,551 unique genes were identified. The contents of 45 flavonoid-related metabolites, and the expression of 15 genes participating in the flavonoid pathway, varied significantly between the two cultivars. Seven anthocyanins (cyanidin 3-O-glucosyl-malonylglucoside, cyanidin O-syringic acid, cyanidin 3-O-rutinoside, cyanidin 3-O-galactoside, cyanidin 3-O-glucoside, peonidin 3-O-glucoside chloride, and pelargonidin 3-O-glucoside) were found to be the major metabolites, with higher abundance in the Queen rose. Thirteen anthocyanin biosynthetic related genes showed an upregulation trend in the mutant flower, which may favor the higher levels of anthocyanins in the mutant. Besides, eight TRANSPARENT TESTA 12 genes were found upregulated in Queen rose, probably contributing to a high vacuolar sequestration of anthocyanins. Thirty transcription factors, including two MYB and one bHLH, were differentially expressed between the two cultivars. Conclusions This study provides important insights into major genes and metabolites of the anthocyanin biosynthetic pathway modulating flower coloration in miniature rose. The results will be conducive for manipulating the anthocyanin pathways in order to engineer novel miniature rose cultivars with specific colors. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03063-w.


Background
Throughout human history, roses have held a high iconic cultural significance. Roses are used as garden ornamental plants and are employed in various industries such as food and cosmetic. There is a wide variation in flower colors of rose cultivars. In ornamental plants, the floral color is the main attribute and is closely related to the distribution of pigment types [1]. Natural pigments such as flavonoids, carotenoids, betalains, and alkaloids are well known in flower color formation [2][3][4].
Particularly, flavonoids have attracted more attention as important secondary metabolites in plants [5]. Flavonoids are known to play significant roles in various aspects of plant physiological processes, including protection from ultraviolet radiation, pigmentation, and plant defense [6]. Flavonoids are the products of the phenylpropanoid biosynthesis pathway which has been well characterized in Arabidopsis and petunia [7,8]. Several genes such as phenylalanine lyase, cinnamic acid hydroxylase, 4-coumarate-CoA ligase, chalcone synthase, chalcone isomerase, flavonoid 3-hydroxylase, flavonol synthase, dihydroflavonol 4-reductase, and anthocyanidin synthase, participate at different steps of the pathways. Anthocyanins are the end-products of the flavonoid pathway, but they are unstable in the cytoplasm and require further glycosylation. Genes of the UDP-glucosyltransferase and glutathione S-transferase families enter into action [9,10]. Conserved structural genes have been identified in the flavonoid biosynthetic pathway, but their regulatory mechanisms vary across plant species [5,7,8,[11][12][13][14].
As the key aesthetic characteristic of roses, flower color has been early investigated by breeders. Up to date, only few flavonoids have been reported in roses. Rose petals are known to contain anthocyanins, such as pelargonidin, cyanidin, and peonidin [15][16][17]. The presence of 3,5-diglycosyl anthocyanidins in association with 3-glycosylated flavonols confer the pink and red colors of roses [18,19]. Three genes (RhGT1, RhGT2, and RhGT3) encoding flavonoid 3-glycosyltransferases were found to be involved in the formation of anthocyanin glycosides [20]. Our understanding of flower coloration mechanisms in roses is far to be complete since several structural genes and as well as transcription factors (TF) regulating the anthocyanin pathways are yet to be reported. Moreover, the large variation existing in modern roses could indicate varying molecular mechanisms of color formation among different roses [21,22]. Highquality genomes of various Rosa species are now available and will facilitate fundamental research on these plants [23][24][25]. Although some genes involved in flower coloration in roses have been cloned, the lack of mutants for comparative studies makes it impossible to clearly understand the regulation of anthocyanin pathways [26][27][28][29][30][31][32][33][34][35][36][37][38].
Understanding the genetic and molecular mechanisms that modify coloration will not only answer a fascinating query involving fundamental rose biology but will also allow the manipulation of rose quality. By integrating metabolome and transcriptome analyses, researchers have identified key metabolites and genes (structural or TF) regulating color formation in plants. In this study, we used a miniature rose cultivar (Neptune King) and its natural color mutant (Queen) and investigated the molecular changes affecting anthocyanin composition in their petals.

Metabolome profiling of petal samples
Fresh petals were collected from the miniature rose cultivar Neptune King (H) having light pink color and its natural color mutant Queen (S) displaying deep pink  (Fig. 1a, b). To compare the flavonoid content, petal samples were analyzed using UPLC-MS/MS. The flavonoid profiles of H and S flowers showed marked differences (Figures S1 and S2, Supplementary Materials). Using a local metabolite database, a total of 190 flavonoid-related metabolites were identified, including 18 anthocyanins, 61 flavonols, 8 dihydroflavones, 3 dihydroflavonols, 12 flavonoid carbonosides, 7 isoflavones, 67 other flavonoids, 11 proanthocyanidins, 1 chalcone, and 2 tannins (Table S1, Supplementary Materials). We used the metabolite quantification data to construct a hierarchical clustering heatmap. The results showed a close relationship between biological replicates, a sign of high quality metabolome quantification (Fig. 2a, b). Moreover, it could be observed a clear separation between H and S petal samples, indicating distinct flavonoid profiles in H and S samples.
Differentially accumulated metabolites (DAMs) between H and S petal samples were determined based on the variable importance in projection ≥ 1 and a fold change ≥ 2 or ≤ 0.5. There were 45 metabolites whose accumulation was significantly different between the compared samples, and these differentially accumulated flavonoids (DAFs) were mainly anthocyanins, other flavonoids, procyanidins, and flavonols (Table 1). In particular, seven anthocyanins, namely, cyanidin 3-O-glucosyl-malonylglucoside, cyanidin O-syringic acid, cyanidin 3-O-rutinoside, cyanidin 3-O-galactoside, cyanidin 3-O-glucoside, peonidin 3-O-glucoside chloride, and pelargonidin 3-O-glucoside displayed significantly higher content in the S cultivar than in the H cultivar. The quantity of cyanidin O-syringic acid was high in the S cultivar but undetectable in the H cultivar. Additionally, the content of cyanidin 3-O-glucosyl-malonylglucoside was down-accumulated in the H cultivar compared to the S cultivar. The results indicate that these anthocyanins may play a crucial role for deep pink color formation in S flowers. Thus, the high accumulation of peonidin, pelargonidin, and cyanidin derivatives in the S cultivar may have led to its color-related changes.

Transcriptome sequencing and analysis
Global gene expression was further profiled in the petal samples of the two cultivars. From the six cDNA libraries (three biological replicates per sample), we obtained a total of 47.52 Gb of clean data, with more than 91% of bases scoring Q30 (  (Fig. 3b). This implies that the DAMs in   the two phenotypes were regulated by differentially expressed genes (DEGs). A total of 4,298 DEGs, including 1,851 upregulated and 2,447 downregulated genes ( Figure S2, Supplementary Materials) were identified between the samples.

Expression of genes involved in the flavonoid biosynthesis pathway
The expression patterns of genes involved in the flavonoid biosynthesis pathways were analyzed. A total of 150 genes were assigned to this pathway (Table S5, Supplementary Materials). Comparison of the H and S data showed 15 flavonoid related genes differentially expressed. All genes involved in the specific anthocyanin biosynthetic pathways, except for ANS, were upregulated in the S cultivar, which may be determinant for colorrelated changes ( Table 2).

Expression patterns of transcription factors
MYB and bHLH are particular transcription factors (TF) known to control the expression levels of structural genes involved in the biosynthesis of anthocyanins [39]. Hence, we searched or the DEGs encoding TF between the two cultivars. Only 30 DEGs encoding TFs were identified. Among these TFs, we identified one bHLH (RchiOBHmChr1g0361301) and   (Table 3). RchiOBHmChr1g0361301 (bHLH) and RchiOBHmChr1g0361191 (MYB) were strongly up-regulated in the S cultivar, suggesting a positive regulation of anthocyanin biosynthetic structural genes. In contrast, RchiOBHmChr1g0373541 (MYB) was not expressed at all in the petals of the S cultivar but was expressed in petals of the H cultivar. Comparison of the sequences of this transcript among the two cultivars did not show any nucleotide change, therefore it is probable that a mutation has occurred in the regulatory regions.

Expression patterns of MATE and ABCC gene family members
It has been reported that members of the multidrug and toxic compound extrusion (MATE) and the ATP-binding cassette transporter sub-family C (ABCC) families play preponderant roles for the sequestration of the anthocyanins from the cytoplasm to vacuoles [40,41]. We therefore expanded our study on the expression patterns of ABCC and MATE genes in miniature rose petals.

qRT-PCR validation of the expression patterns of anthocyanin related genes
To further validate the RNA-seq results, nine DEGs related to anthocyanin biosynthesis were selected, and their expression levels were analyzed in H and S petals at the blooming stage using qRT-PCR. The primer sequences of the genes are shown in Table S7, Supplementary Materials. The results confirmed that anthocyanin biosynthetic and regulatory genes, including DFR (RchiOBHmChr5g0037071, RchiOBHm-Chr6g0281711), CHI (RchiOBHmChr1g0372181), BZ1 (RchiOBHmChr4g0393121, RchiOBHmChr6g0302721), UGT79B1 (RchiOBHmChr5g0009571), MYB (RchiOBH-mChr6g0252211, RchiOBHmChr7g0241861), and bHLH (RchiOBHmChr1g0361301) were upregulated in petals of the S cultivar ( Figure S3, Supplementary Materials). We noticed a complete concordance between qRT-PCR and RNA-seq results, demonstrating that the RNA-seq data and DEG analysis in this study were reliable.

Association analysis of the metabolome and transcriptome
To demonstrate the synthetic and regulatory characteristics of the DAMs and DEGs, subnetworks were constructed to show transcript-metabolite correlations.
Only pairs with a correlation coefficient > 0.8 were considered in the analysis.
Four transcript-metabolite correlation networks were constructed. Among these, Ko00941 consisted of 15 nodes and 31 edges, with 26 pairs showing a positive correlation and 5 pairs showing a negative correlation (Fig. 4a). Ko00942 consisted of five nodes and six edges, with four pairs showing a positive correlation and two pairs showing a negative correlation (Fig. 4b, Supplementary Materials). Ko00943 consisted of nine nodes and eight edges, with four pairs showing a positive correlation and four pairs showing a negative correlation (Fig. 4c). Ko00944 consisted of five nodes and six edges, with three pairs showing a positive correlation and three pairs showing a negative correlation (Fig. 4d). Detailed information on the gene-metabolite pairs involved in flavonoid biosynthesis is listed in Table S8-1-S8-4, Supplementary Materials.
These results indicate that several classic flower pigmentation-related genes were highly correlated with their corresponding metabolites, which are involved in flavonoid biosynthesis, thereby reconfirming the importance of flavonoid biosynthesis in flower color.
Based on our results, we reconstructed the anthocyanin biosynthetic pathway in miniature roses. The identified anthocyanins and their relevant intermediates were rearranged to their corresponding positions (Fig. 5).
The key genes involved in the anthocyanin biosynthesis were significantly upregulated in petals of the S cultivar, which was consistent with the high accumulation of anthocyanins in S. Taken together, we deduce that the coordinated upregulation of CHS, DFR, FLS, BZ1, GT1, UGT78D2, and UGT79B1 together with the high expression levels of MATE genes play a crucial role in stimulating anthocyanin synthesis and accumulation, including cyanidin O-syringic acid, cyanidin 3-O-rutinoside, cyanidin 3-O-galactoside, cyanidin 3-O-glucoside, peonidin 3-O-glucoside chloride, and pelargonidin 3-O-glucoside, conferring the deep pink color in the S cultivar.

Discussion
The aesthetic feature of flowers is of central importance in the ornamental quality and commercial value of orchids. Flower coloration has been the focus of biological studies [42][43][44], which have shown that there are species-specific peculiarities in the regulation of pigmentation. To contribute to the existing knowledge, species-specific studies of flower color formation and Fig. 5 Flavonoid biosynthesis pathway and the differential metabolites of petals in miniature rose cultivar Neptune King (H) and natural mutant Queen (S). The red dots represent the up-regulated metabolites elucidation of specific regulatory mechanisms are necessary. In the present study, we investigated the underlying mechanism of color formation in two distinct flower phenotypes, the miniature rose cultivar Neptune King and its mutant, Queen rose.

Effects of the of anthocyanins distribution on flower coloration in the two rose cultivars
A total of 190 flavonoid-related metabolites were detected and categorized into eight subgroups, including 18 anthocyanins and 11 procyanidins. This provides the most complete view of the metabolites involved in flower coloration in miniature roses. Previous studies [15-17, 45, 46] on the rose coloration process have been limited to a few flavonoid metabolites. The composition of anthocyanidins and flavonols determines the diverse phenotypes of flowers [47,48]. The analysis of the DAMs showed that the mutant Queen was characterized by a significant high accumulation of various flavonoids. In the present study, we detected 45 differential flavonoid-related metabolites, including 7 anthocyanins. The findings indicate that cyanidin O-syringic acid, cyanidin 3-O-rutinoside, cyanidin 3-O-galactoside, cyanidin 3-O-glucoside, peonidin 3-O-glucoside chloride, and pelargonidin 3-O-glucoside accumulated to a greater extent in mutant flowers. Therefore, cyanidin derivatives can be considered the key anthocyanins responsible for the color change. It was reported that cyanidin and peonidin are responsible for orange to red colors [3,[49][50][51][52][53][54].
Anthocyanin contents have been regarded as the major factors endowing transgenic tobacco with a pink or light red flower color [55][56][57][58][59]. Similarly, cyanidin and peonidin have been shown to be responsible for the red color in various plants, such as tree peony, strawberry, apple, Chrysanthemum, Rhododendron and Prunus mume [2,3,[60][61][62][63]. In addition, a combination of several flavonoid metabolites results in a co-pigmentation effect, which leads to deeper colors [64][65][66]. In the present work, the levels of one luteolin derivative, three kaempferol, three quercetin derivatives, one apigenin, were simultaneously increased in S samples, indicating a co-pigmentation effect responsible for the deep pink color observed in the mutant Queen rose.

Major anthocyanin related genes affecting flower coloration in the two rose cultivars
Recently, high-throughput RNA-seq has been extensively used in studies on roses [67][68][69][70]. The pathways controlling flower coloration are properly characterized in plants, including roses [71,72]. However, because of the extensive genetic variations observed in roses, varying mechanisms for color formation are expected.
In this study, the transcriptomes of Neptune King and its color mutant were analyzed, and differentially expressed flavonoid biosynthetic genes affecting the color formation were identified. The results of RNAseq showed a significant increase in the expression of genes involved in the flavonoid biosynthetic pathway in S flowers, which strongly supported the findings of our metabolome analysis. Our findings indicated the upregulated expression of genes such as CHS, DFR, BZ1, UGT79B1, GT1, UGT78D2, FLS, and IF7MAT in S compared to H, implying that the high expression of these genes was the immediate cause for the higher accumulation of anthocyanins in mutant flowers. In previous reports, the expression of early or late biosynthetic genes, usually correlates positively with the anthocyanin content [73][74][75][76][77][78][79][80].
CHS is the primary rate-limiting enzyme that catalyzes the transformation of p-coumaroyl-CoA and malonyl-CoA into tetrahydroxy-chalcone in flavonoid biosynthesis. Weak expression of CHS restrains anthocyanin biosynthesis and leads to color changes in plant tissues. Variation of CHS expression levels confers color polymorphism in crabapple and arctic mustard flowers [81,82]. In the present work, we deduce that the high expression levels of CHS in S flowers would provide sufficient amounts of precursor compounds for the higher anthocyanin content.
Anthocyanidins are extremely unstable and degenerate easily; therefore, glycosylation is crucial to stabilize and transport these compounds to vacuoles, where they can function as pigments. This glycosylation role is played by UFGT genes [7,90]. In the present study, RhGT1 expression was upregulated in S flowers, which may account for the deep pink color. In accordance to our study, Griesser et al. [91] demonstrated that a downregulation of GT1 resulted in of a significant reduction of anthocyanin pigments in ripe strawberry fruits.
MYB and bHLH are the major transcription factors regulating the expression of structural genes involved in the anthocyanin biosynthesis [39]. In this study, we identified three key genes that could be further investigated in order to facilitate the modulation of color in miniature roses.

Downregulation of MATE genes could reduce anthocyanin sequestration in the vacuoles of H flowers
Anthocyanins, particularly the hydrophobic aglycone forms, are toxic molecules for the cell machinery due to their high chemical reactivity. Hence, after synthesis, a mechanism of anthocyanin removal is immediately triggered, leading to a sequestration in the central vacuole. High vacuolar accumulation of anthocyanins confers the observed color phenotypes in plant tissues. MATE and ABCC members have been reported as key genes involved in anthocyanins sequestrations in the vacuoles [40,41]. Debeaujon et al. [40] were the first to discover the gene TRANSPARENT TESTA12 (TT12), a MATE member, necessary for flavonoid sequestration in vacuoles of the seed coat endothelium in Arabidopsis. Later on, TT12 has been identified, isolated and functionally validated in various plant species such as Brassicas, grapevine and blueberry [92][93][94]. In this study, we identified eight TT12 upregulated in petals of the S cultivar, indicating that these genes may contribute to the high accumulation of anthocyanins observed in S.

Conclusions
In this study, higher contents of diverse flavonoids and alterations in gene expression patterns were observed in the rose mutant as compared to its wild type. The mutation induced a significant upregulation of the key genes participating in the anthocyanin synthesis pathway, glycosylation and sequestration, resulting in significantly increased anthocyanin accumulation and darker flower coloration in the S cultivar. Functional characterization of the candidate structural genes as well as the three major transcription factors will continue.

Plant materials and sampling
The miniature rose cultivar Neptune King (H) and its natural mutant Queen (S) were cultivated at the Liaoning Academy of Agricultural Sciences in Shenyang City (41°48′11.75″N, 123°25′31.18″E), China. The formal identification of the plant materials was undertaken by Professor Hongmei Sun. The plants are kept at the Key Laboratory of Protected Horticulture of Education Ministry and Liaoning Province, College of Horticulture, Shenyang Agricultural University, Shenyang, China. The petal samples were collected for the metabolome study, RNA-Seq and qRT-PCR analysis. Three biological replicates were gathered per sample, with 10 blooms flowers randomly collected from 10 plants in the same batch. The petals were frozen in liquid nitrogen, and saved at − 80 °C until further use. The H and S samples used for the metabolome analysis were designated DH and DS, respectively, and those used for the RNA-seq analysis were designated ZH and ZS, respectively.

Sample preparation and extraction
Flavonoid metabolites were extracted following a previous protocol [65], with some modifications. Petal samples were vacuum freeze-dried and ground (1.5 min at 30 Hz) to powder using a grinder. The powder (0.1 g) was weighed and dissolved in a 70% aqueous solution of methanol. To improve extraction efficiency the dissolved sample was swirled three times then stored at 4 °C overnight. The supernatant was filtered with a micropore membrane (0.22 μm pore size) after centrifugation (10,000 × g, 10 min), and stored in a sample injection bottle for UPLC-MS/MS analysis.

Flavonoid identification, quantification, and data analysis
All samples were analyzed by a UPLC-MS/MS system (Shim-pack UFLC SHIMADZU CBM30A and Applied Biosystems 6500 QTRAP). A scheduled multiple reaction monitoring method was used to metabolite quantification.
Based on the self-built MWDB database (Metware Biotechnology Co., Ltd., Wuhan, China) and a public database of metabolite information, the fundamental and secondary MS records had been qualitatively analyzed using Analyst 1.6.3. The identified metabolites with significant difference in content were set with 0.5 ≥ fold change ≥ 0 or a fold change ≥ 2, p-value < 0.05, and VIP ≥ 1 were considered DAFs. Finally, the KEGG pathway database (http:// www. kegg. jp/ kegg/ pathw ay. html) and MWDB were centered on metabolic reactions and concatenated possible metabolic pathways.

RNA extraction, quantification, and sequencing
RNA extraction from petal samples, quantification and sequencing were conducted according to the procedures detailed by Dossa et al. [95,96]. Six libraries representing the collected H and S petal samples (three biological replicates respectively) were prepared following standard procedures of Illumina HiSeq 4000 platform.

Transcriptome data analysis
The raw data were cleaned by removing adaptors and lowquality reads using Fastqc with default parameters. Then, the clean reads were mapped to the R. chinensis genome (https:// lipm-brows-ers. toulo use. inra. fr/ pub/ RchiO BHm-V2/) using HISAT tool. Fragments per kilobase of exon model per million reads mapped values were used for gene/ transcript measurement. Differential expressed gene (DEG) analysis was conducted using the DESeq2 tool. Genes with |log 2 Fold Change|≥ 1 and a p-value < 0.05, were described as DEGs [96]. GO and KEGG enrichment analyses of the DEGs were conducted in the R package 'clusterProfiler' .

qRT-PCR analysis
Nine anthocyanin-related genes were targeted for qRT-PCR. The qRT-PCR was performed based on the same RNA samples used for sequencing [95]. The gene Actin was used as the internal control for transcript expression normalization. The primer pairs were designed using the Primer5 tool. The cDNAs were synthesized using Super-Script ™ III First-Strand Synthesis SuperMix for qRT-PCR (Invitrogen). The qRT-PCR was conducted on an ABI 7500 Fast real-time detection system (Applied Biosystems) The reaction mix contained: 5 μl of FastStart Universal SYBR Green Master (Roche Life Sciences), 2 μl of 2 μM primer mix, 2 μl of a diluted 1:10 cDNA and ddH 2 O was added to a final volume of 10 μl. The cycling conditions were set as follow: 95 °C for 10 min, and 40 cycles of 95 °C for 15 s, and 60 °C for 1 min. The experiment was conducted with three biological replicates and three technical replicates. Data was analyzed according to the 2 −∆∆Ct method [97].