Comparative Metabolite and Gene Expression Analyses in Combination With Gene Characterization Revealed the Patterns of Flavonoid Accumulation During Cistus creticus subsp. creticus Fruit Development

Cistus creticus L. subsp. creticus (rockrose) is a shrub widespread in Greece and the Mediterranean basin and has been used in traditional medicine as herb tea for colds, for healing and digestive hitches, for the treatment of maladies, as perfumes, and for other purposes. Compounds from its flavonoid fraction have recently drawn attention due to antiviral action against influenza virus and HIV. Although several bioactive metabolites belonging to this group have been chemically characterized in the leaves, the genes involved in their biosynthesis in Cistus remain largely unknown. Flavonoid metabolism during C. creticus fruit development was studied by adopting comparative metabolomic and transcriptomic approaches. The present study highlights the fruit of C. creticus subsp. creticus as a rich source of flavonols, flavan-3-ols, and proanthocyanidins, all of which displayed a decreasing trend during fruit development. The majority of proanthocyanidins recorded in Cistus fruit are B-type procyanidins and prodelphinidins, while gallocatechin and catechin are the dominant flavan-3-ols. The expression patterns of biosynthetic genes and transcription factors were analyzed in flowers and throughout three fruit development stages. Flavonoid biosynthetic genes were developmentally regulated, showing a decrease in transcript levels during fruit maturation. A high degree of positive correlations between the content of targeted metabolites and the expression of biosynthetic genes indicated the transcriptional regulation of flavonoid biosynthesis during C. creticus fruit development. This is further supported by the high degree of significant positive correlations between the expression of biosynthetic genes and transcription factors. The results suggest that leucoanthocyanidin reductase predominates the biosynthetic pathway in the control of flavan-3-ol formation, which results in catechin and gallocatechin as two of the major building blocks for Cistus proanthocyanidins. Additionally, there is a decline in ethylene production rates during non-climacteric Cistus fruit maturation, which coincides with the downregulation of the majority of flavonoid- and ethylene-related biosynthetic genes and corresponding transcription factors as well as with the decline in flavonoid content. Finally, functional characterization of a Cistus flavonoid hydroxylase (F3′5′H) was performed for the first time.

Cistus creticus L. subsp. creticus (rockrose) is a shrub widespread in Greece and the Mediterranean basin and has been used in traditional medicine as herb tea for colds, for healing and digestive hitches, for the treatment of maladies, as perfumes, and for other purposes. Compounds from its flavonoid fraction have recently drawn attention due to antiviral action against influenza virus and HIV. Although several bioactive metabolites belonging to this group have been chemically characterized in the leaves, the genes involved in their biosynthesis in Cistus remain largely unknown. Flavonoid metabolism during C. creticus fruit development was studied by adopting comparative metabolomic and transcriptomic approaches. The present study highlights the fruit of C. creticus subsp. creticus as a rich source of flavonols, flavan-3-ols, and proanthocyanidins, all of which displayed a decreasing trend during fruit development. The majority of proanthocyanidins recorded in Cistus fruit are B-type procyanidins and prodelphinidins, while gallocatechin and catechin are the dominant flavan-3-ols. The expression patterns of biosynthetic genes and transcription factors were analyzed in flowers and throughout three fruit development stages. Flavonoid biosynthetic genes were developmentally regulated, showing a decrease in transcript levels during fruit maturation. A high degree of positive correlations between the content of targeted metabolites and the expression of biosynthetic genes indicated the transcriptional regulation of flavonoid biosynthesis during C. creticus fruit development. This is further supported by the high degree of significant positive correlations between the expression of biosynthetic genes and transcription factors. The results suggest that leucoanthocyanidin reductase
As the F3O and PA biosynthetic route is a branch of the flavonoid pathway, it has to be co-ordinately regulated with other flavonoid branches in different tissues and organs and in response to various developmental and environmental cues. Biosynthetic genes are turned on and off by specific transcription factors (TFs), which can simultaneously control several genes, or a single step within the pathway. Among plant TFs, the MYB factor proteins (V-myb myeloblastosis viral oncogene homolog) are most important within the flavonoid pathway (Martin et al., 2001). Following MYB factor binding to specific DNA regulatory elements in the promoter regions of target genes, transcriptional activation is initiated. Unlike many others, the PA-and anthocyanin-specific MYBs also need to partner with basic helix-loop-helix (bHLH) and WD-40 repeat proteins, forming the so-called MBW complex to promote transcription (Feller et al., 2011;Montefiori et al., 2015).
The aim of the present work was to analyze the polyphenolic composition of C. creticus subsp. creticus flowers and fruit, which, to our knowledge, have not been phytochemically characterized before. Species belonging to the genus Cistus primarily develop five-valve fruits (Demoly and Montserrat, 1993), with colors varying from green in juvenile to dark red in ripe fruit. Similarly, as in other fruit species, Cistus fruit development from flower to ripe stage undergoes metabolic alterations regulated by both developmental and hormonal factors. The accumulation of major groups of flavonoids in fruit might be the result of complex changes in the expression of structural/biosynthetic and regulatory genes involved in their metabolism. Therefore, in this work, the expression of structural genes involved in flavonoid biosynthesis and regulatory genes like TFs controlling this pathway was studied in parallel with the content of the major flavonoids during fruit development and senescence. Next, the interplay between ethylene and flavonoid accumulation in Cistus fruit was investigated. Lastly, a Cistus flavonoid hydroxylase (F3 5 H) was functionally characterized.

Plant Material
Cistus creticus L. subsp. creticus plants were grown in the experimental field of the School of Agriculture, Aristotle University of Thessaloniki, located in the area of Thermi, Thessaloniki, Greece (N 40.536247 • and E 22.993830 • ). Whole mature flowers and fruits of three different developmental stages were collected from three groups comprising five individual plants, and each group represented a biological replicate. Fruits of three developmental stages were distinguished based on the size and color of sepals and valves (Figure 2A). Small-sized halfexpanded fruits with green capsules and sepals were marked as stage 1 (S1). Stage 2 (S2) was composed of fully expanded fruits with the sepals starting to change color and detached from valves, while stage 3 (S3) comprised fully colored ripe fruits. Flowers were appointed as the control organ and marked as stage 0 (S0). Upon harvest, all samples were frozen in liquid nitrogen and kept at −80 • C until further use.

UHPLC-MS Orbitrap Identification of Phenolics in Cistus Fruit
For chemical analyses, the frozen plant material was powdered in mortars and pestles and subsequently lyophilized. Around 100 mg of each sample was extracted with 1 ml of 99.8% methanol (AppliChem, United States) by 1 min of vortexing and subsequent sonification in ultrasonic bath (RK100, Bandelin, Berlin, Germany) for 20 min. The samples were centrifuged for 10 min at 10,000 g, and the supernatants were filtered. The samples were kept at 4 • C until analyses. All extractions were performed in triplicate.
FIGURE 2 | Metabolic profiling of Cistus flower (S0) and three fruit stages (S1-S3). (A) HPLC/DAD elution profiles at λ = 260 nm; notable peaks are denoted in the legend on the side. The Y -axis scales are different for each of the chromatograms. (B) Biplot of principal component analysis performed on zero-centered and unit-scaled compound quantity data. The samples are colored by stage (S0-S3) as shown in the legend, and these groups are further encircled by a convex hull of the same color. Variable loadings are indicated by arrows and were scaled by multiplying with five prior graph constructions so they could be seen more readily. (C) Heat map of zero-centered and unit-scaled compound quantity data. The samples (in rows) are arranged according to hierarchical cluster analysis (Ward's method of cluster agglomeration) constructed using Euclidean distances (left tree), and the metabolites (in columns) are arranged according to hierarchical cluster analysis (cluster agglomeration using complete linkage) constructed using Spearman correlation distances (top tree).
Chromatographic separations of compounds in methanol extracts of Cistus flowers (S0) and fruit (S1-S3) were performed using an Accela 600 ultrahigh-performance liquid chromatography (UHPLC) system coupled to a linear ion trap-OrbiTrap hybrid mass spectrometer (LTQ OrbiTrap MS) (Thermo Fisher Scientific, Bremen, Germany). All chromatographic and MS settings [heated electrospray ionization (HESI) and the other MS parameters] were the same as in Banjanac et al. (2017). The injection volume was 5 µl.
The mass spectrometer (MS) was operated in negative ionization mode, and MS spectra were acquired by full-range acquisition covering 100-1,500 m/z. The resolution was set at 30,000 for full-scan (FS) analysis, which was employed to detect the monoisotopic masses of unknown compounds. Fragmentation pathways were proposed by multistage mass spectrometry (MS n ). The ions of interest were isolated in the ion trap with an isolation width of 5 ppm and activated with 35% collision energy (cE) levels. FS analysis was employed, while Xcalibur software (version 2.1) was used for instrument control and data analysis.

UHPLC/DAD/(±)HESI-MS/MS Quantification of Phenolics in Cistus Flowers and Fruits
UHPLC/diode array detector (DAD)/(±)HESI-MS/MS method was developed for the separation, identification, and quantification of targeted phenolic compounds in samples (methanol extracts) of Cistus flowers and fruits. Dionex Ultimate 3000 UHPLC system (Thermo Fisher Scientific, Bremen, Germany) was equipped with a DAD and connected to a triplequadrupole MS (TSQ Quantum Access Max, Thermo Fisher Scientific, Basel, Switzerland). Elution was performed at 40 • C on Syncronis C18 column (100 mm × 2.1 mm) with 1.7-µm particle size (Thermo Fisher Scientific, Bremen, Germany). The mobile phase consisted of (A) water + 0.1% formic acid and (B) acetonitrile, which were applied in the gradient elution previously described in Mišić et al. (2015), and with a flow rate of 0.4 ml min −1 . Acquisition of UV spectra was performed at λ = 260 and 320 nm, and the injection volume was set to 5 µl. The HESI source of MS was operated with a vaporizer temperature of 350 • C, while the ion source settings were as follows: spray voltage, 3,510 V; sheet gas (N2) pressure, 28 AU; ion sweep gas pressure, 0.0 AU, and auxiliary gas (N2) pressure, 4 AU; capillary temperature, 270 • C; and skimmer offset, 0 V. Selected reaction monitoring (SRM) experiment was conducted for the quantification of targeted phenolics, using argon as the collision gas and cE of 30 eV. The parameters of the adopted SRM UHPLC/DAD/(−)HESI-MS 2 method are presented in Supplementary Table 1, including the major MS 2 fragments and λ max values of targeted compounds. Instrument control, data acquisition, and processing analysis were performed using Xcalibur software (version 2.2). The phenolics were quantified based on the calibration curves of commercial standards: gallic acid, caffeic acid, quinic acid, rutin hydrate, myricetin, luteolin, quercetin, naringenin, apigenin, catechin, epicatechin, gallocatechin, epigallocatechin, epigallocatechin gallate, and gallocatechin gallate. Myricetin 3-O-rhamnoside was quantified relatively using the calibration curve of myricetin. The total amount of each targeted compound is expressed as µg 100 mg −1 dry weight (DW).

Ethylene Production and Respiration Rates
Air samples were taken after closing the fruit of all stages in 0.5-L glass jars for 1 h, according to Koukounaras et al. (2010). CO 2 concentration was measured by injecting the gas samples into a stream of N 2 carrier gas flowing through a CO 2 /O 2 analyzer (model Combo 280, David Bishop Instruments, United Kingdom), while ethylene concentration was measured by injecting the gas sample into a Varian 3300 gas chromatographer (Varian Instruments, Walnut Creek, CA, United States) equipped with a flame ionization detector. The ethylene production and respiration rates were expressed as µl C 2 H 4 kg −1 h −1 and mg CO 2 kg −1 h −1 , respectively.

Selection of Candidate Genes and qPCR Analysis
Putative candidate genes involved in the flavonoid pathway were identified by a comparative analysis of C. creticus fruit transcriptome BLAST search against the NCBI public database. Although the BLAST search resulted in several hits for some genes involved in flavonoid and ethylene pathway, only the full-length sequences retrieved from C. creticus fruit transcriptome were selected for this study. Based on the detected sequences, highly specific primer pairs for qPCR were designed using NCBI Primer-BLAST 1 (Supplementary Table 2). The sequences of characterized CcTTG1 and CcSPBPA/B from Cistus transcriptome cDNA library (Falara et al., 2008;Ioannidi et al., 2016) have GenBank accession numbers as follows: CcTTG1: KT892927, CcSPLA: KU145276, CcSPLB: KU041720, and CcF3 5 H: MT707661.
For the analysis of ethylene biosynthesis, total RNA isolation was performed from C. creticus flowers and the three developmental stages of fruit in biological triplicates using Sigma Spectrum Kit (Sigma Aldrich, Germany), with a slight modification of the manufacturer's instructions, namely, RNA was extracted from 50 mg of tissue; the volume of lysis buffer used was 1.5 ml per sample, whereas the volume of binding buffer was 3 ml per sample. To remove traces of genomic DNA, the RNA samples (1 µg) were treated with DNase I (Fermentas, Vilnius, Lithuania) in a final reaction volume of 10 µl. RNA quality was confirmed using NanoDrop 2000C Spectrophotometer (Thermo Scientific, United States), and their integrities were assessed by agarose gel electrophoresis.
First-strand cDNA was synthesized from 1 µg of total RNA using the SuperScript TM III Reverse Transcriptase (Thermo Fisher Scientific) following the manufacturer's instructions. The expression of genes was analyzed by real-time PCR using Light Cycler QuantStudio 3 (Thermo Fisher Scientific) and KAPA SYBR R Fast qPCR Master Mix (2X) Universal (KAPA Biosystems, United States). The general thermocycler conditions were 95 • C for 4 min and 40 cycles of 95 • C for 10 s, 64 • C for 20 s, and 72 • C for 15 s. The relative expression values were normalized against elongation factor-1α (EF062868.1) as the endogenous control and calculated by the 2 − Ct method (Livak and Schmittgen, 2001). The data represent the means ± SD of triplicates.

Phylogenetic Analysis
The phylogenetic analyses were performed as previously described in Aničić et al. (2020). Flavonoid-related protein sequences of MYB and nucleotide sequences of bHLH TFs were aligned using Clustal W within the software package MEGA 6.0. Gblocks v0.91b, with default parameters, was utilized to select highly conserved blocks of alignment positions. The maximum likelihood method within MEGA version 6.0 with default settings (NNI heuristic method, BioNJ initial tree, JTT model for MYB tree, and Tamura-Nei model for bHLH tree and 1,000 bootstrap replicates) was used to conduct the phylogenetic analyses of the conserved blocks.

Functional Analysis of Putative F3 5 H by Heterologous Expression in Saccharomyces cerevisiae
Multiple alignment with ClustalW was used to identify the consensus of flavonoid hydroxylase genes. Known and publicly available characterized protein sequences from various plant species coding for functional F3 H and/or F3 5 H were aligned. The obtained consensus sequence served as a template in local tblastn against RNA sequencing available for C. creticus fruit and leaves to identify candidate genes for hydroxylases. The primers used for the PCR amplification of putative full-length genes from Cistus cDNA were designed in Primer3Plus online tool and were for c15585 forward: ATGGTGGAAACACTGACTCCC and reverse: TTAGGAAACATAAGCACCCGGC. Leaves of the second stage (S2) of C. creticus (Falara et al., 2008) were used for total RNA isolation using Spectrum Plant total RNA kit (Sigma), according to the manufacturer's protocol. Subsequently, cDNA was synthesized with SuperScript III Reverse Transcriptase (Invitrogen) using a 1:1 mix of oligo(dT) 12−18 and random primers. PCR reaction was performed in a 5-µl final volume using Q5 high-fidelity polymerase (NEB, Ipswich, MA, United States), according to the following program: initial denaturation at 98 • C for 10 s, followed by 35 cycles including 98 • C for 10 s, 62 • C for 30 s, and 72 • C for 1 min and a single cycle of a 2-min final extension at 72 • C. pGEM T-Easy vector (Promega, Madison, WI, United States) was used for the cloning of the PCR-amplified product. After confirmation of selected clones by sequencing, the genes were subcloned in pYES2 yeast expression vector (Invitrogen) using the same set of primers, with the addition of EcoRI/XhoI restriction sites in their 5 ends. The resulting plasmids were transformed into yeast strain INVSc1 (Invitrogen) using LiAc method (Gietz et al., 1992), and positive yeast colonies were verified through PCR. A 5-ml SD-URA liquid pre-culture of a single transformed yeast colony was grown in a selective SD-URA plate containing 20 g l −1 glucose and was incubated at 30 • C and 200 rpm for 24 h. The propagated cells were collected and re-dissolved into 30-ml SD-URA medium containing 20 g l −1 galactose as a carbon source and for induction of protein expression. For functional characterization experiments, naringenin or dihydrokaempferol dissolved in dimethyl sulfoxide (DMSO) was separately added into the culture to a final concentration of 5 µM. The yeast culture was incubated at 30 • C for 36 h. A liquid/liquid extraction of the yeast products was performed by addition of 1:1 ethyl acetate (v:v), sonication for 15 min, and centrifugation for 5 min at 11,000 × g for three times. The ethyl acetate fractions were collected and evaporated in an EZ-2 ENVI Genevac (GeneVac, Ipswich, United Kingdom). The crude residue was re-dissolved in 150 µl of 80% methanol (v:v), filtered through 0.22-µm polytetrafluoroethylene membrane filters into 1.5-ml glass vials, and injected directly to ultra-performance liquid chromatography (UPLC)-DAD and/or UPLC-tandem mass spectrometry (MS/MS) instruments.

UPLC-PDA and UPLC-QTOF MS Analyses of Yeast Products
Yeast extracts were analyzed on a 1290 Infinity Binary UPLC (Agilent; Santa Clara, CA, United States) equipped with photodiode array (PDA) using a modified method described previously by Rafique et al. (2016). Separation of the compounds was achieved on a Machery and Nagel Nucleodur C18ec column (4.6 µm, 250 mm × 4 mm) set at 25 • C. The gradient consisted of two solvents containing 1% phosphoric acid: solvent A was water and solvent B was acetonitrile. The analysis was performed following the chromatographic conditions, namely: 100% A to 50% A in 25 min, plateau of 3 min, up to 100% A in 7 min, and final plateau of 5 min with a flow rate of 1 ml min −1 and monitoring at 222 and 280 nm. In order to identify naringenin, eriodictyol, dihydrokaempferol, dihydroquercetin, and dihydromyricetin, external standards were injected in a concentration of 1 mM.
A Waters Acquity UPLC coupled via an electrospray ionization (ESI) interface to a Synapt HDMS QTOF MS (Waters, Manchester, United Kingdom) operating in W-mode and controlled by MassLynx 4.1 was used. Both LC and MS parameters were previously described (Shahaf et al., 2013;Arapitsas et al., 2014). Yeast extracts were chromatographically separated using an ACQUITY UPLC 1.8 µm, 2.1 mm × 150 mm HSS-T3 column (Waters, Manchester, United Kingdom), thermostated at 40 • C. Mobile phase [0.1% (v/v) formic acid in water (A) and 0.1% in methanol (B)] was eluted with a flow rate of 0.28 ml min −1 , adopting the multistep linear gradient previously described by Arapitsas et al. (2014). The injection volume was 5 µl. Mass spectrometric data were collected in negative ESI mode over a mass range of 50-2,000 amu, with a scan duration of 0.4 s in centroid mode. The source parameters, as well as transfer cE and trap cE, were set as previously reported (Arapitsas et al., 2014). Annotation of compounds in yeast extracts was performed by comparing retention times and mass spectra (mass difference less than 5 ppm, isotopic distribution, and minimum three m/z ions) to those of the standards and based on internal database (Shahaf et al., 2013). In cases when standards were not available (pentahydroxyflavanone), tentative identification was made by using spectral features and literature data.

Statistical Analysis
Statistical analysis was performed using R Software (R Core Team, 2018) by applying the package stats for hierarchical clustering (HCA) and principal component analysis (PCA), gplots (Warnes et al., 2016) for heat map generation, corrplot for visualization of correlation matrices (Wei and Simko, 2017), and ggplot2 for data visualization (Wickham, 2016). HCA and PCA of metabolomics data were performed after cantering the data to 0 and scaling to unit variance. The expression data were not scaled or centered prior to performing the mentioned methods since these were already on a relative scale (log2 fold change). HCA was performed based on Euclidean distances with cluster agglomeration using Ward's (Ward, 1963) minimum variance method. When the dissimilarity matrix was based on correlation distances (1 -cor), cluster agglomeration was performed using complete linkage. The details on dissimilarities used for HCA are provided in the legends of the respective figures. Oneway analysis of variance was conducted for both metabolomics (absolute quantitates) and expression data (log2 fold change), followed by a Tukey post hoc test.

Identification of Phenolic Compounds in Cistus Fruit
To our knowledge, the current study is the first attempt to characterize the phytochemicals of the flowers and fruits of any Cistus species. Previous investigations have been primarily focused on the leaves or aerial parts (Vogt et al., 1987;Demetzos et al., 1989Demetzos et al., , 1990Danne et al., 1993;Kreimeyer et al., 1997;Pomponio et al., 2003;Barrajón-Catalán et al., 2011;Tomás-Menor et al., 2013;Papaefthimiou et al., 2014;Maggi et al., 2016). UHPLC-MS Orbitrap analysis identified 53 compounds in the flowers and fruits of C. creticus subsp. creticus, which belong to phenolic acids, flavonoids and their derivatives, and quinic acid ( Table 1). The peak numbers, compound names, molecular formulas, calculated and exact masses ([M−H] − , m/z), retention times (R t , min), mass accuracy errors, as well as major MS 2 , MS 3 , and MS 4 fragment ions are summarized in Table 1.
Qualitative metabolite analysis revealed that all 53 compounds were present in fruit stages S1-S3, while in Cistus flowers (S0) the following six compounds were not recorded: B-type PD dimer-isoform 3 (21), all three isoforms of B-type PC dimers (23, 25, and 32), aromadendrin 7-O-hexoside (24), and astragalin (40). Based on the obtained metabolite data, a flavonoid biosynthetic pathway in Cistus fruit has been proposed (Figure 1). The major flavonoid groups in Cistus fruit are flavones, flavonols, F3Os, and PAs. In Cistus fruit, at least three flavonoid biosynthetic branches starting from the flavanone intermediate naringenin (52) are present. One pathway results in the production of flavones apigenin (53) and luteolin (50), and the other proceeds via commonly occurring dihydroflavonols (dihydrokaempferol-46, dihydromyricetin-41, dihydroquercetin-42), leucoanthocyanidins (leucocyanidin and leucodelphinidin) to the anthocyanidins (cyanidin and delphinidin), which are further converted into F3Os and PAs, but also to some extent to anthocyanins as can be expected by the obvious pigmentation. Dihydroflavonols also give rise to flavonols kaemferol (40), quercetin (51), myricetin (49), and their glycosides in a side branch. It could be presumed that DFR shows distinct substrate specificity toward dihydroflavonols of the cyanidin and delphinidin branches while not converting dihydrokaempferol (46) as found in several plant species (e.g., Petunia), that LAR pathway overrides the ANR pathway in the biosynthesis of F3Os, and that catechin (12) is the predominant extension unit in oligomeric PAs in Cistus flowers and fruits. However, the dominant group of PAs in the analyzed fruit samples are PDs, indicating a tissue-specific expression of the flavonoid 3 ,5 -hydroxylase (F3 5 H). The PA amounts were generally decreasing during Cistus fruit development, as given in the heat map in Figure 1. The amounts of some F3Os (EC, C) and PAs have been previously found to decrease in the progression of fruit ripening in grape (Boss et al., 1996), bilberry (Jaakola et al., 2002), and strawberry (Schaart et al., 2013). Astringent persimmon fruits are rich in PAs even at maturity, while in nonastringent types the content decreases during development (Akagi et al., 2009a).
The peak eluting at R t = 0.46 min and displaying the deprotonated molecule [M−H] − at m/z 191 was identified as quinic acid (1) (Figure 2A and Supplementary Table 1). Although 1 is not a phenolic compound, it was interesting to trace its amount in Cistus fruits because it is involved in the regulation of the biosynthesis of aromatic compounds (Ghosh et al., 2012). This free acid is synthesized via the shikimate pathway and is abundant in a variety of fruits, such as papaya, pineapple, lemon, kiwi, cranberry, lingonberry, blueberry, apple, and orange (Jensen et al., 2002;Chinnici et al., 2005;Albertini et al., 2006;Erk et al., 2009;Hernández et al., 2009;Barboni et al., 2010;Ye et al., 2014). Significant amounts of 1 were found in S0 stage (flower) (∼12 µg 100 mg −1 DW), while it was significantly lower in fruit (stages S1-S3) (Supplementary Figure 4).
The peak visible at R t = 5.87 min, showing [M-H] − at m/z 301, was assigned as flavonol quercetin (51). Its amount was relatively stable in flowers and fruits of all stages, although it is usually further metabolized to different glycosides. The peak eluting at R t = 3.95 min, with pseudo-molecular ion [M−H] − at m/z 609, displayed MS2 fragmentation pattern (Supplementary Table 1) characteristic for the flavonol glycoside rutin (29). The content of 29 was the highest in flowers (S0 stage) but significantly reduced in all stages of fruit (S1-S3) (Supplementary Figure 4) (49), which eluted at R t = 4.76 min, was identified as deprotonated molecular ion [M-H] − at m/z 317. The amount of 49 was the highest in Cistus flowers (S0 stage) and S3 fruit but significantly decreased in S1-S2 fruit (Supplementary Figure 4). The major flavonol in Cistus fruit was myricetin 3-O-rhamnoside (31), reaching around 23 µg 100 mg −1 DW in flowers (S0). As in the case of 49, its amount was severely reduced in fruit (stages S1-S3) (Supplementary Figure 4). Therefore, the tissue-and stage-specific expression of the F3 5 H to yield the precursor dihydromyricetin for flavonol formation and a rhamnosyltransferase involved in the biosynthesis of 31 can be predicted. The peak corresponding to 31 and showing [M-H] − at m/z 463 eluted at R t = 4.01 min.
Apigenin (53) Table 1). The content of 50, 52, and 53 was not significantly different between Cistus flowers (S0) and fruits (Supplementary Figure 4). Some flavonoids, including 49, 51, 53, and kaempferol and their derivatives, are abundant in the leaves and resin of Cretan C. creticus subsp. creticus (Demetzos et al., 1989(Demetzos et al., , 1990). In the present study, kaempferol was not detected in rock-rose flowers and fruits, indicating its efficient and complete conversion by glycosyltransferases. Similarly, kaempferol was not recorded in the methanol extracts of C. creticus and C. monspeliensis leaves (Hickl et al., 2018). Although a few kaempferol glycosides are identified in Cistus flowers and fruits (Table 1), they were present in amounts which were below the limits of quantification of the analytical procedure and were thus not quantified.
The major F3Os in Cistus fruits are GC (5) and C (12), while C-3-O-gallate (37) and EC-3-O-gallate (38) were less abundant. C (12) was visible in the negative ionization mode as adduct with formic acid, which was used as the mobile phase. It showed pseudo-molecular ion [M−H + formic acid] − at m/z of 334 and was eluted at R t = 2.29 min (Supplementary Table 1), and its MS 2 fragmentation pattern was in accordance with some previous studies (Del Rio et al., 2004;Stöggl et al., 2004). Galloylated  Table 1). The contents of targeted F3Os varied during Cistus fruit development. The content of 12, 37, and 38 was the highest in fruits of S1 and S2 stages (Supplementary  Figure 4). The content of 5 was the highest in flowers (S0), with a concentration of ∼40 µg/100 mg −1 DW, and it decreased, during fruit development, down to ∼1.9 µg/100 mg −1 DW in S3 phase (Supplementary Figure 4). F3Os EC, C, EGC, and GC have been previously reported for C. incanus (Riehle et al., 2013).

Quinic acid (1) was very abundant in flowers (Supplementary
The quantitative content of targeted metabolites obviously changed during Cistus fruit development, as supported by PCA analysis (Figure 2B). Not surprisingly, the chemical profile of Cistus flowers (stage S0) was distinctively different from that of fruit of all developmental stages (S1-S3), with PC1 accounting for 43.2% of the total variance. The main contributors to PC1 are 31, 1, 29, and 14. On the other hand, fruit stage S3 segregates from S1 and S2 in PC2, explaining 25.51% of the data variance. Stages S1 and S2 showed no visible separation in PC1 and PC2. The highly contributing compounds to PC2 are 38, 37, 12, 49, 51, and 50. The PCA indicates that stages S1 and S2 are phytochemically closer to each other than to stages S0 and S3.
For a better perception on the phytochemical correlation among flowers (S0) and the three developmental fruit stages (S1-S3), data per compound normalized to 0-1 range, are presented as a heat map ( Figure 2C) with biological replicates (in rows) arranged according to HCA based on Euclidean distances and metabolites (in columns) organized according to HCA based on Spearman correlation distances (1 −cor sp , i.e., 100% positive correlation equals 0 and 100%, negative correlation equals 2). Flowers (stage S0) form a homogenous cluster, while fruits of different developmental stages are clustering together. However, samples belonging to stage S3 form a separate sub-cluster. On the other hand, HCA based on Spearman correlation distances (Figure 2C, top) provides a clear depiction of the targeted compound linkages. Two distinct clusters are visible. The first cluster contains the majority of F3Os (12, 37, and 38) and compounds 14 and 43. The second cluster is divided into two subclusters. The first sub-cluster consisted of flavonol aglycones and metabolites 1, 5, 29, 31, and 49, while all the rest of the compounds (2, 50, 51, 52, and 53) belonged to the second subcluster.
The expression of CcPAL1, CcPAL2, CcCHS1, and CcCHI is down-regulated during fruit development and senescence (Supplementary Figure 5) which is in agreement with the metabolomics data, showing the highest naringenin (52) levels in flowers (S0) and S1 Cistus fruit and its decrease during further fruit development. Comparable data were observed during development of the apple fruit (Henry-Kirk et al., 2012). Conversely, the levels of Cc4CL1 transcripts were relatively high in flowers and fruits of S1 and S2 developmental stages, and their expression decreased in S3 fruit (Supplementary Figure 5). For CcC4H, Cc4CL2, and CcCHS2, no apparent differences between samples were recorded, although the expression was slightly higher in flowers (S0) and S1 fruits. The enzyme 4CL converts p-coumaric acid to p-coumaroyl-CoA. This enzyme is concurrently engaged in controlling the efflux of p-coumaroyl-CoA in divergent branches of the phenylpropanoid pathway as well as in promiscuously converting other hydroxycinnamic acids (caffeic and ferulic acid) in the lignin biosynthesis. Therefore, the slightly different correlation pattern of 4CL2 expression in comparison with the other structural genes of flavonoid biosynthesis could indicate its additional role in lignification of Cistus fruit. The transcript levels of CcF3H1 and CcF3 5 H were relatively stable in flowers (S0) and in S1 and S2 fruits, and they significantly decreased in S3 fruit (Supplementary Figure 5). The expression of CcF3H2 is the highest in flowers (S0) and is down-regulated during fruit development, with the lowest transcript levels recorded in S3 fruit. A similar trend was observed for CcFLS gene, which is responsible for converting dihydroflavonols into flavonols. Thus, both CcF3H1 and CcF3H2 are active in flowers, while in fruit of early developmental stages only CcF3H1 is involved in the synthesis of flavonols and PAs. The expression of CcDFR, which reduces dihydroflavonols to leucoanthocyanidins, was detected throughout fruit development, and the transcript levels of both CcDFR1 and CcDFR2 reached their maximum in flowers (S0) and decreased during fruit development (Supplementary Figure 5). Two pathway branches are involved in the synthesis of F3Os, the LAR and ANR branches. In Cistus fruit, CcANR and CcLAr1 and also the intermediate CcANS expressions are relatively stable in flowers and S1 and S2 fruits but are decreased in S3 fruit. The transcript levels of a second LAR candidate, CcLAR2, are highest in flowers (S0) and slightly drop during fruit development and ripening (S1 and S2), reaching significantly lower amounts in ripe S3 fruit (Supplementary Figure 5).
In view of the above-mentioned condition, it could be presumed that flavonoid metabolism is differently modulated in flowers and fruits in a way to complement their morphology, physiology, and function. Higher transcript levels of LBGs (CcDFR1, CcDFR2, CcANR, CcANS, CcLAR1, and CcLAR2) and, at the same time, generally lower F3Os content in flowers, when compared to fruits, led to the belief that the flavonoid pathway in flowers is most likely directed toward anthocyanin biosynthesis, the main pigments of pink Cistus flowers, and/or toward PAs, both not quantified in this tissue within the present study. In fruits (S1-S3), the transcript levels of CcANR, CcANS, CcLAR1, and CcLAR2 followed the trend of F3Os content during their development and ripening. Based on metabolomic and transcriptomic data, it could be presumed that the LAR pathway is the predominant one so that the metabolic flux is directed toward the synthesis of 2R,3S-trans-flavan-3-ols (GC, C, and C-3-O-gallate). Although galloyled F3Os, such as GC, CG, and ECG, are recorded in Cistus fruits, the genes directly involved in their biosynthesis (CcECGT and CcCGT) have not been identified in the transcriptome of C. creticus fruit (stage S2) and were therefore not analyzed within the present study. During the maturation and ripening of bilberry (Zifkin et al., 2012), blackberry (Chen et al., 2012), and peach (Zhou et al., 2015), a decrease in the expression of ANR and LAR genes was observed similarly to the data obtained here. Furthermore, genes specific for the PA pathway, LAR and ANR, have been functionally characterized in a variety of fruit crops, such as grapevine (Bogs et al., 2005), persimmon (Ikegami et al., 2007;Akagi et al., 2009a,b), apple (Han et al., 2012;Henry-Kirk et al., 2012), strawberry (Schaart et al., 2013), and peach (Ravaglia et al., 2013). Besides the known function of LAR in the conversion of, e.g., leucocyanidin to (+)-catechin (Tanner et al., 2003), a new role in regulating the oligomerization and extension of PAs in Medicago truncatula has been proposed .
To explain the regulatory background of Cistus flavonoid metabolism during fruit development, we studied the expression profile of MYB, bHLH, and WD40 TFs. MYB proteins function as direct activators of structural genes and as activators of the gene(s) encoding bHLHs (Schaart et al., 2013). Regulation of F3O and PA biosynthesis may be conditioned by the feedback interactions of MYB and bHLH components of the MBW activation complex (Chezem and Clay, 2016). In Arabidopsis thaliana, the MBW complex formed by MYB, bHLH, and TTG activates the genes DFR, ANS, and ANR, the products of which coordinate the production of PAs in the seed coat (Nesi et al., 2001;Debeaujon et al., 2003). The regulation of F3O and PA biosynthesis and accumulation has previously been studied in fruits of several species, including persimmon (Ikegami et al., 2007;Akagi et al., 2009a,b), grape (Czemmel et al., 2009;Terrier et al., 2009), and apple (Henry-Kirk et al., 2012). In grape berries development, several MYBs (VvMYBPA1, VvMYBPA2, VvMYB5a, and VvMYB5b) specifically regulate PA synthesis (Deluc et al., 2008;Bogs et al., 2007;Terrier et al., 2009). PAspecific MYB regulators (DkMYB2 and DkMYB4) have also been described from persimmon fruit, which has unusually high PA levels (Akagi et al., 2009b(Akagi et al., , 2010. DkMYB4 was found to be the specific activator of DkANR but not of DkLAR (Akagi et al., 2010). To date, PA-related MYB activators have been identified in fruits of various plant species, such as PpMYBPA1 in nectarine (Ravaglia et al., 2013), MdMYB9/MdMYB11 in apple (Gesell et al., 2014;An et al., 2015), PpMYB7 in peach (Zhou et al., 2015), and PbMYB9 in pear (Zhai et al., 2015). Transgenic tomato lines expressing AtMYB12 TF of Arabidopsis under constitutive promoter exhibited an enhanced accumulation of flavonols in fruits, accompanied with the elevated expression of phenylpropanoid pathway genes involved in flavonol biosynthesis (Pandey et al., 2015). Generally, MYB factors are involved in primary and secondary metabolism; they regulate many physiological processes in plants such as cell fate and identity, development, hormone signal transduction, and response to environmental stresses (Dubos et al., 2010). To examine the patterns of PA regulation by MYBs, we traced the expression of eight MYB candidates identified in Cistus fruit transcriptome (CcMYB1, CcMYB2, CcMYB3, CcMYB4, CcMYB5, CcMYB6, CcMYB12a, and CcMYB12b). The expression profiles of these MYBs followed the same trend, being relatively stable in flowers (S0) and in S1 and S2 stage fruits, whereas these decreased significantly in S3 fruit (Supplementary Figure 5). To propose a putative specific function of identified C. creticus MYBs, a phylogenetic analysis was constructed for MYB TFs of different plants ( Figure 3A). CcMYB4 clusters closely with PA-related MYB activators, such as AtMYB123 (AtTT2) in Arabidopsis (Nesi et al., 2001), MdMYB9/MdMYB11 and MdMYB6 in apple (Gao et al., 2011;Gesell et al., 2014;An et al., 2015), and others. Anthocyanin and PA-related MYB activators require specific bHLH co-activators to work, and they contain a bHLH-binding domain in the N-terminal R3-MYB repeat (Ma and Constabel, 2019). The lignin-and flavonol-activating MYBs do not have this domain. Since Cistus CcMYB4 contains bHLH-binding domain, it could be presumed that it is most likely anthocyanin and PA biosynthesis activator. C. creticus CcMYB5 and CcMYB6 contain no bHLH-binding domain and are presumably involved in the activation of lignin and/or flavonol biosynthesis. They cluster close to the MYBs of the flavonol clade ( Figure 3A). The same goes for CcMYB12a and CcMYB12b which, according to some previous studies on A. thaliana and Solanum lycopersicum (Mehrtens et al., 2005;Ballester et al., 2010;Pandey et al., 2015), might play an important role in regulating the flavonol pathway. Cistus CcMYB1 and CcMYB2 are clustered closely to grapevine VvMYBC2L-1 Cavallini et al., 2015), which falls into the group of MYB repressors. However, based on the present study, CcMYB1 and CcMYB2 expression in Cistus flowers and fruits follows the decreasing expression patterns of biosynthetic genes during the fruit's development. Similarly, CcMYB3 is similar to apple MdMYB16 (Xu et al., 2017) and other R2R3-MYB repressors involved in the regulation of general phenylpropanoid and lignin biosynthetic pathway but is down-regulated during Cistus fruit development. Within the MYB phylogeny, most MYB repressors belong to the subgroup of R2R3-MYBs (Ma and Constabel, 2019), which is separated into two clusters, both containing the bHLH domain: (1) a general phenylpropanoid and lignin MYBs (CcMYB3) and (2) flavonoidrelated MYBs (CcMYB1 and CcMYB2). Only the anthocyanin and PA repressors interact with bHLH proteins. It has been suggested that, unlike most MYB flavonoid activators, MYB repressors seem to affect the biosynthesis of multiple flavonoids due to their bHLH-binding activity (Ma and Constabel, 2019). Activator MYBs are more specific than corresponding repressor MYBs. The mechanism of action of activator and repressor MYBs is largely determined by their competition for corresponding cofactor or DNA cis-element. Binding of flavonoid MYB repressors with bHLH co-activators can interfere with the MBW complex (Ma and Constabel, 2019). The relative abundance of MYB activator and repressor proteins in a given cell determines the incident of promoter activation and repression events. Functional characterization of identified Cistus MYBs would suggest their regulatory role within the flavonoid biosynthetic pathway in fruits.
The MBW complexes might include different classes of MYBs and bHLHs with specific functions in regulating the transcription of flavonoid biosynthetic genes. Within the present study, we also examined several Cistus bHLH candidates, and the majority of them (CcbHLH3, CcbHLH4, and CcbHLH5) showed a relatively stable expression in flowers and during the early development of fruits (S1 and S2 stages) and a decreasing trend in late fruit developmental stage (S3) (Supplementary Figure 5). The exceptions are CcbHLH1 and CcbHLH2 since their expression did not change during fruit development. It is not strange that some bHLH genes show observable different transcript behavior since they are regulating different branches of the flavonoid biosynthetic pathways and possibly other metabolic pathways. Some previous studies showed that the expression of 113 FabHLH genes in strawberry is dependent on the variety, fruit tissues/organs, and developmental stages (Zhao et al., 2018). By comparing the nucleotide sequences of Cistus bHLHs and bHLHs from other species, the phylogenetic tree was constructed ( Figure 3B). According to Montefiori et al. (2015), the bHLH TFs have been divided into two major groups, bHLH2/AN1/TT8 and bHLH/JAF13/EGL3 clades. Many members of both clades have been shown to regulate different branches of the flavonoid pathway (Montefiori et al., 2015). Cistus bHLH candidates CcbHLH1 to CcbHLH3 belong to the bHLH/JAF13/EGL3 clade (Figure 3B), together with grapevine VvMYC1 (Hichri et al., 2010), Arabidopsis bHLH proteins AtGL3 and AtEGL3 (Payne et al., 2000), apple MdbHLH33 (Espley et al., 2007), lychee LcbHLH1 (Lai et al., 2016), and petunia PhJAF13 (Quattrocchio et al., 1998), which are involved in the production of different flavonoids. The grape bHLH TFs VvMYC1 and VvMYCA1 were found to induce anthocyanin and PA production through the interaction with corresponding MYBs and consequent activation of the promoters of biosynthetic genes (Hichri et al., 2010). In lychee, LcMYB1 has been identified as the key regulator of anthocyanin biosynthesis (Lai et al., 2014). On the other hand, CcbHLH4 belongs to the bHLH2/AN1/TT8 clade, clustering close to strawberry FabHLH3 and putative negative regulator FabHLH3 (Schaart et al., 2013), apple MdbHLH3 (Espley et al., 2007), VvMYC1 (Hichri et al., 2010), Arabidopsis AtTT8 (Nesi et al., 2000), lychee LcbHLH2 (Lai et al., 2016), etc. Strawberry FabHLH3 and FabHLH3 are found to be involved in PA biosynthesis through the interaction with four MYBs (Schaart et al., 2013).
Similarly as MYBs and the majority of bHLHs, Cistus TTG candidate (CcTTG) showed a decreasing trend of expression during fruit development and ripening (Supplementary Figure 5). It was previously suggested that, in planta, TTG1 could regulate both the specific activity (i.e., interactions with other proteins or DNA) and the quantity (e.g., stability and localization) of the MBW complexes (Xu et al., 2015) and is therefore essential for flavonoid biosynthesis. Apart from its role in flavonoid biosynthesis, it is suggested that MYB-BHLH-TTG1 complex is modified by squamosal promoter-binding protein (SPBPs) to control trichome development and patterning in C. creticus leaves (Ioannidi et al., 2016). Overexpression of CcSPBPA/B Cistus genes in Arabidopsis leaves affected a broad range of plant developmental processes, most probably by directly binding to TTG1 (Ioannidi et al., 2016). Moreover, previous reports have shown that the elevated levels of SPBPs negatively modify TTG1-dependent physiological processes, including anthocyanin accumulation and trichome distribution of floral stems and leaves (Gou et al., 2011;Ioannidi et al., 2016). The expression patterns of the two SPBP candidates identified, CcSPBP1 (syn. CcSPLA) and CcSPBP2 (syn. CcSPLB), showed opposite trends in flowers (S0) and fruits (S1-S3). The expression of CcSPBP1 was increased in stages S1-S3 when compared to S0, while the level of CcSPBP2 transcripts gradually decreased during Cistus fruit development and was lowest in S3 stage (Supplementary Figure 5). These results indicate that CcSPBP1 is most likely an isoform involved in the regulation of flavonoid metabolism in fruit, while CcSPBP2 predominates in flowers. Negative correlations observed between CcTTG and CcSPBP1 further support the proposed role of CcSPBP1 as the MBW antagonist influencing flavonoid biosynthesis in rockrose fruit.

Ethylene Biosynthesis and Respiration During Cistus Fruit Development
Distinct hormone signaling pathways are known to be implicated in the regulation of the flavonoid pathway. In this context, the regulatory roles of abscisic acid (ABA) (Lacampagne et al., 2010;Zifkin et al., 2012), jasmonate (JA) Delgado et al., 2018), ethylene , brassinosteroids (Zhou et al., 2018), and other hormones have been proposed.
The highest ethylene production rate was observed in Cistus fruits of S1 and S2 stages, which is followed by a significant reduction in S3 stage. Respiration rates gradually decreased during fruit development, with the lowest value in S3 senescent fruits (Figure 4). The signals that trigger ripening in fruits are not clear; the ethylene production rates follow a similar pattern with the respiration rates, which suggests that fruits might behave as non-climacteric. Non-climacteric fruits, such as grape and strawberry, do not exhibit a sharp peak in respiration, although a rise in ethylene is sometimes observed . Obviously, non-climacteric fruits have minimum capacity to synthesize ethylene, which may influence some physiological and molecular events during development and ripening of this class of fruit . By contrast, in climacteric fruits, such as tomato, banana, or apple, the onset of ripening is marked by an obvious respiratory burst linked to ethylene action, which is the primary cue mediating and controlling most aspects of climacteric fruit ripening at the physiological, biochemical, and molecular levels Giovannoni et al., 2017).
Ethylene is synthesized from methionine (Figure 4) in a simple three-step biosynthetic pathway involving enzymes such as S-adenosyl methionine synthetase (SAMS), 1-amino cyclopropane-1-carboxylate synthase (ACS), and 1-amino FIGURE 4 | Ethylene production and respiration in Cistus fruits of three developmental stages (S1-S3). Bars with different letters are significantly different (p < 0.05) according to ANOVA post hoc Tukey's test. Ethylene is synthesized from methionine (A) in a simple three-step biosynthetic pathway involving the enzymes S-adenosyl methionine synthetase (SAMS), 1-amino cyclopropane-1-carboxylate synthase (ACCS), and 1-amino cyclopropane-1-carboxylate oxidase (ACO). S-adenosyl methionine, synthesized from methionine in a reaction catalyzed by SAMS, is converted by ACCS into ACC, and then ACC is oxidized by ACO to form ethylene. Ethylene further induces a signaling cascade which activates ethylene-responsive genes.
cyclopropane-1-carboxylate oxidase (ACO). S-Adenosyl methionine, synthesized from methionine in a reaction catalyzed by SAMS, is converted by ACS into 1-amino-cyclopropane-1-carboxylic acid (ACC) that then is oxidized by ACO to ethylene . In this study, we utilized available transcriptomic information from Cistus to identify ethylenerelated genes and known controllers of fruit development and ripening. Their expression in flowers and three stages of fruit development was investigated, aiming to better understand the maturation and ripening process in relation to flavonoid metabolism. We analyzed the expression patterns of genes involved in endogenous ethylene biosynthesis (CcSAMS, CcACS, CcACO1, and CcACO2) or in ethylene perception (ethylene response sensor, CcERS). The expression levels of CcSAMS were relatively stable during fruit development, with no significant changes from stages S0 to S3. Two ACO genes, CcACO1 and CcACO2, displayed different expression profiles, in which the expression of CcACO1 decreased in S3 stage, while the expression of CcACO2 was relatively stable (Supplementary Figure 5). This negative correlation between the expressions of the two CcACO genes was also depicted in the correlation analysis of gene expression shown in Figure 5. On the other hand, the expression of CcACS gradually decreased during fruit development. In tomato fruit, the expression of ACO and ACS increased during tomato fruit ripening, which has been linked to ripening-associated ethylene production in climacteric fruits (Nakatsuka et al., 1998).
Ethylene induces a linear transduction pathway, which starts with the hormone perception by a specific ethylene receptorethylene response sensor, which then initiates a signaling cascade by releasing the block exerted by serine-threonine protein kinase constitutive triple response 1 (CTR1) on ethylene insensitive 2 (EIN2). This further activates a transcriptional cascade involving ethylene insensitive 3 (EIN3)/ethylene insensitive-like 1 (EIL1) as the primary TF and then ethylene response factors (ERFs) which, in turn, regulate genes underlying ripening-related traits, such as color, firmness, aroma, taste, and post-harvest shelf-life Ju et al., 2012;Cheng et al., 2013;Liu et al., 2015). In the absence of ethylene, CTR1 inhibits the activity of EIN2 (Kieber et al., 1993). Cistus CcERS expression slightly increased in S3 fruit (Supplementary Figure 5).
Transcription factors play critical roles in relaying ripeninginducing signals and controlling ethylene biosynthesis and signaling. TFs such as the SPBP, colorless non-ripening, ERF, as well as various MADS-box genes regulate downstream ripening genes (Vrebalov et al., 2002;Liu et al., 2015). In the TF networks that affect fruit ripening, MADS-box family TFs serve as the central regulators of ripening and are involved in early fruit development, maturation, and pre-ethylene ripening events in both climacteric and non-climacteric fruits (Yao et al., 1999;Boss et al., 2001;Tadiello et al., 2009;Seymour et al., 2011;Liu et al., 2015;Giovannoni et al., 2017). Within the present study, we traced the expression of four putative MADSbox genes identified in Cistus fruit transcriptome (CcMbox1, CcMbox2, CcMbox3, and CcMbox4). The expression of CcMbox1 is significantly increased in S1 and S2 stages when compared to S0, and it decreased in S3 stage (Supplementary Figure 5). The expression of CcMbox2, CcMbox3, and CcMbox4 is relatively stable during the stages S0-S2, and it significantly decreased in S3 stage (Supplementary Figure 5). This is in accordance with the proposed regulatory role of MADS-box TFs during maturation and pre-ethylene events in fruits.

Chemometric Analysis of Gene Expression Data
A PCA combining gene expression studies revealed a clear separation between Cistus flowers and fruits ( Figure 5A). PC1 and PC2 cumulatively explained 92.6% of total variance. Stages S0-S2 are clearly separated along PC1 (87.1% variability) from S3. Furthermore, flowers (S0) and S3 fruits are separated from S1 and S2 along PC2 (5.5%). In general, stages S1 and S2 are very similar in terms of gene expression profiles. The genes highly contributing to the diversification of Cistus samples along PC1 are CcDFR2, CcPAL1, CcLAR2, CcSPBP2, CcFLS1, CcDFR1, and CcF3 5 H (Figure 5B). On the other hand, CcACO1, CcSPBP1, CcMbox3, and CcACO2 are the major variables contributing to the diversification of fruit developmental stages along PC2.
To analyze the transcriptional association among flowers (S0) and the three developmental stages of fruits (S1-S3), the data are also presented as a heat map of log2 fold changes (Figure 5C), arranged row-wise (samples) according to HCA constructed using Euclidean distances (Figure 5C, left) and column-wise (genes) according to HCA constructed using Pearson correlation distances (Figure 5C, top). Based on the presented tree ( Figure 5C, left), samples from stage S3 form a separate homogenous cluster, while samples from stages S0-S2 are clustered together. As for the targeted gene expression linkages, ethylene-related biosynthetic genes (CcERS and CcSAMS) and TFs CcBHLH1 and CcSPBP1 form a separate cluster (A). The expression of these genes is stable or slightly increases during fruit ripening. All the other genes cluster together (B) and are separated into two sub-clusters (Figure 5C, top). The first sub-cluster (b1) comprises CcACO1 and CcMbox1 FIGURE 6 | Correlation analysis of gene expression data. The correlation matrix was constructed using Spearman's rank-order correlation. Spearman's correlation coefficients are shown in the lower triangle, while the significance of the association samples using a two-sided test (p < 0.05) is denoted by asterisks in the upper triangle. The genes are ordered according to hierarchical cluster analysis (cluster agglomeration using complete linkage) constructed based on 1-corsp dissimilarity matrix (left tree).
to CcMbox4 TFs, while all the other genes belong to the second sub-cluster (b2). Within sub-cluster b2, two groups of genes could be diversified: (a) containing one EBG-Cc4CL2 and (b) comprising all the other genes divided into two subgroups: (1) with CcACO2 and (2) containing all the rest of the genes. Interestingly, the majority of EBGs (Cc4CL1, CcCHI, CcPAL2, CcC4H, and CcCHS2) cluster together with one LBG-CcANS and CcMYB5 TF. On the other hand, the majority of LBGs (CcF3H1, CcF3H2, CcFLS, CcLAR1, CcLAR2, CcANR, CcDFR1, and CcDFR2) cluster close to two EBGs (CcPAL1 and CcCHS1), five MYBs (CcMYB1-4 and CcMYB6), CcbHLH2 to CcbHLH5, CcTTG, and CcMbox4 gene. The presented associations between genes indicate that EBGs are most likely regulated by MYB4 to MYB6, while LBGs are under the control of MYB1 to MYB3 and corresponding bHLH TFs (CcbHLH2 to CcbHLH5).
A pairwise correlation analysis of the expression patterns between flavonoid biosynthetic genes and TFs was carried out to identify TFs that are co-expressed with the biosynthetic genes (Figure 6). Positive correlation has been observed between all the analyzed biosynthetic genes (EBGs and LBGs) and TFs (MYBs, bHLHs, and TTG), indicating the transcriptional regulation of flavonoid biosynthesis in Cistus fruit. The exception is the expression of CcbHLH1, which was negatively correlated with the analyzed biosynthetic genes, especially with Cc4CL2 and CcCH4. The positive correlation in the expression levels between MYB, bHLH, and TTG TFs indicated the role of MBW activation complex in the regulation of PA synthesis in fruits. MBW complexes might be involved in the developmental regulation of the flavonoid pathway at the transcriptional level, especially in the LBGs (DFR, ANR, ANS, LAR, F3H, F3 5 H, and FLS). EBGs (PAL, C4H, 4CL, CHS, and CHI) are probably less affected by the activities of these complexes. The exception again seems to be CcbHLH1, the expression of which was negatively correlated with the expression of other MYB, bHLH, and TTG TFs. Thus, it appears that PA biosynthesis in Cistus fruits is controlled by a network of MYBs (activators and repressors) and bHLHs, which could allow for a more subtle regulation of the transcriptional mechanism through different combinations of TFs within the MBW complex.
Among the analyzed ethylene-related genes, CcACO1, CcACS, and, to some extent CcACO2 were positively correlated with most of the flavonoid biosynthetic genes and TFs, with the exception of CcSPBP1 and CcbHLH1 (Figure 6). On the other hand, CcERS and CcSAMS were significantly negatively correlated with all the biosynthetic genes and majority of TFs but were positively correlated with CcSPBP1 and CcbHLH1. The only significant positive correlation of CcSAMS was with CcSPBP1.

Correlation Analysis Between Metabolomics and Gene Expression Data
The expression of structural genes is coordinately regulated and well correlated with metabolite pools, supporting the hypothesis that the biosynthesis of flavonoids during Cistus fruit development is controlled at the transcriptional level ( Figure 7A). All the structural gene transcripts were more abundant in flowers (S0) and during the first stages of fruit development (S1 and S2) when the majority of targeted flavonoids accumulated and the transcripts of the TFs were at a high level.
The exceptions are correlations between caffeic acid (14), quercetin 3-O-rhamnoside (43), and luteolin (50) content on one side and analyzed gene expressions on the other ( Figure 7A). Of the five bHLH candidates analyzed, CcbHLH1 was the only one positively correlated with the amount of 14 and 50, while both CcbHLH1 and CcbHLH2 were positively correlated with 43. The above-mentioned compounds were also positively correlated with CcSPBP1 and with some ethylene-related genes (CcERS and CcSAMS). However, such correlations should be considered with caution, especially in complex metabolic pathways in which the accumulation of the intermediates is balanced between the rate of biosynthesis and further utilization and conversion to other molecules. For example, the amount of naringenin (52) resulted from the activity of biosynthetic genes upstream of the pathway (PAL, 4CL, C4H, and CHI) but is also influenced by the rate of its conversion into other flavonoids as mediated by the activity of downstream enzymes (Figure 1). As stated before, 52 is considered the branching point for the biosynthesis of different flavonoid classes. In the same way, the detected concentrations of other intermediates of the pathway, including several detected aglyca, should be considered.
A correlation analysis of flavonoid-related biosynthetic genes, ethylene-biosynthetic genes, and TF expression patterns on the one hand and ethylene production/respiration on the other was carried out to identify possible associations between flavonoid and ethylene metabolism ( Figure 7B). Ethylene production rate and respiration were positively correlated with the expressions of CcACS and CcACO1 and negatively correlated with CcACO2, CcERS, and CcSAMS expressions. On the other hand, the expression of the majority of flavonoid biosynthetic genes and TFs analyzed was positively correlated with the level of ethylene production and respiration ( Figure 7B). Significant positive correlations were observed for endogenous ethylene production and CcDFR1, CcANS, CcLAR1, CcCHI, and CcPAL2 transcript levels in Cistus fruits. In strawberry, FabHLH98 is responsive to abiotic stress with the implementation of ABA and ethephon, which, according to Zhao et al. (2018), suggests its involvement in fruit ripening. Ethylene has strong effects on persimmon fruit ripening (Nakano et al., 2003) and on the removal of astringency resulting from high PA contents . ERFs induce a decrease of PAs in persimmon via the ethylene pathway . As previously reported, ethylene signaling is involved in the regulation of VvCHS, VvF3H, VvDFR, VvLDOX, and VvUFGT expression in grape (El-Kereamy et al., 2003). The transcription of VvCHI1, VvCHS3, VvDFR, VvLDOX, VvUFGT, and VvMYBA1 was up-regulated by ethephon, a commercial growth regulator that is quickly converted to ethylene, which resulted in the activation of PA and anthocyanin biosynthesis in grape berries (Liu M. Y. et al., 2016). Since ethylene production and respiration are positively correlated with the expression of the majority of flavonoid biosynthetic genes and related TFs, there are indications that ethylene biosynthesis and signaling might be involved in the biosynthesis of flavonoids in Cistus fruits. Further studies are needed in order to find connections between these two processes.

Isolation and Heterologous Expression of Flavonoid Hydroxylase (F3 5 H) in Yeast
Hydroxylation appears as one of the crucial steps defining the final metabolite pattern in Cistus flowers and fruits. Besides 3 ,4 -hydroxylated flavonols and F3Os, such as quercetin, catechin, and epicatechin derivatives including the respective PAs (i.e., procyanidins), 3 ,4 ,5 -hydroxlated counterparts (myricetin, gallocatechin, epigallocatechin, and its derivatives and related prodelphinidins) were found as major metabolites in the investigated tissues. These findings indicate the significant involvement of active hydroxylases responsible for the decoration of the basic flavonoid structure. To identify putative genes coding for functional F3 5 H proteins from C. creticus, a blast similarity search was performed. Functionally characterized genes from 42 plant species were selected and retrieved from NCBI public database (Supplementary Table 3). The local tblastn against RNA-seq from C. creticus fruit revealed two sequences for putative flavonoid hydroxylases, named c15585 and c24086, respectively. To assign a putative function to these genes, the retrieved candidate genes were used to generate a phylogenetic tree (Figure 8). The tree demonstrated that the putative contig c15585 is more closely related to F3 5 H genes from Theobroma cacao and Gossypium hirsutum and was clustering within the group of functionally characterized F3 5 Hs from various plant species. The candidate c24087 grouped clearly within the F3 H genes, branching again close to the same species (T. cacao and G. hirsutum). However, according to Seitz et al. (2006), a prediction of the in vivo function of flavonoid hydroxylases from phylogenetic sequence analysis could be, at least for candidates from Asteraceae family, misleading, and an enzymatic characterization of the encoded protein is suggested. Therefore, heterologous expression of the putative flavonoid 3 ,5 -hydroxylase in yeast and in vivo bioconversion of selected substrates were performed to determine the exact function of this important step in determining the substitution pattern of major Cistus metabolites. The yeast strain INVSc1 was used to transform and overexpress the candidate c15585, cloned into the expression vector pYES2. The induced transformed yeast cells were fed with two potential substrates, naringenin and dihydrokaempferol, dissolved in DMSO (Figure 1). The UPLC-PDA analysis of yeast FIGURE 7 | (A) Pairwise correlation analysis of gene expression and phenolics quantitative data in flowers (S0) and fruits of three developmental stages (S0-S3). Since different biological samples were used for metabolomics and expression measurements, the correlation was performed on averaged data by developmental stage. (B) Pairwise correlation analysis of gene expression and ethylene production and respiration (CO 2 ) data in Cistus fruits of three developmental stages (S1-S3). For both (A,B), the correlation matrix was constructed using Spearman (rank) correlations between compound quantity and relative gene expression. Asterisks denote significant (p < 0.05) association between paired samples using a two-sided test.
FIGURE 8 | Protein sequences for F3 H and F3 5 H in various plant species available in GenBank were selected, including those of Cistus creticus. Sequences were aligned using ClustalW and phylogenetic distance tree analysis followed in Geneious Tree Builder, with Jukes-Cantor genetic distance model and neighbor-joining method. A consensus tree was built using bootstrap resampling tree method, with 10,000 replicates and majority greedy clustering. The scale bar indicates 0.09 amino acid substitutions per site. Bootstrap values are shown in nodes. The accession numbers displayed are given in Supplementary Table 3. culture extracts after 24, 48, and 72 h of incubation revealed that the expression of c15585 gene produced a functional enzyme with the predicted activity, validated by the presence of the respective 3 ,4 ,5 -hydroxylated products pentahydroxyflavanone and dihydromyricetin (Figure 9). Due to a lack of authentic standard, pentahydroxyflavanone was tentatively identified by its retention time behavior and UV profile, which is characteristic for flavanones. In all extracts, the first expected reaction products, eriodyctiol and dihydroquercetin, were also detected. Furthermore, the identity of all six metabolites including pentahydroxyflavanone was confirmed by means of MS analysis using UPLC coupled via an ESI interface to a Synapt HDMS QTOF MS (Figure 9-C and Supplementary Figure 6). The mass accuracy error was, for all compounds, below 5 ppm. Compared with the respective standard, as available, the isotope similarity was above 70%, which allows the unambiguous identification of the given compounds (Figure 9-C). Detailed fragmentations of the compounds at two different energy levels are given in Supplementary Figure 6. These results clearly identify c15585, deposited in GenBank with number MT707661, as a functional F3 5 H involved in the specific hydroxylation of flavanones and/or dihydroflavonols, determining the final hydroxylation pattern of flavonols, F3Os, and PAs and most probably also of anthocyanins, which were not analyzed in this study.
Furthermore, the identified activity is in line with the phylogenetic clustering of the putative flavonoid hydroxylases' protein sequences. Therefore, we can postulate that the F3 5 H enzyme encoded by c15585 is involved in converting the metabolic flux toward biosynthesis of myricetin and its derivatives and through the delphinidin branch within the flavonoid biosynthesis network in C. creticus. Detection of products in methanolic extracts of yeast culture feeding assay with naringenin after biotransformation with F3 5 H gene. Ultra-performance liquid chromatography (UPLC)-photodiode array (PDA) chromatograms of (a) the substrate naringenin and the product eriodictyol and their corresponding absorption (b) and (c); (d) yeast culture transformed with CcF3 5 H flavonoid hydroxylase after 72 h of feeding with 5 µM naringenin and the alignment of the identified absorption of the products pentahydroxyflavanone (e) and eriodictyol (f) and the substrate naringenin (g) (line in blue) with the external standards from the diode array detector (DAD) library (line in red); (h) negative control yeast culture transformed with CcF3 5 H after 72 h of feeding with 5 µM naringenin without induction with galactose. (B) Detection of products in methanolic extracts of yeast culture feeding assay with dihydrokaempferol after biotransformation with F3 5 H gene. UPLC-PDA chromatograms of (a) the substrate dihydrokaempferol and the products dihydroquercetin and dihydromyricetin and their corresponding absorption (b), (c), and (d); (e) yeast culture transformed with CcF3 5 H flavonoid hydroxylase after 72 h of feeding with 5 µM dihydrokaempferol and the alignment of the identified absorption of the products dihydroquercetin and dihydromyricetin and the substrate dihydrokaempferol (h) (line in blue) with the external standards dihydroquercetin, dihydromyricetin, and dihydrokaempferol from the DAD library (line in red); (i) negative control yeast culture transformed with CcF3 5 H after 72 h of feeding with 5 µM dihydrokaempferol without induction with galactose. (C)

CONCLUSION
Being a rich source of flavonols, F3Os, and PAs, C. creticus subsp. creticus fruit represents a convenient model system to study the metabolism of these compounds in relation to development and ripening processes. Patterns of changes in the flavonoid content and expression profiles of biosynthetic genes and TFs during Cistus fruit development suggest a coordinated regulation of gene expression at the transcriptional level. The results propose a significant role of flavonoid hydroxylases in determining the content and ratio of major flavonoid compounds belonging to the groups of flavanones, flavones, dihydroflavonols, and flavonols and indirectly of F3Os and PAs. The function of the most important gene belonging to this group, CcF3 5 H, which is responsible for conducting the metabolic flux through delphinidin branches of the biosynthetic pathway, is confirmed. The involvement of a network of MYBs (activators and repressors) and bHLHs in the subtle regulation of the transcriptional mechanism through different combinations of TFs within the MBW complex is proposed. Although the results suggest the involvement of ethylene biosynthesis and signaling in the biosynthesis of flavonoids in non-climacteric Cistus fruit, further studies are needed in order to confirm the connections between these two processes.

DATA AVAILABILITY STATEMENT
The Data has been deposited to GenBank: GenBank MT707661.