Complex evolution of novel red floral color in Petunia.

Red flower color has arisen multiple times and is generally associated with hummingbird pollination. The majority of evolutionary transitions to red color proceeded from purple lineages and tend to be genetically simple, almost always involving a few loss-of-function mutations of major phenotypic effect. Here we report on the complex evolution of a novel red floral color in the hummingbird-pollinated Petunia exserta (Solanaceae) from a colorless ancestor. The presence of a red color is remarkable because the genus cannot synthesize red anthocyanins and P. exserta retains a nonfunctional copy of the key MYB transcription factor AN2. We show that moderate upregulation and a shift in tissue specificity of an AN2 paralog, DEEP PURPLE, restores anthocyanin biosynthesis in P. exserta. An essential shift in anthocyanin hydroxylation occurred through rebalancing the expression of three hydroxylating genes. Furthermore, the downregulation of an acyltransferase promotes reddish hues in typically purple pigments by preventing acyl group decoration of anthocyanins. This study presents a rare case of a genetically complex evolutionary transition toward the gain of a novel red color.

Introduction

Whether adaptation involves single genes of large phenotypic effect or proceeds through many genes with small individual effects is a crucial question in evolutionary biology (Orr and Coyne, 1992; Orr, 2005; Chevin and Beckerman, 2012; Barton et al., 2016; Boyle et al., 2017). In his theoretical work, Fisher predicted that mutations of small effect were most likely to produce phenotypes with increased fitness (Fisher, 1918, 1930). However, contemporaries of Fisher as well as more recent theory suggest that a single bout of adaptation could involve loci with a distribution of effect sizes, with both large-effect and small-effect mutations (Orr, 1998; Orteu and Jiggins, 2020). Today, there is abundant experimental evidence that the genetic basis of natural variation in individual traits can be extremely complex with many loci involved (Atwell et al., 2010; Chan et al., 2010; Turchin et al., 2012; Kooke et al., 2016; Guo et al., 2018; Sohail et al., 2019). At the same time, evidence for adaptation proceeding via few loci of large effect also abounds (Doebley, 2004; Hoekstra et al., 2006; Nadeau et al., 2016; Todesco et al., 2020). The relative importance of such large-effect mutations, however, remains contentious (Rockman, 2012) highlighting the need to examine genetic mechanisms with a critical eye.

A prime example of the relevance of mutations of large phenotypic effect has been pollinator-mediated selection on floral traits. Adaptation to shifts in pollinator availability is widely accepted to be a driving force in the rapid diversification of the angiosperms (Johnson, 2006; Sapir and Armbruster, 2010; Schiestl and Johnson, 2013; van der Niet et al., 2014). Color is a trait that can be easily be quantitated and monitored in different tissues during development. Many studies have demonstrated the importance of flower color for pollinators (Yuan et al., 2013a) and some have directly linked single genes to pollinator preference (Hoballah et al., 2007; Hopkins and Rausher, 2011; Yuan et al., 2013b; Sheehan et al., 2016; Kellenberger et al., 2019).

Anthocyanins are the major floral pigments in the angiosperms and are produced by the flavonoid pathway (Winkel-Shirley, 2001). The two most common branches of the flavonoid pathway are the anthocyanins (red, purple, and blue pigments responsible for visible color) and flavonols (responsible for ultraviolet [UV] color). Flavonoids are synthesized as a part of the complex metabolic network of phenylpropanoids, which includes a large variety of primary and secondary compounds, such as lignins, volatile signals, developmental regulators, and defense compounds (Winkel, 2006; Yang et al., 2017). Extensive knowledge of flavonoid pathway biosynthesis and regulation provides a foundation to study process of floral color evolution.

In floral color adaptation there are three general types of phenotypic transitions as follows: loss of color, shifts in color hue, and gain of color (Rausher, 2008). Losses of pigmentation tend to be relatively simple, with loss-of-function mutations in single genes of major phenotypic effect. Between species, these mutations typically occur in transcription factors which leads to downregulation of anthocyanin biosynthetic genes (Quattrocchio et al., 1999; Schwinn et al., 2006; Hoballah et al., 2007; Lowry et al., 2012; Streisfeld et al., 2013; Yuan et al., 2013b; Esfeld et al., 2018). Shifts in color hue, which are typically from blue–purple to red, always involve a change in anthocyanin hydroxylation (in the case of blue–purple to red, a decrease in hydroxylation) and often also include transcription factors. Thus the genetic mechanisms of color shifts are more diverse, through both loss-of-function mutations and expression changes in biosynthetic and regulatory genes (Zufall and Rausher, 2004; Streisfeld and Rausher, 2009; Des Marais and Rausher, 2010; Hopkins and Rausher, 2011; Smith and Rausher, 2011; Wessinger and Rausher, 2015).

It has been suggested that floral color gains are frequent at macroevolutionary scales (Smith and Goldberg, 2015). Land plants have the ability to synthesize flavonoid and anthocyanin pigments (Campanella et al., 2014); given that floral color is an evolutionarily labile trait, the likelihood of several instances of floral color gains in the evolutionary history of a genus is considerable (i.e. Armbruster, 2002). The rarest case of floral color evolution is the regain of floral color in a lineage which has already lost color. Obviously, it is easier to break something than to fix it and with increasing time, it becomes more difficult to re-evolve a complex trait that is lost, as first stated in Dollo’s law (Dollo, 1893; Gould, 1970). Indeed, demonstrations of molecular–genetic mechanisms underlying interspecific floral color gains are not well represented in the literature (but see Cooley et al. (2011)), and most documented color gains are within-species polymorphisms (i.e. Streisfeld et al., 2013).

A genus with losses, shifts, and regain of floral color is Petunia (Solanaceae). The garden petunia, Petunia hybrida, is a horticultural hybrid with a remarkable diversity of colors and color patterns and a long history of research on the chemistry and genetics of anthocyanin biosynthesis (Figure 1; Koes et al., 2005; Quattrocchio et al., 2006a; Tornielli et al., 2009; Bombarely et al., 2016). The naturally occurring species of the genus Petunia are native to South America and have undergone several shifts in pollination system, resulting in two main clades, the short-tube and long-tube clades, so named based on their floral tube length (Reck-Kortmann et al., 2014). The species in the short-tube clade are bee-pollinated and have the ancestral purple flowers, whereas the species of the long-tube clade are diverse in flower color, scent, and morphology and are visited by different pollinators (Sheehan et al., 2012; Hermann et al., 2013; Reck-Kortmann et al., 2014). Petunia axillaris presents white UV-absorbent flowers pollinated by hawkmoths (Stehmann et al., 2009), P. secreta purple UV-reflective flowers pollinated by solitary bees (Stehmann and Semir, 2005; Rodrigues et al., 2018a,b), and P. exserta red UV-reflective flowers pollinated by hummingbirds (Lorenz-Lemke et al., 2006). These three species grow in the Serra do Sudeste region in the south of Brazil, where both P. secreta and P. exserta are strict endemics.

Figure 1

Petunia exserta uses purple pigments to create a red flower color. A, Schematic representation of the Petunia flavonoid biosynthetic pathway, highlighting the end product groups flavonols and anthocyanidins, and the general chemical structure of the six common anthocyanidins. Flavonoids are classified into three groups according to their B-ring hydroxylation: monohydroxylated, dihydroxylated, and trihydroxylated, highlighted in yellow in figure diagram. The dihydroflavonols DHK, DHQ, and DHM are common precursors to both flavonols and anthocyanidins (anthocyanin aglycones). Enzymatic modification of phenylpropanoid precursors by F3H creates anthocyanidins and flavonols that are monohydroxylated on the flavonoid B-ring, while F3′H (encoded by HT1) activity creates dihydroxylated and F3′5′H (encoded by HF1 and HF2) trihydroxylated pigments; F3′5′H can also act upon DHK directly. FLS modifies dihydroflavonols into the flavonols kaempferol, quercetin, and myricetin, in increasing B-ring hydroxylation order. Anthocyanins are produced from dihydroflavonols by the sequential action of DFR and ANS. Further modifications by a suite of six enzymes yield five anthocyanin types, in increasing B-ring hydroxylation order: the brick-red to orange pelargonidin, the red to magenta cyanidin and peonidin, as well as the blue to purple delphinidin, petunidin, and malvidin. The enzyme FLS has low activity on DHM in Petunia, triggering reduced presence of myricetin. Increased hydroxylation of anthocyanidins through the action of F3′H/HT1 and F3′5′H/HF shift color from red toward purple–blue. Glycosylation, acylation, and methylation of the anthocyanins has a similar effect. We group the biosynthetic genes/enzymes in the schematic into the following: middle biosynthetic genes (hydroxylating HT1/F3′H, HF1/F3′5′H, and HF2/F3′5′H, anthocyanin-specific DFR and ANS, flavonol-specific FLS) and late biosynthetic genes (decorating 3GT, ART, AAT, 5GT and methylating MT/3′AMT, MF1/3′5′AMT, MF2/3′5′AMT). Anthocyanidins are differentially substituted on the anthocyanidin B-ring as follows: Pelargonidin R1=R2=H, Cyanidin R1=OH, R2 =H, Delphinidin R1= R2 =OH, Peonidin R1=OMe, Petunidin R1=OMe, R2 =OH, Malvidin R1= R2=OMe. B, Anthocyanidins (anthocyanin aglycones) and flavonol concentrations of P. exserta petal hydrolyzed extracts, with visible and UV photos of P. exserta. C, The LC–UV chromatogram of hydrolyzed extract of P. exserta anthocyanidins overlaid with a chromatogram of a mixture of reference standard anthocyanidins, indicates the presence of delphinidin, cyanidin, petunidin, peonidin, and malvidin as a shoulder peak, but no pelargonidin is detected; n = 9, error bars ± sd. D, Light microscopy of adaxial epidermal peel of P. exserta petal limbs (top and bottom left panels, overhead and side views respectively) shows vacuoles containing anthocyanins but no carotenoid-containing chromoplasts. In contrast, an epidermal peel of the adaxial (inner) floral tube of P. exserta shows cells with yellow carotenoid-containing chromoplasts clustered around the edges of the cells, with anthocyanins in the vacuole from a single anthocyanin-pigmented vein cell. E, The species tree of the four focal Petunia species as estimated in Esfeld et al. (2018). F, In Petunia, the enzyme DFR does not accept DHK as a substrate, thus preventing the biosynthesis of pelargonidin. Amino acids from the region of the DFR protein thought to be involved in substrate specificity as defined by Johnson et al. (2001) are in bold, with the most important residue 143D starred; all Petunia species have identical sequences in this region. A full and detailed alignment against the crystal structure is shown in Supplemental Figure S1. Abbreviations: For genes and enzymes that have different names, we state the Petunia gene name after a forward slash. DHK, dihydrokaempferol; DHQ, dihydroquercetin; DHM, dihydromyricetin; F3H, flavanone 3-hydroxylase; F3′H/HT1, flavonoid-3′-hydroxylase, encoded by HT1; F3′5′H/HF, flavonoid-3′5′-hydroxylase, encoded by HF1 and HF2; FLS, flavonol synthase; DFR, dihydroflavonol-4-reductase; 3GT, anthocyanin-3-glucosyltransferase; ART, anthocyanin rhamnosyltransferase; AAT, anthocyanin-3-rutinoside acyltransferase; 5GT, anthocyanin 5-glucosyltransferase; 3′AMT, 3′-anthocyanin methyltransferase encoded by MT; 3′5′AMT, 3′5′-anthocyanin methyltransferase encoded by MF1 and MF2. Chemical abbreviations: Pel, pelargonidin; Cy, cyanidin; Del, delphinidin; Peo, peonidin; Pet, petunidin; Mal, malvidin.

The most recent phylogeny robustly places the white P. axillaris as sister to the two colored species P. exserta and P. secreta (Figure 1E;Esfeld et al., 2018). The phylogenetic history of AN2 and MYB-FL—the key R2R3-MYB transcription factors that determine the spatial and temporal expression of the anthocyanin and flavonol biosynthetic pathways in Petunia—further supports the notion that P. secreta and P. exserta regained floral pigmentation from a colorless ancestor (Sheehan et al., 2016; Esfeld et al., 2018). Notably, the independent regains of color have resulted in two different colors: the ancestral purple in P. secreta and a red color unique to the genus in P. exserta. At the molecular level, P. secreta regained floral color by a compensatory deletion in the AN2 coding sequence that restored the AN2 reading frame. This gain-of-function mutation to a transcription factor is an “easy fix” to a difficult problem: to restore the synthesis of the ancestral purple color without additional changes to the essential biosynthetic genes (Esfeld et al., 2018).

Compared to the surprisingly simple pseudogene resurrection in P. secreta, reacquisition of color in P. exserta must be inherently more complex. First, all P. exserta accessions studied (Esfeld et al., 2018) contain one or more frameshifts in the AN2 coding region, and thus are unlikely to encode full-length functional proteins. Consequently, a substitute transcription factor must have been recruited. Second, additional modifications are required to obtain the distinct red color. Therefore, P. exserta represents both a regain in floral color and a transition from purple to red, a new color hue in the genus. Red color can be achieved through different biochemical means, the main strategies being changes in anthocyanins, synthesis of carotenoids, alterations to vacuolar pH, or through any combination (Ng and Smith, 2016b). In P. hybrida, the key anthocyanin biosynthetic enzyme DFR has lost the ability to synthesize the orange–red pelargonidin anthocyanins (Johnson et al., 2001). This may also be true for its wild ancestors. Previously published pigment analyses do not discuss carotenoids and disagree on the composition and relative abundance of anthocyanins present in P. exserta flowers. Specifically, it is unclear whether the two major types of anthocyanins consist of pelargonidin and cyanidin (Griesbach et al., 1999), or cyanidin and delphinidin (Ando et al., 1999; Ando et al., 2000). These issues need to be resolved to understand how P. exserta makes red flowers. Here we present a chemical analysis of the floral pigments of P. exserta and an in depth molecular–genetic analysis, as well as functional validation. These led us to three candidate transcription factors and at least three candidate biosynthetic genes, which were validated by functional analysis. We conclude that the evolution of red floral color in P. exserta flowers proceeded through multiple and subtle genetic alterations.

Results

Hydroxylation, methylation, and acylation modifications to P. exserta anthocyanins create the red color

To characterize the biochemical aspects of the red color of P. exserta, we analyzed the flavonoid pigment profile in hydrolyzed petal limb extracts. Acid hydrolysis removes the O-glycoside and acyl decorations, but not B-ring hydroxylation or methylation, to reveal the flavonoid backbone (Harborne, 1998), here resulting in anthocyanidins (anthocyanin aglycones) and flavonol aglycones. Cyanidin and delphinidin were the most abundant anthocyanidins (44% cyanidin, 40% delphinidin, 8% petunidin, 4% malvidin, 3% peonidin; Figure 1, A and B). Notably, the orange–red pelargonidin was not detected (Figure 1C). Flavonol concentrations were low, as previously described in Sheehan et al. (2016). As fully decorated cyanidin and delphinidin typically produce magenta to blue hues, we proceeded to identify the anthocyanin compounds (anthocyanin aglycone and attached decorative moiety) in nonhydrolyzed extracts of P. exserta limbs. We identified four major compounds as follows: a delphinidin diglycoside, cyanidin diglycoside, petunidin diglycoside, and peonidin diglycoside, with traces of additional malvidin anthocyanins (all likely rutinosides, 6-O-α-l-rhamnosyl-d-glucose; Supplemental Figure S1). When UV–Vis (ultraviolet to visible) spectra were observed in the LC–MS (liquid chromatography – mass spectrometry) experiment (at low pH), the absorption maxima of these anthocyanins was between 519 and 522 nm, that is, red-shifted relative to purple-shifted maxima of 530–550 nm (Supplemental Figure S1B; Harborne, 1958; Mabry et al., 1970; Markham, 1982). In comparison, purple-colored P. secreta and P. inflata produce several types of acylated anthocyanins (acyl moieties are aromatic or aliphatic acids, such as p-coumaric acid; Supplemental Figure S2). In addition, the purple Petunia species produce methylated anthocyanidins petunidin and malvidin almost exclusively (Esfeld et al., 2018). We conclude that P. exserta does not produce monohydroxylated pelargonidin anthocyanins and that the anthocyanins produced are severely deficient in acylation and methylation.

Petunia exserta does not produce yellow–orange carotenoid copigments in the petal limb epidermis, although it does produce them in the inner floral tube epidermis (Figure 1D). Thus, P. exserta has the biosynthetic machinery to synthesize carotenoids but they do not contribute to the red hue of the petal limb. Neither does P. exserta have an especially low petal homogenate pH compared to sister taxa (Supplemental Figure S3), nor to those reported in various wild-type P. hybrida varieties (de Vlaming et al., 1983; Quattrocchio et al., 2006b). Although we cannot rule out an interaction between particular anthocyanins and vacuolar pH, we note that P. exserta and P. secreta have approximately the same petal homogenate pH. Given the major color difference between the two species (red versus purple), we conclude that the major biochemical basis of color differentiation between the two species is that of anthocyanin composition rather than vacuolar pH. Thus, P. exserta’s red color is due to its relatively simple flavonoid composition, highlighting changes in anthocyanin hydroxylation, methylation, and acylation. This indicates that protein function and/or expression of F3′H/HT1, F3′5′H/HF1, F3′5′H/HF2, AAT, 3′AMT/MT, and 3′5′AMT/MF might be compromised (many of the Petunia genes have specific names; we have added these after the forward slash, abbreviations in Figure 1A).

No evidence for deficiencies in most anthocyanin biosynthetic protein sequences

To assess potential functional divergence of proteins and genes, sequences of biosynthetic genes F3′H/HT1, F3′5′H/HF1, F3′5′H/HF2, FLS, DFR, ANS, 3GT, ART, AAT, 5GT, 3′AMT/MT, 3′5′AMT/MF1, and 3′5′AMT/MF2 were compared. With the exception of MF1 and MF2, all flavonoid pathway genes in the four Petunia species in this study had low pairwise divergence, with ≤1.5% nucleotide divergence among the long-tube clade species (P. exserta, P. axillaris, P. secreta). Protein coding sequences showed no loss-of-function mutations in the four species examined, with two exceptions. Petunia exserta 3′5′AMT/MF1 harbors a frameshift mutation in the fifth exon leading to an early stop codon. Petunia axillaris features a different frameshift mutation in 3′5′AMT/MF1, and the gene cannot be found in P. secreta. Although P. exserta has a seemingly functional 3′5′AMT/MF2, both P. axillaris and P. secreta have nonsense mutations in 3′5′AMT/MF2. The low abundance of malvidin in P. exserta could be explained by an MF1 pseudogene, with residual product from 3′AMT/MT which can have anthocyanin 5′-O-methylation activity (Provenzano et al., 2014).

The complete absence of the orange–red pelargonidin suggests that the P. exserta DFR enzyme cannot reduce the monohydroxylated precursor DHK, as is the case in P. hybrida (Johnson et al., 2001). Indeed, the DFR amino acid sequences in the defined region of substrate specificity are identical between P. hybrida and the wild species, including P. exserta (Figure 1F;Supplemental Figure S4A; Johnson et al., 2001; Petit et al., 2007). We additionally examined conserved motifs and known active sites previously described in the literature for F3′H/HT1, F3′5′H/HF1, F3′5′H/HF2, and AAT (Nakayama et al., 2003; Seitz et al., 2007). All four Petunia species had identical active site sequences in these regions for each hydroxylating gene, altogether suggesting that P. exserta does not contain unique mutations in flavonoid biosynthetic loci that affect enzyme function. There are, however, significant amino acid changes in substrate recognition sites 1 and 2 between the F3′5′H genes HF1 and HF2 (Supplemental Figure S4B). This suggests a potential functional difference between these two F3′5′H paralogs. It follows that differential expression (DE) of the two HF copies could result in a different flavonoid composition. In summary, the in silico analysis of the protein sequences of all but one of the biosynthetic genes makes it plausible that they encode active proteins with conserved function in Petunia; the duplicate F3′5′H/HF proteins may recognize different substrates.

Anthocyanin biosynthetic gene expression is restored in P. exserta

With the exception of the MF genes, all of the flavonoid biosynthetic genes encode functional proteins. We therefore compared gene expression in developing petal limbs between red P. exserta, white P. axillaris, and purple P. secreta with reverse transcription quantitative PCR (RT-qPCR; Figure 2A;Supplemental Table S1) and added additional purple species P. inflata with RNAseq (Supplemental Figure S5 and Table S2). Expression levels of phenylpropanoid and flavonoid early biosynthetic genes were generally lower in P. exserta than in P. axillaris, and similar to expression levels in P. secreta and P. inflata (Supplemental Figure S5 and Table S2). The first two committed anthocyanin biosynthetic genes, DFR and ANS, were higher in P. exserta than in P. axillaris but remained lower than in P. secreta, 2.4-fold for DFR and 1.4-fold for ANS (Figure 2A). We observed only moderate values (<0.65) of allele-specific expression (ASE) as measured in the F1 hybrid of P. axillaris × P. exserta for most of the flavonoid pathway biosynthetic genes, indicating that the DE observed is at least partly due to transregulatory effects (Figure 2B).

Figure 2

Expression of P. exserta flavonoid biosynthetic pathway genes differ in comparison to the white P. axillaris and the purple P. secreta, and ASE in the F1 hybrids between P. axillaris and P. exserta. A, Relative expression as obtained by RT-qPCR of the biosynthetic genes of the flavonoid pathway conducted on stage 4 petal limb tissue. Expression is represented as the mean ratio of the gene relative to expression of the reference gene SAND. Box plots show medians of relative expression from three biological replicates, error bars ± sd. Statistics were calculated using aligned ranks transformation ANOVA with Tukey post hoc comparisons; letters indicate significantly different groups and full statistics are given in Supplemental Table S1. Stage 4 values for P. axillaris and P. secreta were previously published in Esfeld et al. (2018). B, ASE conducted on three biological replicate F1 hybrids between P. axillaris and P. exserta (representative photo of F1 hybrid shown in legend). Genes only shown if ASE could be calculated (presence of SNPs). For ASE, mapped RNAseq reads were realigned and averaged over SNPs between the parental alleles and counted (see “Materials and methods”), bars depict the mean read counts over all SNPs per gene, per species; points depict mean read counts per individual biological replicate. ASE, calculated as the major allele frequency ratio, is written above bars, error bars ± se. We use a strict threshold of ASE values greater than or equal to 0.75 (indicating at least 3:1 ratio of allelic imbalance) to indicate cis-regulation of the alleles, ASE values less than this indicates regulation in trans (e.g. ASE of 0.50 indicating a 1:1 ratio of allelic expression and equal activation by a trans factor). **P < 0.01, ***P < 0.001

Changes in flavonoid hydroxylation gene expression associated with shift to novel red color

Expression of F3′H/HT1 was very high in white, UV-absorbing P. axillaris, in agreement with the high levels of the dihydroxylated flavonol quercetin in this species (Figure 2A;Supplemental Figure S5). Expression of F3′H/HT1 is reduced in the UV-reflective P. exserta and P. secreta related to P. axillaris due to the absence of MYB-FL activity as well as cis-acting mutations (Figure 2, A  and B; Supplemental Figure S5; Sheehan et al., 2016; Esfeld et al., 2018). Comparison of F3′H/HT1 expression between P. exserta and P. secreta showed a 4.5-fold higher expression in P. exserta, promoting a shift from trihydroxylated toward dihydroxylated anthocyanins (e.g. cyanidin).

The F3′5′H/HF1 and F3′5′H/HF2 proteins synthesize trihydroxylated delphinidin, petunidin, and malvidin. That P. exserta requires F3′5′H at all is in line with the high concentration of delphinidin. We observed an unexpected pattern in the HF loci. Expression of both genes is low in P. axillaris and high in P. secreta. However, whereas F3′5′H/HF1 remained low in P. exserta, F3′5′H/HF2 was elevated to the expression level of P. secreta (Figure 2A;Supplemental Figure S5C). Interestingly, in slightly older petal limbs differences became more pronounced; P. secreta increased F3′5′H/HF1 expression while P. exserta increased F3′5′H/HF2 expression (Supplemental Figure S6). Thus, P. exserta preferentially expresses HF2 and P. secreta HF1. Together with divergence in the substrate recognition sites between HF1 and HF2 (Supplemental Figure S4B), this suggests that F3′5′H/HF1 and F3′5′H/HF2 have different properties. In summary, the higher F3′H/HT1 expression as well as the relative difference in expression between F3′5′H/HF1 and F3′5′H/HF2 between red P. exserta and purple P. secreta could lead to the difference in phenotype of dihydroxylated versus trihydroxylated anthocyanins.

Reduced expression of anthocyanin late biosynthetic genes prevents acylation and most methylation

Given the lack of acylated anthocyanins present in P. exserta (Supplemental Figure S1), we reasoned that one or more downstream anthocyanin modifying genes could be downregulated. The 3GT gene is prerequisite to obtain simple glucosylated anthocyanins (e.g. 3-O-glucosides), ART to obtain rhamnosylated anthocyanins, and AAT to obtain acylated anthocyanins. Levels of expression for 3GT and AAT in P. exserta were comparable to those of P. axillaris but importantly, much lower than in P. secreta (2.8× and 13.5× less, respectively; Figure 2; Supplemental Figures S5 and S6). We observed only moderate ASE for AAT (0.62) in the P. axillaris × P. exserta F1 hybrid, indicating regulation primarily in trans. Petunia exserta and P. axillaris have identical 3GT sequences, so no ASE could be calculated. However, given that the anthocyanins observed were glycosylated yet distinctly lacking in acylation, we conclude that the downregulation of AAT is an essential aspect of the P. exserta phenotype.

As anthocyanin biosynthesis in Petunia is thought to proceed in a sequential fashion due to strict substrate specificity, downregulation of AAT should prevent the majority of flux from downstream modifications (i.e. methylation; Tornielli et al., 2009). However, 8% of P. exserta anthocyanidins are petunidin, which is singly methylated on the B-ring (Figure 1A) by 3′AMT/MT, indicating some nonlinearity in the pathway. Furthermore, we observed moderate upregulation of 3′AMT/MT in P. exserta (Figure 2A) as well as strong ASE (0.98, biased to P. exserta copy) in the P. axillaris × P. exserta F1 hybrid indicating cis-regulatory effects. Further methylation to produce malvidin is hindered as 3′5′AMT/MF1 is a pseudogene in P. exserta (as well as P. axillaris), and although expression of 3′5′AMT/MF2 is curiously high in P. exserta compared to sister species, P. exserta only makes small amounts of malvidin (3%) (Figures 1B, 2A; Supplemental Figures S5 and S6). Thus, the overall contribution of the methylating loci to the P. exserta phenotype is small. Taken together, these data point to 3GT and AAT as promising anthocyanin modifying candidate genes warranting further investigation.

Identifying candidate replacement MYB transcription factors for anthocyanin biosynthesis

We previously showed that a frameshift mutation in the P. exserta MYB-FL gene underlies a major QTL (quantitative trait locus) for loss of UV-absorbing flavonols (Sheehan et al., 2016). A new analysis revealed that this QTL is also responsible for a gain of anthocyanin production (Figure 3A;Supplemental Table S3). MYB-FL belongs to MYB family subgroup (SG7) and is not an anthocyanin biosynthesis activator. Therefore, the loss of MYB-FL function in P. exserta may enable an increase in anthocyanin production by a trade-off between the anthocyanin and flavonol branches of the flavonoid biosynthesis pathway. Three minor QTLs on chromosomes 1, 3, and 7 exclusively affected anthocyanin production (Supplemental Tables S3 and S4).

Figure 3

QTL analysis and flavonoid pathway-related MYB SG6 and “G20”: gene tree, expression in stage 4 petal limbs as obtained by RNAseq, and ASE in the F1 hybrids between P. axillaris and P. exserta. A, QTL analysis of anthocyanin and flavonol content in an F7 RIL (recombinant inbred lines) population of P. axillaris × P. exserta. QTLs are shown with allelic effects; LOD values greater than zero indicate P. exserta alleles (“BB”) affect the phenotypic variance and below zero indicate P. axillaris alleles (“AA”) affect the phenotypic variance; this calculated by multiplying the LOD value by the additive effect divided by the absolute value of the additive effect. A major QTL is located on chromosome 2, with minor QTL on chromosomes 1, 3, and 7. Percent variance estimated (PVE) listed per chromosome. B, Phylogenetic relationship of MYB proteins belonging to SG6 and the PH4 clade (“G20”) from P. exserta, P. axillaris, and P. secreta. The complete MYB tree for P. exserta and P. axillaris is in Supplemental Figure S7. MYB SG6 and G20 form monophyletic groups within the MYB tree. Support values shown are based on 1000 bootstraps from a RAxML maximum likelihood analysis. AN4 has two canonical copies in long-tube clade species, denoted as either “1” or “2”; P. exserta AN4-2 is a pseudogene (denoted b). Additional AN4-like genes (denoted a) were discovered in P. exserta but not in either P. axillaris or P. secreta; one copy was excluded due to an early nonsense mutation yielding a protein of only seven amino acids (Peex113Ctg05333g00001.1). None of the additional AN4-like genes have measurable expression in stage 4 limbs; gene IDs and accessions given in Supplemental Table S4. Branch length represents substitutions per site. B, Normalized read counts of candidate MYB from three RNAseq experiments; L2FC values provided in Supplemental Table S2, data provided in Supplemental Dataset S5. In each RNAseq experiment, a single-colored Petunia species (P. exserta, P. secreta, P. inflata) was compared to the white P. axillaris as control, with three biological replicates each, depicted as dots; see Supplemental Figure S8, error bars ± 2sd. C, Of the candidate MYBs, ASE in P. axillaris and P. exserta could be measured on SNPs in DPL, PHZ, and PH4 only; highly significant ASE was detected in PHZ and PH4. ASE bars depict the mean read counts over all SNPs per gene, per species; points depict mean read counts per individual biological replicate. Number of SNPs indicated in parentheses, ASE calculated as major allele frequency indicated above bars, error bars ± se. *P < 0.05, **P < 0.01, ***P < 0.001

Anthocyanin biosynthesis is activated by the MBW complex, comprising MYB, bHLH, and WD40 transcription factors. MYB family SG6 members specialize in anthocyanin biosynthesis and control spatial and temporal specificity (Dubos et al., 2010; Feller et al., 2011) and are easily identified by their specific amino acid signatures (Stracke et al., 2001; Zimmermann et al., 2004; Lin-Wang et al., 2010; Hichri et al., 2011). The long-tube clade of Petunia has four known paralogs belonging to this SG that target different aspects of anthocyanin patterning: AN2 (petal limb), AN4 (petal tube, present in duplicate copies), DPL (petal vein), and PHZ (light-induced petal blush and vegetative tissues), whose functions have been demonstrated in P. hybrida (Quattrocchio et al., 1993; Quattrocchio et al., 1999; Albert et al., 2014).

Given the main MYB transcription factor that activates anthocyanin biosynthesis in petal limbs, AN2, is a pseudogene in P. exserta (Esfeld et al., 2018), an obvious hypothesis is that a closely related MYB substitutes for AN2 in the MBW complex. This led us to focus first on MYB SG6, to which AN2 belongs. To define the MYB SGs, we identified 256 and 246 MYBs from the P. exserta and P. axillaris genomes, respectively, and constructed a maximum likelihood phylogenetic tree with annotated Arabidopsis thaliana MYBs from Dubos et al. (2010). MYB SG6 formed a well-supported clade with A. thaliana PAP1, PAP2, MYB113, and MYB114. SG6 contained additional P. exserta members (Supplemental Figure S7A). To determine whether the presence of additional SG6 MYBs is unique to P. exserta or any species with anthocyanin-based color, we searched for SG6 MYBs in sister species P. secreta (out of a total of 177 MYBs) and analyzed a SG6-specific phylogeny (Figure 3B;Supplemental Table S5). Petunia axillaris, P. secreta, and P. exserta each have single copies of AN2, PHZ, DPL, and two copies of AN4 (AN4-1, AN4-2). Six additional AN4-like sequences in the P. exserta genome are pseudogenes (Figure 3B). This absence of novel or functional duplicated SG6 MYBs in P. exserta suggests a potential shift in expression or tissue specificity of a current SG6 MYB.

To identify de novo upregulation of SG6 MYBs in petal limbs, we compared P. exserta to the white, anthocyaninless P. axillaris. We subsequently compared P. exserta expression patterns to purple species P. secreta and P. inflata to characterize whether any MYB expression in P. exserta is unique, or whether patterns simply reflect petal limbs with anthocyanin biosynthesis (Figure 3C;Supplemental Figure S8A). AN2 expression was elevated in P. exserta, but the encoded protein is nonfunctional in both P. exserta and P. axillaris (Esfeld et al. 2018). PHZ has no obvious changes in protein function, and was equally and weakly expressed in P. axillaris, P. exserta, and P. secreta, and even more reduced in P. inflata. The AN4 genes were not expressed in petal limbs (Figure 3C, Supplemental Figure S8A). In contrast, DPL encodes a functional protein and was more highly expressed in P. exserta than in P. axillaris (Figure 3C;Supplemental Tables S2 and Table S6). Furthermore, both P. secreta and P. inflata (which have functional AN2) expressed DPL less than P. axillaris in petal limbs. DPL has weak ASE (Figure 3D;Supplemental Figure S8C), suggesting that it is mostly regulated in trans. DPL is located on chromosome 7, which has a minor QTL for anthocyanin content (Figure 3A;Supplemental Table S4). We examined the DPL protein sequence and found only a single amino acid replacement unique to P. exserta, which is not in a functionally annotated area of the protein (Supplemental Figure S9). Taken together, DPL is a paralog of AN2 that is upregulated in P. exserta petal limbs and therefore a promising AN2-replacement.

Broadening the search for potential transcription factors

Given that SG6 MYB proteins operate in a complex with additional protein partners (bHLH, WD40, other MYBs), we then expanded our search to identify any additional candidate genes from an entire transcription factor dataset (2,284 proteins from Yarahmadov et al. (2020)). We used the P. axillaris × P. exserta RNAseq dataset, filtering for statistically significant DE (1,388 genes), a base mean of ≥25 read counts (1,130 genes), and a log2 fold change (L2FC) of at least ±1.5 (153 genes). Along with ASE calculations, we assessed predicted functional variants (i.e. missense, nonsense, frameshift mutations, see Materials and methods).

Seven additional MYB transcription factors were differentially expressed, three of which were more highly expressed in P. exserta than in P. axillaris (Supplemental Table S6). Of these, two MYB transcription factors are known to influence or interact with floral color: PH4, which regulates vacuolar pH and influences floral color hue using the same bHLH and WD40 partners as AN2 (AN1/JAF13 bHLH and AN11 WD40) in SG “G20” (Quattrocchio et al., 2006b) (Figure 3C, Supplemental Figure S8A and Table S6) and MYBx, an R3 repressor of anthocyanin biosynthesis (Albert et al., 2014). The R2R3-MYB activator PH4 is more highly expressed in sister taxa P. exserta and P. secreta compared to P. axillaris, with the magnitude of difference greater in P. exserta (L2FC 1.85 P. exserta to P. axillaris, L2FC 0.89 P. secreta to P. axillaris; Supplemental Figure S8A; Supplemental Tables S2 and S6). PH4 additionally has moderate ASE (Figure 3D;Supplemental Figure S8C), suggesting both cis- and trans-regulation. The PH4 protein sequence contained one amino acid replacement which was shared between P. inflata and P. exserta (Supplemental Figure S10). The enhanced expression of PH4 in P. exserta qualifies PH4 as a valid non-SG6 candidate MYB.

Put in context with P. secreta and P. inflata expression, MYBx appears to be expressed similarly in P. exserta and P. secreta (L2FC 4.13 P. exserta to P. axillaris, L2FC 4.24 P. secreta to P. axillaris), but because it is a repressor that is more highly expressed in all of colored species rather than in white P. axillaris, it does not fit the candidate gene profile. Additionally, none of the MYBs expressed more highly in P. axillaris appear to be repressors.

Six bHLH and six WD40 transcription factors were differentially expressed. None of the WD40 functional annotations suggested an obvious interaction or role in floral color (Supplemental Table S7). Of the bHLH genes, JAF13 was less expressed in P. exserta than in P. axillaris (L2FC of 2.50; Supplemental Figure S8B and Table S8). As JAF13 is redundant with its more highly expressed paralog AN1 in the MBW complex (Quattrocchio et al., 1998; Spelt et al., 2000; Tornielli et al., 2009; Montefiori et al., 2015), it was not further considered. Of the remaining transcription factors, we considered only those transcription factors that have DE as well as ASE of ≥0.75 in order to locate genes with causal (in cis) mutations (Supplemental Table S9). Eleven genes were more highly expressed in P. exserta than P. axillaris, and ten genes more highly expressed in P. axillaris than P. exserta. Based on the detailed examination of the entire transcription factor dataset, we selected MYB-FL, DPL, and PH4 for functional validation.

Transcription factor candidate gene validation

To assess the contributions of each candidate gene to the floral pigmentation phenotype, functional validation was performed using virus-induced gene silencing (VIGS). MYB-FL silencing produced purple, anthocyanin-containing sectors in P. axillaris, demonstrating a negative association between flavonol and anthocyanin concentrations in the naturally occurring species, although it is not clear whether this is due to regulatory competition between the branches or dihydroflavonol substrate competition. DPL silencing produced white sectors; importantly, this phenotype was observed only in P. exserta and not in the other species tested (Figure 4A). Analysis of flavonoid concentrations in P. exserta showed a significant reduction of total anthocyanidins in DPL-silenced petal limbs (Figure 4B). Therefore, inactivation of DPL specifically interfered with anthocyanin production in P. exserta and had no effect in the purple species P. secreta and P. inflata.

Figure 4

VIGS phenotypes and chemistry in Petunia. A, Three candidate transcription factors and four candidate flavonoid biosynthetic pathway genes were silenced in the Petunia species using VIGS. Flowers where the visible phenotype was altered after gene silencing are indicated with a star. VIGS of DPL yielded a loss of color in P. exserta only, whereas VIGS of PH4 induced a loss of color in both P. exserta and P. secreta and a color shift in P. inflata. A color shift rather than a complete color loss (i.e. white sectors) was observed for the red P. exserta when either HT1 or HF genes were silenced. In the purple species, HT1 silencing did not yield a phenotype whereas a color loss (i.e. white sectors) were observed for HF silencing. VIGS of AAT yielded pink sectors in purple P. secreta and P. inflata. B, Total anthocyanidins (compounds shown as stacked bars in order from top to bottom: cyanidin, peonidin, delphinidin, petunidin, malvidin) were lower in P. exserta whole petal limbs (n = 5) where DPL was silenced compared to the control (pTRV2, n = 4), (Kruskal–Wallis X2 = 6.00, df = 1, P = 0.029). Individual values for total anthocyanidins shown as points. C, Total anthocyanidins were lower in P. exserta and P. secreta whole petal limbs where PH4 was silenced (n = 3) compared to the control (pTRV2; same control as in DPL-VIGS), but were not significantly different in silenced versus control P. inflata (n = 7, 2 respectively; P. exserta Kruskal–Wallis X2 = 4.50, df = 1, P = 0.034, P. secreta Kruskal-Wallis X2 = 3.857, df = 1, P = 0.049, P. inflata Kruskal-Wallis X2 = 1.191, df = 1, P = 0.275); individual values shown as points. D, Although silencing of HT1 in P. exserta petal limbs produces a slight decrease in cyanidin and delphinidin, a more drastic reduction of delphinidin is observed in HF silencing. Individual values for total anthocyanidins shown as points. E, For silencing of AAT in P. secreta, pooled petal limb sectors instead of entire petal limbs from three individuals per treatment were analyzed. Pink sectors contained lower concentrations of anthocyanidins than in purple sectors. Qualitatively, pink sectors also produced more delphinidin and less malvidin than did purple sectors, as well as small but detectable amounts of cyanidin in P. secreta. *P < 0.05, **P < 0.01, ***P < 0.001, error bars ± sd.

Silencing of PH4 produced a shift from purple to blue in P. inflata (Figure 4A). This same phenotype was observed in P. hybrida ph4 mutants (Quattrocchio et al., 2006b). Extensive studies of P. hybrida have shown that PH4 is an activator of vacuolar P-ATPase genes and that ph4 mutants have a reduced uptake of protons into the vacuole. When vacuoles become more basic, the same anthocyanins that appear red in acidic conditions appear more blue (Brouillard, 1988; Yoshida et al., 1995, 2003; Quattrocchio et al., 2006b). To our surprise, we observed white sectors in P. exserta as well as in sister species P. secreta. Indeed, PH4-silenced petals in P. exserta and P. secreta produced significantly lower levels of anthocyanidins, whereas PH4 silencing in P. inflata did not significantly lower the amount of anthocyanidins (Figure 4C).

We conclude that both DPL and PH4 are key transcription factors that are required for pigmentation in P. exserta. However, DPL inactivation exclusively affects P. exserta, whereas PH4 has a function in all three colored species.

Anthocyanin biosynthetic candidate gene validation

Biosynthetic genes F3′H/HT1, F3′5′H/HF1, F3′5′H/HF2, 3GT, and AAT were identified as promising candidates involved in the establishment of the red color. Silencing of the dihydroxylating F3′H/HT1 produced a visible phenotype in P. exserta (Figure 4A) with lighter-colored sectors. This subtle but highly reproducible phenotype could not be confirmed by chemical analysis of whole limbs, emphasizing the sensitivity of sector analysis (Figure 4D), and raising the possibility of redundant or compensatory 3′ B-ring hydroxylation by F3′5′H/HF. No change in phenotype was observed in P. secreta and P. inflata, which was to be expected as these species produce trihydroxylated purple pigments and do not express HT.

Next, we silenced the F3′5′H/HF trihydroxylation genes. Although HF1 and HF2 are not identical, the high sequence similarity prevented specific targeting of either one via VIGS-silencing (which uses siRNA molecules of 20–27 bp); thus, both HF1 and HF2 gene copies were simultaneously targeted. HF-silenced sectors in P. exserta were a light shade of red. When comparing the concentrations of anthocyanidins, the amount of delphinidin was reduced in the silenced sectors (Figure 4D). Further, the amount of cyanidin in F3′5′H/HF-silenced petals was significantly higher than in F3′H/HT1-silenced petals, demonstrating a shift in anthocyanidin hydroxylation. Thus, both F3′H/HT1 and F3′5′H/HF1/2 (and thereby both dihydroxylated and trihydroxylated anthocyanins) contribute to the red color of P. exserta. In contrast, F3′5′H/HF silenced sectors in P. secreta and P. inflata were completely white, marking the absence of anthocyanidins. This result indicates that the F3′5′H/HF genes are responsible for the purple pigmentation of these species whereas F3′H/HT1 activity is inconsequential.

The anthocyanin-modifying genes 3GT and AAT yielded contrasting results. Silencing 3GT produced no visible phenotype in any species, suggesting that it is not essential to either red or purple anthocyanin production (Figure 4A). Silencing AAT had no effect on the color phenotype in P. exserta, in line with its low expression (Figure 2; Supplemental Figures S5D and S6). In contrast, silencing AAT in the purple species P. secreta and P. inflata produced pink sectors (Figure 4A). The proportion of delphinidin increased in the pink sectors of both P. secreta and P. inflata, and additionally introduced dihydroxylated cyanidin in P. secreta pink sectors (Figure 4E). These results indicate not only that AAT is a key step in acylation of anthocyanins, but also that AAT-modified anthocyanins may act as substrate for methylated anthocyanins (i.e. petunidin, malvidin, peonidin) (Tornielli et al., 2009). Thus, both acylation and methylation strategies contribute to dark purple hues in Petunia. We conclude that P. exserta modifies expression of the hydroxylating enzymes F3′H/HT1 and F3′5′H/HF to produce its red hue, and that the low expression of downstream decorating enzyme AAT prevents a purple hue.

DPL, but not PH4, restores anthocyanin biosynthesis in white P. axillaris

To determine whether DPL can activate anthocyanin biosynthesis, transgenic lines expressing DPL under the control of the CaMV35S promoter were generated in the P. axillaris background. Of the nine independent lines, two representative lines were chosen for further analysis. Transgenic and nontransgenic siblings were compared from each line. The 35S plants were intensely pigmented in all visible parts of the plant including vegetative tissue (Figure 5A, Supplemental Figure S11). Petunia axillaris is white-flowered (genotype an2 an4 MYB-FL; Figure 5A) and 35S expression complements the an2 and an4 mutations, restoring anthocyanin biosynthesis to the corolla (petal limb and petal tube; Figure 5A;Supplemental Figure S11). Total anthocyanidin concentration increased in P. axillaris 35S petal limbs (Figure 5B).

Figure 5

Stable overexpression of 35S in P. axillaris. A, Overexpression of DPL yields purple flowers instead of white flowers in P. axillaris. B, The purple color is due to an increase in delphinidin anthocyanins with a small increase in cyanidin. The flavonol composition only changes slightly, with an increase in the monohydroxylated kaempferol. Petunia axillaris n = 3, 35S = 6 (three each from two independent lines). C, Gene expression analysis of transcription factors and biosynthetic genes for flavonoid biosynthesis in 35S transgenic lines (abbreviated as OE) compared to transgene negative P. axillaris (abbreviated as AX). Expression is represented as the mean ratio of the gene relative to expression of the reference gene SAND. Statistics were calculated using two sample or Welch t tests, **P < 0.01, ***P < 0.001. The full statistical analysis is given in Supplemental Table S11. D, VIGS of PH4 of in 35S lines induced white sectors compared to the control (pTRV2). E, Transient overexpression of 35S and 35S in wild-type P. axillaris, shown together with buffer and wounding control. Box plots depict medians of relative expression from three biological replicates, error bars ± sd. **P < 0.01, ***P < 0.001.

As expected, DPL transcript levels were highly increased in the transgenic line. Anthocyanin-specific biosynthetic genes were upregulated (DFR, ANS, 3GT, ART, AAT, 5GT, 3′AMT/MT, but not MF; Figure 5C). Expression of genes dedicated to flavonol biosynthesis, namely, MYB-FL, FLS and F3′H/HT1, were not significantly different in the overexpression lines compared to P. axillaris. Both F3′5′H/HF1, and F3′5′H/HF2 remained lowly expressed (Figure 5C). Thus, similar to P. exserta (Figure 2A), DPL overexpression in a P. axillaris background upregulated the genes of the anthocyanin pathway. However, DPL overexpression with the CaMV35S promoter appeared to upregulate all of the anthocyanin-modifying genes whereas this did not occur in P. exserta (low expression of 3GT and AAT; Figure 2A). We conclude that DPL can activate anthocyanin biosynthesis, but additional genetic variation in P. exserta as well as further degeneration in P. axillaris anthocyanin transcriptional network likely contributed to the red versus purple color.

We next investigated the potential role of PH4 in the transcriptional activation of the biosynthetic genes. Inactivation of PH4 by VIGS in the P. axillaris 35S transgenic background induced white sectors, as it did in the wild-type P. exserta and P. secreta backgrounds (Figure 5D, 4A). Thus, PH4 is indispensable for color formation even in a background where DPL is constitutively active. PH4 has been shown not to induce biosynthetic gene expression in P. hybrida (Quattrocchio et al., 2006b), but it could be an activator of DPL expression in P. exserta. However, transient expression of 35S in the wild-type P. axillaris background did not induce anthocyanin pigmentation (Figure 5E), which strongly argues against PH4 as an upstream regulator of DPL.

Discussion

In a scenario of shifts in pollinator competition and availability, selection for floral color is likely to be strong, favoring genetic changes of large phenotypic effect. Indeed, there is a rich and diverse literature on such large-effect genes. However, almost all documented cases of floral color changes involve loss-of-function mutations or severely reduced expression of the identified genes. In Petunia, we identified two cases of reacquisition of color from a colorless ancestor. Reacquisition of the ancestral purple color in P. secreta was by a simple 2-bp compensatory deletion in the R2R3-MYB transcription factor AN2 (Esfeld et al., 2018). Acquisition of red floral color in P. exserta turned out to be far from simple.

Balanced shift in activities of three hydroxylation enzymes

Documented shifts to red floral color are caused either by addition of carotenoids or by redirecting anthocyanin synthesis to the orange-red pelargonidins. Our new LC–MS analyses clearly show that P. exserta contains neither carotenoids nor pelargonidins (Figure 1; Supplemental Figure S1), but instead cyanidin and delphinidin are the major pigments. The presence of delphinidins is surprising as they tend to be blue/purple (Holton and Cornish, 1995). How then can P. exserta flowers be red?

The high abundance of dihydroxylated cyanidin and trihydroxylated delphinidin in P. exserta implies balanced activities of F3′H/HT1 and the two F3′5′ hydroxylases, F3′5′H/HF1 and F3′5′H/HF2. Compared with its sister species P. secreta, P. exserta displays a substantially increased expression of F3′H/HT1, in line with high cyanidin concentrations. In the case of the two F3′5′ hydroxylases, expression of F3′5′H/HF1 is decreased and expression of F3′5′H/HF2 is increased. Detailed analysis of the genetic interactions of these three genes in P. hybrida has shown that HF1 is fully epistatic over both HT1 and HF2; that is, in the presence of HF1 activity, anthocyanins are always trihydroxylated irrespective of the allelic state of the two other genes (Wiering and De Vlaming, 1984). Genotypes without HF1 but with HT1 and HF2 produced both dihydroxylated and trihydroxylated anthocyanidins. Therefore, we propose that reduction of F3′5′H/HF1 activity was a prerequisite for the accumulation of dihydroxylated cyanidins by F3′H/HT1. In addition, the Petunia-specific duplication of F3′5′H and subsequent functional divergence of HF1 and HF2 specificity was also vital to the delphinidin-based contribution to the red color in P. exserta.

Anthocyanin acylation and methylation are further contributors to the purple-to-red continuum

Shifts in floral color from blue to red across the angiosperms always involve a decrease in anthocyanin hydroxylation. But to the best of our knowledge, P. exserta is the only red-flowering species in the Solanaceae and possibly in the eudicots that retains delphinidin production and still produces a red hue (Berardi et al., 2016; Ng and Smith, 2016a, 2016b; McCarthy et al., 2017; Larter et al., 2018; Ng et al., 2018). To resolve this conundrum, we focused our attention on the sugar, acyl or methyl modifications of the anthocyanidin backbone. Without acylation, the simply glycosylated versions of cyanidin and delphinidin appear redder in situ than do their acylated versions (Curaba et al., 2019; Tasaki et al., 2019). In support of this, when anthocyanin acyltransferase AAT was silenced in purple P. secreta or P. inflata, the silenced sectors appeared pink (Figure 4).

We propose that in addition to essential changes in anthocyanin hydroxylation, P. exserta is red because it lacks decoration of the anthocyanidin backbone by acylation (AAT), and as a consequence color hue shifts toward red (Fukui et al., 1998; Slimestad et al., 1999; Hashimoto et al., 2002). Indeed, there may be many undetected cases where important anthocyanin modifying genes, such as AAT contribute to shifts in floral color, and considering the full complexity of the anthocyanin biosynthetic pathway could be rewarding.

Regain of color from a colorless ancestor

Unlike sister species P. secreta, P. exserta did not resurrect the AN2 pseudogene, raising the question of which transcription factor(s) induce the anthocyanin biosynthetic genes in floral tissue. Genetic analysis identified a strong QTL on chromosome 2, with opposing effects on visible and UV color (Figure 3A). We previously identified transcription factor MYB-FL as the gene underlying the UV color QTL (Sheehan et al., 2016). Transposon insertions in MYB-FL yielded UV-reflective pink sectors in the white UV-absorbing background, most likely due to substrate competition between flavonol and anthocyanin biosynthesis. It is highly plausible that MYB-FL underlies the visible color QTL on chromosome 2, as well. UV-absorbing flowers stand out against the foliage background and are an important guide for night-active hawkmoths, whereas UV color has little effect on day-active pollinators (Chittka et al., 1994; White et al., 1994; Raguso and Willis, 2005). Thus, MYB-FL activity is dispensable for daytime pollination by hummingbirds, and its loss is necessary to free up substrate for anthocyanin biosynthesis. It is a rare example of gain of function by a loss-of-function mutation.

DPL reinstates anthocyanin biosynthesis

AN2 being nonfunctional, we searched for substitute transcriptional activators that intensify visible flower color by activating the anthocyanin biosynthetic genes (Figure 2B). Within the R2R3-MYB SG6 clade to which AN2 belongs, we found a single candidate, DPL. DPL is moderately upregulated relative to P. axillaris, whereas it is low in purple P. secreta and P. inflata (Figure 3C). Silencing of DPL strongly reduced color in P. exserta but not in P. secreta and P. inflata (Figure 4, A and B), providing functional validation of its unique role in P. exserta. We further demonstrated that DPL is capable of reinstating anthocyanin biosynthesis in the anthocyaninless species P. axillaris (Figure 5).

In P. hybrida Mitchell (which is mostly P. axillaris-like) DPL is responsible for vein pigmentation in the petal tube but has no role in the limb (Albert et al., 2011). We conclude that DPL has shifted from activating anthocyanin biosynthesis in the veins of the floral tube to activating anthocyanin biosynthesis in the limb. Duplication and diversification of MYBs is pervasive in floral color diversification (Des Marais and Rausher, 2008; Yuan et al., 2014; D'Amelia et al., 2018). The duplication of the ancestral SG6 MYB presumably took place in the Petunia ancestor given the presence in all Petunia lineages and phylogenetic relatedness of the four anthocyanin-MYB paralogs (AN2, DPL, PHZ, AN4; Figure 3B). The upregulation of DPL in the P. exserta lineage is an essential step in the regaining of intense anthocyanin pigmentation in its flowers, and is an example of the gain of a new expression pattern (regulatory neofunctionalization) rather than partitioning of ancestral expression patterns (subfunctionalization) (Moore and Purugganan, 2005).

A novel role of PH4 in the Petunia long-tube clade

PH4 is 3.6-fold upregulated in P. exserta compared to P. axillaris and has moderate ASE (0.64), suggesting it may be a candidate AN2 substitute (Figure 3C). Silencing the gene in the short-tube species P. inflata yielded the expected shift to blue, presumably through disruption of vacuolar pH acidification. In contrast, PH4 silencing results in white sectors in P. exserta as well as P. secreta, indicating a new function specific to the long-tube clade (Figure 4, A and C). Could this new function be transcriptional activation of anthocyanin biosynthesis, either directly or indirectly? PH4 does not induce anthocyanin biosynthetic genes in P. hybrida (Quattrocchio et al., 2006b) and neither does it induce anthocyanin biosynthesis in P. axillaris (Figure 5E). Moreover, silencing of PH4 causes white sectors even when DPL is constitutively expressed (Figure 5D). Therefore, PH4 does affect pigment accumulation, but not by direct transcriptional activation of the biosynthetic genes nor indirectly by activating DPL. We conclude that PH4 induces anthocyanin accumulation by an independent mechanism. What could this be? In P. hybrida, PH4 is required for acidification of the vacuole as well as for volatile transport out of the vacuole (Quattrocchio et al., 2006b; Cna'ani et al., 2015). We speculate that PH4 activates a transporter that exports anthocyanins from the cytoplasm into the vacuole.

A complex molecular mechanism in the shift to hummingbird pollination

Most transitions in floral color studied at the molecular level appear to be relatively simple. The most complex transition described to date is the blue to red shift in Iochroma (Solanaceae), which involved modifications of DFR, F3′H, and F3′5′H (Smith and Rausher, 2011; Larter et al., 2018). That P. exserta retains delphinidin production, yet presents a red color, qualifies it as more complicated than a transition that goes from delphinidin (blue–purple) to pelargonidin (brick red). In the case of P. exserta, our major conclusions are that the evolution of red color in P. exserta involved (1) inactivation of a competing biosynthetic pathway by loss of function of its specific transcriptional activator MYB-FL; (2) upregulation of DPL to replace the ancestral function of AN2 in petal limbs; (3) reshuffling of the expression patterns of the hydroxylating genes F3′H/HT1, F3′5′H/HF1 and F3′5′H/HF2; and (4) absence of anthocyanin acylation due to low expression of the responsible enzyme, AAT (Figure 6).

Figure 6

Schematic for the mechanism of red flower color in P. exserta. The evolution of red flower color in P. exserta involved the loss of function in the flavonol biosynthetic transcription factor MYB-FL and a loss of function in the ancestral anthocyanin transcription factor AN2. To compensate for the loss of AN2, related MYB transcription factor DPL is upregulated in petal limbs to activate anthocyanin biosynthesis. Expression of the transcription factor PH4 is essential for anthocyanin pigmentation in P. exserta, although its function is not directly related to biosynthesis but potentially anthocyanin transport. Together with re-activation of anthocyanin biosynthesis, other essential changes in gene expression are necessary to achieve red pigmentation relative to sister taxa: moderate expression of F3′H/HT1 (in bold), upregulation of F3′5′H/HF2 (in bold) instead of F3′5′H/HF1. Last, the downregulation of AAT expression prevents cyanidin and delphinidin anthocyanin pigments from acylation which would impart a blue–purple hue. Transcription factors are indicated in ovals, flavonoid structural genes/enzymes without a bounding shape, and biochemical products in boxes. Items in gray indicate inactivity or low expression (in the case of genes) or low concentration (in the case of products).

Finally, there must be additional levels of complexity. First, overexpression of DPL behind the strong 35S promoter in P. axillaris induces high expression of all late biosynthetic genes, including 3GT and AAT which are not expressed in P. exserta (Figures 5 and 2). This hints at additional genetic diversity between the P. axillaris and P. exserta backgrounds. It is also possible that the P. exserta DPL expression level (Figure 3C) is insufficient to induce certain target genes. In that case, a precise intermediate level of DPL expression would be functionally relevant. Second, ASE analysis shows that DPL is largely activated in trans (Figure 3D). A plausible candidate upstream activator was PH4, but we have confidently ruled it out (Figure 5, D and E). Broadening the search to the entire transcription factor data set yielded 21 potential candidates. That none of them is annotated as even remotely related to flavonoids does not necessarily disqualify them as candidates. Indeed, it has been well established—also in plants—that many different transcription factors bind to many different promoter elements (Franco-Zorrilla et al., 2014; Taylor-Teeples et al., 2015; Gaudinier et al., 2018), potentially creating large and complex transcriptional networks. Mutational effects in such networks result in a continuum of large to small effects in a given phenotype and may diffuse through the network causing subtle changes in unrelated pathways (Boyle et al., 2017), in line with theoretical models (Orr, 1998, 2005). To what extent such a scenario is relevant for the evolution of red color in P. exserta case remains an interesting question.

Materials and methods

Plant material and growth conditions

Wild Petunia accessions have been previously described (Segatto et al., 2014; Turchetto et al., 2014, 2016; Sheehan et al., 2016). The reference accession P. axillaris N is from the Rostock Botanical Garden (Germany) and is registered in the Amsterdam collection under the designation of P. axillaris S26; P. inflata S6 (hereafter referred to as P. inflata) was provided by R. Koes (University of Amsterdam, the Netherlands); P. exserta from R.J. Griesbach (USDA, Beltsville, USA). Petunia secreta was collected in its natural habitat and maintained through laboratory crossings (Esfeld et al., 2018). Plants were grown in a growth chamber under a light:dark regime of 15:9 h using 400-W Clean Ace Daylight metal halide lamps (Eye Lighting International; two bulbs per square meter giving 200–250 µmoL m−2 s−1), at 22°C:17°C at 60%–80% relative humidity, in commercial soil (70% Klasman substrate, 15% Seramis clay granules, 15% quartz sand), and fertilized once a week.

Color and UV images

Color images were recorded using a Panasonic DMC-TZ10 camera. UV pictures were recorded using a Nikon 60-mm 2.8-D microlens with a Nikon D7000 SLR camera converted to record UV light using a UV-specific filter (blocking visible and infrared light).

DNA sequencing and reference genomes

Reference draft genomes of P. axillaris and P. exserta used were taken from DNAZoo (https://www.dnazoo.org/assemblies/Petunia_axillaris and https://www.dnazoo.org/assemblies/Petunia_exserta) and modified, detailed in Supplemental Methods S1, to reflect correct chromosome arrangement and names. The reference genome for P. inflata was v1.0.1 from (Bombarely et al., 2016), assembly and annotation files accessible on the SolGenomics website (https://solgenomics.net/organism/Petunia_inflata/genome). We additionally completed a draft genome assembly and annotation of P. secreta for this article. This Whole Genome Shotgun project has been deposited at NCBI GenBank under the accession JAFBXY000000000, project PRJNA674325. The version described in this article is version JAFBXY010000000, and annotation is deposited at Dryad (https://doi.org/10.5061/dryad.jsxksn083). Further methods on DNA sequencing and assembly are detailed in Supplemental Methods S1.

RNA sequencing

To detect differential gene expression and ASE in P. exserta, P. axillaris, the P. axillaris × P. exserta F1, and P. secreta, petal limb tissue was harvested via dissection of petal limbs from stage 4 buds (P. axillaris/P. secreta: 22–30 mm, P. exserta/P. axillaris × P. exserta F1, 25–34 mm) from plants grown under controlled conditions. For each species or hybrid, petal limbs from three biological replicates (individual plants) were collected. RNA was extracted using the Qiagen RNeasy Plant Mini Kit. RNA was prepared and sequenced in the Lausanne Genomic Technologies Facility (Lausanne, Switzerland) in two separate experiments. Experiment 1: P. axillaris, P. exserta, P. axillaris × P. exserta F1; Experiment 2: P. axillaris, P. secreta. Quality of RNA was checked using a Fragment Analyzer (Advanced Analytical). Libraries were prepared using Illumina TruSeq PE Cluster Kit v3 and each experiment was sequenced on a single lane of an Illumina HiSeq 2500 as single-end 100-nt reads. Sequencing data were processed using the Illumina Pipeline Software v.1.82. Reads were uploaded to NCBI SRA database under PRJNA674380.

Transcriptomes from Illumina data were assembled for P. exserta, P. axillaris, and P. secreta using Trinity v2.4.0 in genome-guided mode with P. axillaris v3.0.4 as the reference genome (Grabherr et al., 2011; Haas et al., 2013) for the purposes of examining predicted transcripts.

DE analysis

Although we collected material for RNAseq data from Petunia species from the same tissue type at the same developmental stage (petal limbs from bud stage 4) and under controlled growth conditions, samples were not sequenced in the same experiment. Thus, to avoid any experimental or platform-based biases, we did not analyze all datasets together. Each colored species (P. exserta, P. secreta, P. inflata) was always sequenced with the white P. axillaris. We carried out DE analysis for each of the three experiments (P. exserta, P. secreta, P. inflata compared with P. axillaris) in parallel, and statistical analysis was only carried out within experiments (colored species compared with P. axillaris). The P. inflata/P. axillaris experimental files are NCBI SRA accession PRJNA300613.

All read data were processed with Trimmomatic v.0.36 (Bolger et al., 2014) to remove Illumina adaptor sequences and trim low-quality regions. These preprocessed reads were mapped against the draft reference genome of P. axillaris (v.3.0.4, described above) using STAR v.2.6.0c in two-pass mode, with splice junctions—sjdbOverhang 100 and ignoring reads that map more than 20 times in total (Dobin et al., 2013). Reads were counted per gene using featureCounts v.1.5.2 (Liao et al., 2014).

DE analysis was performed with DESeq2 v.1.26.0 (Love et al., 2014) in R v.3.6.0 (R Core Team, 2019) using RStudio v.1.3.1073. Counts were normalized using rlog-transform in DESeq2 and for each gene mean counts were computed over the sample replicates. For each comparison, the colored species (P. exserta, P. secreta, P. inflata) was compared with the P. axillaris from the same experiment.

ASE in the P. axillaris × P. exserta F1 hybrids

We performed ASE analysis as in Esfeld et al. (2018) and Yarahmadov et al. (2020). The bamfiles of three P. axillaris × P. exserta F1 hybrids were used to detect variants according to GATK Best Practices for RNA-seq data (Van der Auwera et al., 2013). In short, after duplicate marking with Picard-tools v2.18.11 (http://broadinstitute.github.io/picard/) and splitting reads with N in their CIGAR string, local realignment around indels was undertaken and base quality scores were recalibrated, using a set of high-quality SNPs determined by an initial run of the HaplotypeCaller (DePristo et al., 2011). Variants were filtered using hard thresholds, selecting for biallelic alleles and removing variants matching DP < 10, AF < 0.75, QD < 2.0, MQ < 40.0, FS > 60.0. Clustered SNPs with more than three occurrences in a window of five were also removed.

Allelic coverage for variant positions was detected with ASEReadCounter implemented in GATK v3.5.0 (Van der Auwera et al., 2013; Castel et al., 2015) using -mmq 50 and -minDepth 20 filters on rounded read counts averaged over the biological replicates. Analyses of allelic imbalance were conducted in R (v3.6.0) with the package MBASED (v.1.20.0) (Mayba et al., 2014). Parameters of read count over-dispersion were estimated with a custom R script provided by the author of the MBASED package. P-values of ASE and heterogeneity were corrected for multiple comparisons using the p.adjust function with the Benjamini–Hochberg method in R base package “stats” (v3.6.0). ASE scripts are made available on https://github.com/Kuhlemeier-lab/Exserta-red/.

Analysis of functionally relevant SNPs in the coding region of candidate transcription factors

To identify functional SNPs between P. exserta and P. axillaris, variant files of the parental species obtained from ASE analysis described above were scanned using SNPeff v.4.3T (Cingolani et al., 2012). Detected mutations classified as high-impact, frameshift, loss and gain of stop codon, and missense variants were considered for further analysis.

RIL population, GBS, and genetic map

An F7 recombinant inbred line (RIL) mapping population of 195 progenies were bred from selfed progeny of an F2 population of the parental species P. axillaris × P. exserta. These individuals were sequenced and genotyped using a GBS (genotyping by sequencing) protocol (Elshire et al., 2011); details further described in Supplemental Methods S2. Raw sequence reads are available at the NCBI SRA database under PRJNA704924,

Genetic maps were constructed using packages R/qtl v1.46-2 (Broman et al., 2003) and ASMap v.1.0-4 (Taylor and Butler, 2017) in R v3.6.0 and RStudio v1.3.1073 (R Core Team, 2016). Individuals with <70% genotypes were removed. Markers in initial draft genetic maps with identical genotypes over all individuals were binned for a second round of map construction using a cutoff of P-value of 110. A genetic map of 714.1 cM was generated with Kosambi mapping function and was used to group 1,409 markers into seven linkage groups. The finalized genetic map is summarized in Supplemental Table S10. The genetic map generation script is made available on https://github.com/Kuhlemeier-lab/Exserta-red/.

Phenotyping and QTL analysis of P. axillaris × P. exserta F7 RIL population

Anthocyanin and flavonol absorbances were measured in the P. axillaris × P. exserta F7 RIL population using the same protocol as described by Sheehan et al. (2016). For each plant, five flowers were phenotyped. For each flower, a disk 8 mm in diameter of petal limb tissue was sampled and placed into 1 mL of MAW extraction buffer (2:1:7 methanol:acetic acid:water). After 48 h in the dark, absorption spectra were measured on a spectrophotometer (Ultraspec 3100 pro, Amersham Biosciences, England, UK). Anthocyanin values represent summed absorbance values over 445–595 nm and flavonol values represent summed absorbance values over 315–378 nm.

Anthocyanins and flavonols were evaluated for the identification of QTLs based on their phenotypic distribution using the “scanone” function in R/qtl v1.46-2 (Broman et al., 2003). Anthocyanins were normally distributed and analyzed with the “hk” Haley–Knott regression method. Flavonols showed a bimodal distribution reflecting groups with high and low concentrations and were analyzed using the “2part” method. Genome-wide significance level was established using 10,000 permutations using the “n.perm” argument for each of the anthocyanin and flavonol phenotypes at α = 0.05. The R/qtl script is made available on https://github.com/Kuhlemeier-lab/Exserta-red/.

Identification and filtering for candidate MYB, bHLH, WD40, and global transcription factors

MYB transcription factors were detected in the P. exserta, P. axillaris, and P. secreta genomes using HMMER v3.1b2 (Eddy, 2009); P. secreta MYBs were only used to identify SG6 MYBs. The alignment of A. thaliana MYBs from Dubos et al. (2010) was used to create a MYB hmm profile, and then searched against each genome. A cutoff of E-value of 0.01 for the full sequence was used. WD40 and bHLH transcription factors were identified in P. exserta and P. axillaris only. Identification of WD40 proteins was performed with HMMER as described for MYBs above, but the hmm profile was constructed from WD40 proteins from the Solanaceae downloaded from pfam database (https://pfam.xfam.org/family/WD40#tabview=tab7). The bHLH proteins were previously identified in Yarahmadov et al. (2020) for P. exserta and P. axillaris.

Once identified, we extracted the MYB, bHLH, and WD40 transcription factors present in the DE gene dataset from the RNAseq experiment of P. exserta, P. axillaris, and their F1 hybrid, and then applied strict filters (DE, filtering for ±1.5 LFC, and only considering genes that had a baseMean of least 25 normalized read counts).

We additionally scanned a Petunia global transcription factor dataset (2,284 proteins from Yarahmadov et al. (2020)) to extend the search outside of potential MBW candidates to any possible transcription factor (thus including other major transcription factor families, such as WRKY or Zinc fingers). We applied the same DE filters but additionally applied a filter for an ASE value of 0.75; this would detect mutations in cis. This approach served to narrow down the candidate gene list and would capture additional candidate transcription factors.

Phylogenetic analysis: MYB tree estimation

Two gene trees were estimated: one for all MYBs to define the MYB SGs, and a second tree focused on MYB SG6. For the entire-MYB tree, MYBs from P. exserta, P. axillaris, and A. thaliana were aligned using Clustal and MAFFT algorithms in Geneious v9.1.8 (Biomatters Ltd.). We included the recently discovered SG6 MYB ASR3 from P. inflata (Genbank accession MF623311) and used the C-myb protein from Danio rerio as an outgroup (Genbank accession AAH59803) as in Gates et al. (2018). A neighbor-joining tree was generated in Geneious using default parameters to serve as a starting tree for phylogenetic estimation. A Perl script for model selection available with RAxML estimated the best protein evolution model as PROTGAMMAVT. The phylogeny was estimated in RAxML v8.2.10 (Stamatakis, 2014) by estimating the best-scoring maximum likelihood tree with 1,000 bootstrap replicates (using the standard RAxML bootstrapping algorithm, -f b) with raxmlHPC. Petunia MYBs were classified into the canonical MYB SGs if they formed discrete clades with known A. thaliana proteins, as well as by using motifs defined in Stracke et al. (2001).

The second gene tree focused on MYB SG6 and the PH4 clade (“G20”), and included protein sequences from P. exserta, P. axillaris, and P. secreta. MYBs from SG6 were identified from the larger MYB tree by their well-known motifs: the bHLH interacting motif [D/E]Lx2[R/K]x3Lx6Lx3R, the ANDV motif, and the C-terminal SG6-defining motif [R/K]Px[P/A/R]x2[F/Y/L/R] (Stracke et al., 2001; Zimmermann et al., 2004; Lin-Wang et al., 2010; Hichri et al., 2011). The PH4 sequences were identified using the same bHLH interaction motif as well as two “G20” conserved motifs (Quattrocchio et al., 2006b). The full protein sequences of these genes were extracted, realigned using the MAFFT algorithm, and then a neighbor-joining tree was constructed in Geneious v9.1.8 (Biomatters Ltd.). The gene tree was estimated with raxmlHPC program using the rapid bootstrap approach (-f a) with 1,000 bootstraps, PROTGAMMAAUTO model selection (JTT was selected), and with the neighbor-joining starting tree.

Fasta alignments and newick tree files for both trees are provided as Supplemental Datasets S1–S4.

Coding sequences of flavonoid biosynthetic loci and transcription factors

We used known sequences from P. hybrida as well as those previously identified in Supplemental Note 7 from Bombarely et al. (2016) and in Esfeld et al. (2018), and used BLAST to find orthologs in P. axillaris, P. inflata, P. secreta, and P. exserta. Genes were aligned first with the MAFFT algorithm and then by codon in Geneious v9.1.8 (Biomatters Ltd). We manually examined the coding sequences of flavonoid biosynthetic loci HT1 (F3′H), HF1 (F3′5′H), HF2 (F3′5′H second copy), FLS, DFR, ANS, 3GT, ART, AAT, 5GT, MT (3′AMT), MF1 (3′5′AMT), and MF2 (3′5′AMT second copy) for any loss-of-function mutations and estimated sequence divergence using the “Distances” function in Geneious. Protein alignments are available on Dryad.

The DFR, HT1, HF2, and AAT coding sequences were translated to amino acid sequences and aligned using the MAFFT algorithm in Geneious v9.1.8 (Biomatters Ltd.). The DFR protein alignment of P. inflata, P. axillaris, P. secreta, P. exserta, and P. hybrida (Genbank KC140107.1) were re-aligned to the sequence and crystal structure of Vitis vinifera DFR (Petit et al., 2007). Protein active and binding sites are known for the hydroxylating genes F3′H and F3′5′H and AAT (Nakayama et al., 2003; Seitz et al., 2007). We examined these sites in amino acid alignments of HT1, HF1, HF2, and AAT for any amino acid changes unique to P. exserta.

Coding sequences for DPL and PH4 in P. exserta, P. axillaris, and P. secreta were cloned during the Gateway cloning procedure (described below) and verified by Sanger sequencing. Exons and introns were defined by aligning the resulting sequences to the respective genome drafts, as well as manually examining Illumina read mapping and assembled transcripts (transcriptomes described above).

Virus induced gene silencing

VIGS was performed as described in (Spitzer-Rimon et al., 2013) with the Tobacco rattle virus (TRV). Briefly, pTRV1, pTRV2-MCS (multiple cloning site), and pTRV2-NtPDS (phytoene desaturase from Nicotiana tabacum, which causes photobleaching) plasmids were obtained from the Arabidopsis Biological Resource Center (ABRC accessions CD3-1039, CD3-1040, CD3-1045, respectively). We selected coding region fragments in the targeted genes that are the most specific to avoid off-targets (∼200–350 bp). These fragments were either synthesized by Genewiz (Leipzig, Germany) or amplified by PCR from cDNA using forward and reverse primers containing BamHI and EcoRI restriction sites (summarized in Supplemental Tables S12 and S13). Fragments were cloned into the pTRV2-MCS plasmid and transformed into A. tumefaciens strain GV3101. pTRV1 was mixed with each pTRV2 derivative in a 1:1 ratio prior to infection of Petunia plants (at approximately a six-leaf stage). Infection was carried out by removing the shoot apex and applying 40 μL inoculum to the cut surface of the stem. Only flowers from the branches arising from infected meristems were phenotyped.

For each species/gene combination, three plants each were infiltrated with pTRV2 derivative, empty vector (pTRV2-MCS) as a negative control, and one pTRV2-NtPDS positive control (for VIGS efficiency). At the time of flowering, at least 25 consecutive flowers per plant (at least ∼75 total flowers per gene/species) were phenotyped 2 d postanthesis at the same time every day. Silenced sector phenotypes were determined when the same phenotype arose in all three biological replicates per species and produced in the majority (>50%) of flowers for P. axillaris, P. exserta, and P. inflata. Petunia secreta is more susceptible to TRV and tended to produce silenced sectors at a lower, but consistent rate (∼25%–30%).

Overexpression constructs and transformation in P. axillaris

The construct used for stable transformation, pNWA12 (35S) from Albert et al. (2011), was kindly gifted by Nick Albert. Stable transgenic P. axillaris lines were generated by leaf disk transformation with Agrobacterium tumefaciens strain LBA4404 following an adapted version of the protocol described by Conner et al. (2009). Nine independent lines from the T1 generation were established with two independent lines selected for characterization. Further methods are described in Supplemental Methods S3.

The P. exserta ORFs for DPL and PH4 were amplified and cloned into the Gateway-compatible binary vector pGWB402 (Addgene plasmid #74796) with a CaMV35S promoter. Transient transformation of petals using A. tumefaciens strain GV3101 was performed as in Van Moerkercke et al. (2011). One-day open flowers were syringe-infiltrated at each of the five petal midveins, and 3 d later phenotypes were observed and photographed.

Reverse transcription-quantitative PCR

Expression levels of biosynthetic genes were measured for stages 4 and 6 petal limb material in P. axillaris, P. exserta, and P. secreta. Primer sequences not previously reported in Esfeld et al. (2018) are included in Supplemental Table S12. Stage 4 expression levels of P. axillaris and P. secreta are published in Esfeld et al. (2018), but were conducted in the same experiment as presented in this manuscript, including bud stage 4 data from P. exserta as well as bud stage 6 data.

Petal limbs were collected from buds stages 4 and 6 from three different plants of each species, representing three biological replicates (n = 12); P. axillaris and P. secreta (bud length stage 4: 22–30 mm; bud length stage 6: 40–50 mm), P. exserta (bud length stage 4: 25–35 mm; bud length stage 6: 50–65 mm). Samples were frozen in liquid nitrogen and stored at −80°C until extraction. RNA was extracted with the RNeasy Plant Mini Kit (Qiagen), replacing β-mercaptoethanol by dithiothreitol. RNA was treated with DNase I (Sigma-Aldrich), and quantified with a Nanodrop ND-1000 (Thermo Fisher); RNA was additionally analyzed on a Bioanalyzer 2100 (Agilent Technologies) prior to cDNA preparation.

Random hexamer primers were used for first-strand synthesis using the Transcriptor Universal cDNA Master (Roche) kit according to the supplier. Quantitative RT-PCR experiments were performed using the LightCycler 96 Real-Time PCR System (Roche) with KAPA SYBR_FAST qPCR Kit optimized for LightCycler 480 (KAPA Biosystems), according to the manufacturer’s recommendations. Technical replicates were obtained by running each sample in triplicate, means of these were taken as the value for the biological replicate. Cycle of quantification (Cq) thresholds and normalization calculations were determined by the LightCycler 480 Software (v.1.1.0.1320; Roche). Standard curves were used to determine the PCR efficiencies. The SAND reference gene (Mallona et al., 2010) and controls with no-reverse transcription were included for each sample. For each of the flavonoid genes of interest, pairwise comparisons of the mean relative gene expression values were made among the three species and bud stages in R (v3.6.0) using the TukeyHSD function.

For RT-qPCR expression analyses in the stable transgenic lines (35S in P. axillaris), the experiment was conducted as described above with the following adaptations. Petal limbs were collected from buds at stage 4 only from three stable transformant plants and from three negative transformants (n = 6); bud length 22–30 mm. The experiment was conducted as described above with the following adaptations. First strand synthesis was performed with random hexamer primers using the qScriber cDNA Synthesis Kit (HighQu) according to the supplier. qPCR experiments were performed using the LightCycler 96 Real-Time PCR System (Roche) with ORA qPCR Green ROX L Mix, 2X (HighQu), according to the manufacturer’s recommendations but reducing the reaction volume by half. For each of the genes, comparisons of the mean relative gene expression values were made using two sample t tests (if heterogeneity of variances) or Welch two sample t tests (if heterogeneity of variances) in R (v.3.6.0).

Extraction and identification of flavonoid pigments

To quantitate and identify anthocyanidins (anthocyanin aglycones) and flavonol aglycones, acid hydrolysis procedure was followed as described in Esfeld et al. (2018) based on Harborne (1998). Briefly, fresh petal limb tissue was dissected from three flowers from three individual plants per species (or per VIGS treatment: either whole limbs or sectors were sampled, sectors when possible; on average, n = 9 for each treatment presented) was weighed and placed overnight at 4°C in 1 mL of 2 M HCl in a darkened 2-mL screw-cap tube. Flavonoid extraction was performed via acid hydrolysis for 90 min with phase separation using ethyl acetate and isoamyl alcohol, resulting in two fractions per sample: one containing flavonols and the other containing anthocyanidins. Final extracts were dried at 4°C and stored at −20°C until analysis. Samples were resuspended in 75% MeOH with 1% formic acid (v/v) and the two flavonoid fractions per sample were combined for analysis, except for P. axillaris samples which were analyzed separately due to the disparate concentrations of flavonols and anthocyanidins. Further analyses are described in Supplemental Methods S4.

To qualitatively identify intact anthocyanins, three second-day open flowers each of P. exserta, P. secreta, and P. inflata were sampled and combined per species. Care was taken to avoid or remove pollen from the fresh petals used. Fresh petals were ground on ice in 1 mL of methanol 1% HCl (v/v), diluted, and immediately analyzed using LC–MS. Further methods are described in Supplemental Methods S4.

Crude petal pH measurements

To measure the pH of crude petal homogenates, we followed the protocol described in Spelt et al. (2002). The petal limbs of two flowers (P. exserta, P. axillaris, P. secreta) or three flowers (P. inflata; approximately three P. inflata corollas equals the weight of two flowers for the other species) were ground with a pestle and mortar in 4 mL of distilled water, and the pH of the homogenate was measured immediately (to avoid alkalinization by uptake of atmospheric carbon dioxide) with a pH electrode that corrects for temperature (Hanna Instruments HI991001).

Accession numbers

Sequence data from this article can be found under the accession numbers listed in Supplemental Table S14. Additional Genbank/GenPept accession numbers for reference proteins used in this work are as follows: P. inflata ASR3 (MF623311.1), P. hybrida DFR (KC140107.1), Vitis vinifera DFR (P51110.1), P. hybrida Mitchell DPL (ADW94950.1), P. hybrida R27 PH4 (AAY51377.1), and Danio rerio C-MYB (AAH59803.1). The P. inflata/P. axillaris experimental files are NCBI SRA accession PRJNA300613. Raw sequence reads for the F7 mapping population of P. axillaris × P. exserta are available at the NCBI SRA database under PRJNA704924. The Whole Genome Shotgun project has been deposited at NCBI GenBank under the accession JAFBXY000000000, project PRJNA674325. The version described in this paper is version JAFBXY010000000, and annotation is deposited at Dryad (https://doi.org/10.5061/dryad.jsxksn083).

Supplemental data

Identification of four main anthocyanins in P. exserta methanolic extracts.

Identification of an acylated anthocyanin present in purple P. secreta and P. inflata methanolic extracts.

Determination of crude petal extract pH in four Petunia species.

Protein alignments of Petunia DFR, HT/HF1/HF2, and AAT against known crystal structures, binding sites, and motifs.

RNAseq expression of phenylpropanoid and flavonoid biosynthetic genes.

Expression of flavonoid biosynthetic pathway genes in petals from bud stages 4 and stage 6.

Maximum likelihood (ML) phylogenetic tree of MYBs and SG relationships.

RNAseq and ASE in the F1 hybrid between P. axillaris and P. exserta of known flavonoid transcription factors.

Protein alignment of DPL.

Protein alignment of PH4.

Overexpression of DPL under the control of the CaMV35S promoter in P. axillaris.

RT-qPCR statistics for Petunia biosynthetic gene expression.

L2FC values and statistics for three RNAseq experiments.

Summary of QTL analysis of anthocyanin and flavonol quantity in petal limbs.

Chromosomal locations of phenylpropanoid, flavonoid pathway structural genes and transcription factors based on the P. axillaris and P. exserta draft genomes.

Gene IDs for MYB SG6 and PH4 clade phylogeny.

Filtering for MYB candidate transcription factors.

Filtering for WD40 candidate transcription factors.

Filtering for bHLH candidate transcription factors.

Filtering for transcription factors with ASE.

Summary statistics of the genetic map for F7 RILs population of P. axillaris × P. exserta.

RT-qPCR statistics for 35Spro:DPL experiment biosynthetic gene expression.

Primers used for RT-qPCR, PCR, and cloning.

Gene fragment sequences for cloning.

List of deposited data and accessions.

DNA sequencing and reference genomes.

RIL population generation, GBS sequencing, and genetic map estimation.

Stable Agrobacterium-mediated transformation of P. axillaris.

LC-UV and LC-MS chromatographic methods.

Also available on Dryad: https://doi.org/10.5061/dryad.jsxksn083.

Fasta alignment of MYB SG6 and PH4 clade.

Newick tree file for MYB SG6 and PH4 clade.

Fasta alignment of entire MYB tree.

Newick tree file for entire MYB tree.

Normalized read counts for flavonoid and phenylpropanoid gene expression across Petunia species.

Click here for additional data file.