Dr. Mark Bassett, University of Florida, developed a large set of backcross introgression lines using 5-593 (PI 608674; black seed, purple flower) as the recurrent parent and donor lines with specific alleles that affect flower and seed coat color and pattern (Bassett, 2007). 5-593 carries the dominant allele for all but two of the genes. Its genotype for these genes is: T P [C r] J G B V Rk Gy sal. Introgression and donor lines used in this research are listed in Table 1. The nomenclature for the lines begins with the introgressed allele, followed by the level of backcrossing to 5-593. For example, g v BC₃ 5-593 means the recessive g and v alleles were introgressed to the BC₃ generation with 5-93 as the recurrent parent. Other genotypes used for resequencing are described in Supplemenatry Table S1. An F₂ population [g v BC₃ 5-593 (brown seed, white flowers) x Black Magic black seed, purple flower; n = 120] was developed and scored for flower color.

See Supplementary Table S2 for chemical sources.

Seed coat powders were prepared as described previously (Hart et al., 2020). Freeze-dried flower samples were ground into powder. The powders were weighed and initially solubilized with a 50% solution of methanol in water and dried using a Labconco vacuum concentrator (Kansas City, MO). Samples were reconstituted with 2 ml (flower) and 3 ml (seed coat) of 50% dimethyl sulfoxide solution in water. Standard curves were also prepared in 50% DMSO.

Samples were analyzed via reversed phase UPLC on a Waters (Milford, MA) Acquity H-Class system equipped with a QDa single quadrupole mass spectrometer and an Acquity (Waters) 2.1 × 100 mm bonded ethyl hybrid (BEH) C₁₈ column packed with 1.7 μm particles. The QDa was calibrated across the m/z range from 30 to 1,250 using the onboard calibrant. Mobile phase A consisted of water, phase B was 1% formic acid in water, and phase C was methanol. A flow rate of 0.5 ml min⁻¹ was used and mobile phase B proportion remained at 10% throughout all gradients. Infusion experiments of individual polyphenol solutions into mobile phase flow to the QDa were used to determine optimal source and polarity settings for each compound, which were then assigned to either a positive or negative polarity panel. For each panel, selected ion recording channels corresponding to the optimal m/z value and determined source settings were established for each constituent polyphenol. Panels and settings information are described in Supplementary Tables S3, S4. The positive ion panel used a gradient of 5 to 30% C in 3.5 min, 30 to 60% C in 5.5 min, 60 to 90% C in 0.5 min, then a return to initial conditions in 0.5 min followed by a 4 min equilibration period. The negative panel utilized a gradient of 5 to 30% C in 3.5 min, 30–40% C in 1.9 min, 40–90% C in 0.35 min, then a return to initial conditions in 1 min, followed by a 3.75 min equilibration period. For the positive ion panel, a four-point standard curve was prepared with each point containing an equimolar mixture of all polyphenols across the range from 0.108 to 108 pmol μL⁻¹. For the negative ion panel, a five-point standard curve was prepared with each point containing an equimolar mixture across the range from 0.100 to 286 pmol μL⁻¹. 1 μl of all standard and bean seed coat and flower samples were injected.

At least 2 technical replicate injections were used to determine a mean concentration for all reported values. For samples with measured concentrations above the upper limit of quantitation, the corresponding sample was diluted to bring the measured concentration within range of the standard curve and the reported response for such values represents that determination multiplied by the dilution factor. Linear, non-weighted fits were used for all standard curves. The instrument control and quantitation functions were performed using Waters Empower 3 software.

DNA was isolated from either leaf or embryo tissue using the Mag-Bind Plant DNA Plus kit (Omega Bio-Tek).¹ The structure of the F3′5′H gene models from G19833 (Phvul.006G018800) and UI 111 (PvUI111.06G022100) were used to develop primers (Supplemenatry Table S5) for PCR amplifications. PCR fragments were amplified in a 25-μl volume using an amplification protocol with 45 cycles and annealing temperatures specific to each primer pair. DNA fragments were extracted with the NEB Monarch Gel Extraction Kit,² and the fragments were Sanger sequenced by Eton Bioscience Inc.³

Thirty genotypes with various seed coat and flower color were selected for draft genome sequencing and assembly (Supplemenatry Table S1). For linked read sequencing, DNA was isolated from leaf tissue using the GE Healthcare Illustra DNA preparation kit (RPN8510). DNA was labeled with barcodes using the 10X Chromium Controller. The labeled DNA was sequenced on an Illumina HiSeq X Ten sequencer. A draft assembly for each genotype was generated from ~275,000,000 linked reads using the Supernova assembler with default settings (Weisenfeld et al., 2017).

For each species, the most recent reference genome represented in the Phytozome 13 database (accessed March 3, 2021)⁴ was searched with the key words “K13083” and “K05280,” the KEGG⁵ identifiers for F3′5′H and F3′H proteins, respectively. All protein sequences annotated with these identifiers were downloaded. Secondly, an exhaustive blastp search was performed using the P. vulgaris 5-593 F3′5′H protein sequence as a query at the NCBI protein database. Each Angiosperm order defined by APG IV (2016) was searched. To capture non-angiosperm sequences, all higher order taxa within the Viridiplantae⁶ were used as the database. All sequences identified by the blastp analysis were screened for the critical P450 heme-binding domain site (PFGAGRRICAG) and substrate recognition site 6 (SRS6). F3′5′H and flavonoid 3′ hydroxylase (F3′H) proteins vary at SRS6 position 8: F3′5′H = A or S; F3′5′H proteins = T (Seitz et al., 2007). Two F3′5′H Asterales proteins evolutionarily derived from F3′H genes (Seitz et al., 2015) were also included. The sequences were further limited to those without a deletion in either region of the other P450 essential domains and SRSs. Only 14 genes were excluded, and in each case, the excluded gene model was a duplicate of a gene with all functional domains in that species. Finally, all protein sequences were reannotated relative to KEGG nomenclature using the BlastKOALA web server (https://www.kegg.jp/blastkoala/). Only those sequences classified with the identifier “K13083” were used for the phylogeny analysis (Supplementary Table S6).

A phylogenetic network was constructed using the median network approach with the nucleotide sequences for each V haplotype as implemented in SplitsTree (Huson and Bryan, 2006). The MUSCLE algorithm (Edgar, 2004) as implemented in the MEGA 7 package (Kumar et al., 2016) was used to align the full protein sequences. MEGA 7 was used to construct an unrooted 50% consensus maximum likelihood tree. The initial tree was defined by the Neighbor-Joining/BioNJ algorithms, and tree construction utilized the Jones–Taylor–Thornton substitution model. Evolutionary rates among sites were modeled as a Gamma distribution with five categories. The tree construction consisted of 500 bootstrap replicates.

The 3D structure of the 5-593F3′5′H protein was modeled using MODELLER (Webb and Sali, 2017), a comparative protein modeling program implemented at ModWeb.⁷ The transmembrane protein was predicted using the MINNOU method (Cao et al., 2006) available at http://minnou.cchmc.org/.

Multiple sequence alignments were developed using the T-Coffee server,⁸ and boxshade displays were created using the BoxShade server.⁹ Critical CYP450 domains (Du et al., 2016) were extracted from the multiple sequence alignment. The WebLogo frequency bit scores (Crooks et al., 2004) were displayed graphically using the WebLogo WWW server.¹⁰

Genotypes with a dominant V allele are purple-flowered, and with a dominant B or dominant G and B alleles, express black-seeded beans. Homozygous recessive individuals have white (v) or pink (v^lae) flowers, and their seeds lack dihydromyricetin-derived flavonoids. V was originally mapped to chromosome Pv06 in the 7–15 F₂ population (McClean et al., 2002; Figure 2A). RAPD marker OD12₈₀₀, located at position 9,288,093 bp in the extremely low recombination heterochromatic region of the chromosome, where the physical/genetic distance ratio is 2.1 Mb/cM (Figure 2B), co-segregated with V. Fifteen indel markers, located in the low recombination region, were polymorphic between the 7–15 parents (Moghaddam et al., 2014), and all co-segregated with V (Figure 2A). With such a low recombination rate in this region of Pv06, an alternative approach to fine-mapping V was necessary.

A detailed phytochemical analysis found that 5-593 (V) flowers and seeds contained myricetin 3-glucoside and delphinidin 3-glucoside (Table 2), two compounds that require F3′5′H activity. These compounds were absent from pink [PI 608680 (v^lae)] and white-flowered [PI 608670 (v)] introgression lines, and greatly reduced or absent in the seed of those lines, respectively. With this background knowledge, a search for flavonoid pathway gene models within the low recombination Pv06 region of the reference Andean genome G19833 was undertaken. Two gene models (Phvul.006G015400 and Phvul.006G018800; Figure 2A) with the F3′5′H KEGG identifier K13083, were found in the interval. The collective genetic and physical evidence that placed V in the low recombination region of Pv06 turned the focus to an analysis of natural variants of these two gene models in bean lines with different V alleles.

Primers were designed to amplify each exon of these two models, and the protein structure was determined from the CDS sequence. The Phvul.006G015400 amplicons from the two parents of the 7–15 F₂ mapping population (McClean et al., 2002) were identical, and the eighth amino acid in the SRS6 domain was threonine which defines F3′H, not F3′5′H, function (Seitz et al., 2007). Therefore, Phvul.006G015400 was excluded as a V candidate.

The 5-593 (V) Phvul.006G018800 CDS was assembled from amplicon sequences, translated, and modeled as a member of the CYP450 protein family. While family members perform many biological functions, their sequences can be as little as 20% identical. Yet, their three-dimensional topology is highly conserved (Graham and Peterson, 1999). The topology consists of two neighboring α-helices clusters that interact to form a structural core to which the heme binds, and β-sheets 1 and 2 that provide the substrate access point. Four conserved motifs (I-helix, K-helix, ERR-triad, and heme-binding) are associated with these structural domains. Six substrate recognition sites (SRS1-6) are nearly universally found in CYP450 proteins and are important for substrate binding and the enzymatic reaction (Gotoh, 1992). Sequence analysis of the 5-593 Phvul.006G018800 protein identified the six conserved substrate recognition sites (SRS) and the four functional motifs typical of CYP enzymes (Figure 3). These domain sequences, except for the transmembrane domain, were nearly identical to those found in functional F3′5′Hs of other legumes (Supplemenatry Figure S1). The CDS from three black-seeded, purple-flowered Andean genotypes (Black Wonder, RH No. 6, W6 9644; Table 1) were identical to the Middle American 5-593.

The CDS sequences of Phvul.006G018800 from 5-593 and the reference G19833 differed by seven SNPs and a single cytosine insertion. The G19833 insertion led to sequence divergence of 28 amino acids beginning at amino acid 358, and an early stop codon that truncated the protein by 125 amino acids. This truncation eliminated the k-helix, SRS5, ERR-triad, heme-binding domain, and SRS6 (Figure 3). The Phvul.006G018800 sequences for multiple homozygous dominant purple-flowered plants of the 7–15 population were identical to 5-593, while the white-flowered plants sequence was identical to the ICA-Calima sequence. The Phvul.006G018800 sequence variation between 5-593 and G19833 was further associated with V by genetically testing three polymorphic KASP markers located in the coding region (Supplementary Table S7) on a F₂ population developed by crossing the purple-flowered Black Magic parent, whose Phvul.006G018000 sequence is identical to 5-593, and white-flowered introgression line g v BC₃-5-593. All three markers co-segregated with flower color.

Phvul.006G018800 was next sequenced from V0491 (PI 527745) and v^lae BC₃ 5-593, two genotypes with pink flowers and mineral brown seeds. Their sequence differed from 5-593 by two SNPs and one amino acid substitution (Figure 3). The R428M amino acid substitution is located near the critical PERF motif that stabilizes the heme in P450 enzymes (Du et al., 2016). M0056, the v donor for the v BC₃-593 introgression line, with mineral brown seed coats and pink flowers, differed from 5-593 by one SNP and a 14-nucleotide deletion in exon three that changed 14 amino acids beginning at amino acid 416 followed by the introduction of a frameshift stop codon that resulted in the deletion of the terminal 97 amino acids (Figure 3, Supplementary Table S8). The deletion altered the ERR-triad and deleted the heme-binding domain and SRS6.

A blastp analysis identified PvUI111.06G022100 as an ortholog of Phvul.006G018800 in the high quality, long-read race Durango pinto UI 111 common bean reference assembly.¹¹ This pinto genotype has white flowers and a brown mottle seed coat pattern on a cartridge buff background and was previously determined to carry the recessive v allele (Prakken, 1977). A nucleotide (T) deletion at position 1,323 relative to 5-593 was detected in the UI 111 CDS that altered the protein sequence starting at amino acid position 441 and introduced an early stop codon 16 amino acids further downstream that eliminated the essential heme-binding domain and SRS6 (Figure 3). Most recently, Labor Ovalle, a race Guatemala landrace with black seeds and purple flower was sequenced with long-read technology.¹² Gene model PvLabOv.06G020400 was found to be identical to the 5-593 protein sequence. This result reiterates the close relationship between the Guatemala and Mesoamerica common bean races (Tobar Pinon et al., 2021).

Draft genome assemblies were developed for 30 P. vulgaris genotypes, with varying flower and seed coat phenotypes, by sequencing 10X linked read libraries and assembling scaffolds. The assemblies ranged from 470 Mb to 570 Mb in size, with an average of 523 Mb. The contig N50 ranged from 26 Kb to 88 Kb with an average of 70.29 Kb. The scaffold N50 ranged from 32 Kb to 4.5 Mb with an average of 1.1 Mb (Supplementary Table S1). Each genotype with white flowers and seeds was homozygous for a recessive p allele. The gene model associated with the V gene was determined for each genotype by homology to the G19833 Phvul.006G018800 and UI 111 PvUI111.06G022100 gene models. Four additional alleles were identified. White-flowered Mex235 contained a 14 nt deletion, white-flowered G2858 contained a 15 nt deletion, pink-flowered MDRK contained an additional four nucleotide insertion relative to V0491, and pink-flowered G5686 had an additional SNP relative to MDRK. Among all genotypes screened, a total of nine P. vulgaris F3′5′H CDS haplotypes were discovered (Figure 3, Supplementary Table S8), and the alleles were annotated with a superscript designation of the reference genotype.

A phylogenetic network shows a direct mutational relationship between all haplotypes with the dominant V^[5-593] as the hub haplotype (Figure 4). This is expected since this haplotype is shared by the two gene pools. Three recessive white-flowered haplotypes were derived directly from V^[5-593]. The v^lae-[V0491] haplotype links the Middle America gene pool to the Andean gene pool haplotypes v^lae-[MDRK], v^lae-[G5686], and v^[G19833] through SNPs and a single nucleotide insertion. Nucleotide and amino variants for each haplotype are found in Supplementary Table S8.

The phylogenetic network revealed an important relationship between the F3′5′H coding region of the v^[Mex235], v^lae-[MDRK], and v^[M0056] haplotypes which provides evidence that V indeed encodes F3′5′H. The explanation begins with the F₂ population from the cross of 5-593 (V^[5-593] sal) and M0056 (v^[M0056] Sal) described previously (Bassett et al., 1990). In this population, individuals homozygous recessive for Sal (dominant Sal is epistatic to all V alleles), segregated 3:1, purple to white flowers. Therefore, for white flowers to appear, M0056 must possess a non-functional v gene. When the F3′5′H coding regions for the v^[Mex235], v^lae-[MDRK], and v^[M0056] haplotypes are aligned, an intragenic crossover event is observed between the v^lae-[MDRK] and v^[Mex235] haplotypes that generated the v^[M0056] haplotype (Figure 4). This natural recombinant haplotype contains the 5′ end of the v^lae-[MDRK] haplotype, and the 3′ end of the v^[Mex235] haplotype with the 14 nt deletion that eliminates the heme-binding domain and SRS6. Therefore, a natural common bean recombinant variant, generated by an intragenic crossover event between a chromosome that expressed white flowers and a chromosome that expressed pink flowers, contains a deletion that eliminated the F3′5′H function of the pink haplotype. This is equivalent to gene editing, where if a portion of a candidate gene is deleted and a phenotypic change occurs, it is concluded that the candidate gene is indeed the gene of interest. Therefore, the V gene encodes the F3′5′H protein and provides the important function of enabling the production of dihydromyricetin-derived flavonoids.

V orthologs in other legumes were discovered by a blastp search. The classic soybean (Glycine max L.) W1 flower color gene, that imparts purple color and encodes a F3′5′H (Zabala and Vodkin, 2007), was the top ortholog. The white flower w1 allele contains a tandem repeat that introduces an early stop codon that eliminates the critical C-terminal SRS6 domain. The light purple flowers of the G. soja w1-lp allele have greatly reduced levels of dihydromyricetin-derived flavonoid glucosides (Takahashi et al., 2010). A V210M substitution in the third amino acid of SRS2 in the w1-lp protein occurs in a residue that is invariant among all legume F3′5′H proteins that produce purple flowers. The B locus in pea (Pisum sativum L.) encodes F3′5′H (Moreau et al., 2012), and the recessive pink-flowered b allele contains a 23 nt deletion that introduces a premature stop codon that eliminates part of SRS1 and all other functional domains. The mutant flowers lack delphinidin and petunidin derivatives found in the wild-type genotype. A single nucleotide change was detected in a second pink-flowered b mutant that resulted in a single amino acid change (G → E) of glycine residue five in SRS1 that is invariant in all legume F3′5′H proteins (Supplementary Figure S1). The importance of this amino acid was further demonstrated in sweet pea (Lathyrus odoratus) where a pink-flowered mutant of the classic A1 gene (Punnett, 1923), which also encodes F3′5′H (Xue and Cronk, 2017), contains the same G → E SRS1 substitution observed in the P. sativum mutant. Finally, the P. vulgaris 5-593 F3′5′H protein and its ortholog in the black-seeded tepary bean (P. acutifolius A. Gray; Moghaddam et al., 2021) were 96% identical with only a single amino acid substitution in the SRS1 domain (Supplementary Figure S1).

To assess the distribution, structural variation, and phylogeny of F3′5′H across the breadth of land plants, V gene orthologs were mined from the Phytozome and NCBI genome databases. The protein sequences were selected based on functional features repeatedly shown to be necessary for the F3′5′H enzymatic reaction (Graham and Peterson, 1999). This whole-genome, computational approach, based on domains experimentally proven to be critical to the function of a gene, are preferred to simple homology searches (Restrepo-Montoya et al., 2020). This exhaustive search identified 177 F3′5′H protein coding sequences from 129 land plant species. The gene was not detected in Chlorophyta. Within the Embryophyta, the gene was not detected in ferns or bryophytes. The presence of F3′5′H in the gymnosperms (Figure 5A) dates the gene back to the appearance of Spermatophyta (348 Ma; Morris et al., 2018). The F3′5′H gene is found in the Nymphaeales basal angiosperms and sporadically throughout other angiosperm orders. The gene was absent in the Magnoliales. A search of the 129 species with a whole-genome reference assembly and annotation in the Phytozome database that represents a broad taxonomic collection of the Viridiplantae, found only 53 species contained a F3′5′H gene model with all essential domains. The gene was discovered in all monocot species and nine of the ten Fabales species (except peanut, Arachis hypogaea). It was absent in all 28 Phytozome members of the Brassicaceae family. When the NCBI and Phytozome sequences were combined, F3′5′H was found in the genome of 29 of the 64 Angiosperm orders (APG IV, 2016) and three of the six Gymnosperm orders (Figure 5A; Forest et al., 2018).

The 177 F3′5′H proteins were aligned, and the highly conserved CYP450 motifs and SRSs were visualized (Figures 6A,B). SRS1 is the most variable region with 86 unique sequences (Figure 6B). The amino acids 245H, in SRS3, and 302N in SRS4, previously shown to be under positive selection (Jia et al., 2020), were nearly invariant among all species. SRS4, which extends across helix I and is an important component of the functional pocket of P450s and associated with oxygen binding, is well-conserved. Among all species, the SRS4 I-helix motif was invariant. Only 10 SRS5 sequences and 21K-helix motif variants were observed. This conservation is not unexpected since amino acid variants in these two adjacent features can alter P450 hydroxylation (Richardson and Johnson, 1994) and stereoselectivity in humans (Ellis et al., 1996), as well as substrate binding orientation in plants (Schalk and Croteau, 2000). Similarly, a single substitution, Ala/Ser vs. Thr, at amino acid position nine in SRS6, defines the 3′5′ versus 3′ hydroxylation of dihydrokaempferol (DHK), respectively (Seitz et al., 2007).

A 50% consensus maximum likelihood gene tree was constructed with the full F3′5′H protein sequence collection (Figure 5B). A single clade with a bootstrap support value (BSV) = 78 revealed a relationship between Gymnosperms, Nymphaeles, and a subset of Poales proteins. Within this clade, strongly supported subclades for each phylogenetic group were detected with BSV values of 97, 99, and 99, respectively. Unlike this group of land plants, a single eudicot clade was not detected. Rather, all species in multiple orders (Cornales, Fabales, Fagales, Gentianales, Lamiales, Malvales, Rannuncales, Rosales, and Solanales) formed a clade. For the Fabales species, which are all members of the Papilionoideae subfamily, the F3′5′H gene tree is consistent with a genus-level legume species tree based on multilocus sequence data (Koenen et al., 2020).

The F3′5′H duplication history was considered for the 25 Phytozome reference genomes assemblies which have multiple copies of the gene. Duplicated copies could be traced to polyploidization events, such as with switchgrass (Panicum virgatum), where duplicates were located on chromosome 9 of the K and N subgenomes (Lovell et al., 2021). With Eucalyptus grandis, the duplicates on chromosomes 2 and 10 were the result of the whole-genome duplication, while an additional copy on chromosome 2 appeared via tandem duplication. These types of duplication events are representative of many E. grandis genes (Myburg et al., 2014). Interpreting the duplications found in the Poales is more challenging. Except for Panicum virgatum, as noted above, all Poales/PACMAD clade species only contain a single F3′5′H gene copy. By contrast, all Poales/BOP species contain multiple copies. These copies were found in two well-supported clades, Poales/BOP I and II (Figure 5B). The copies in Poales/BOP II, the smaller of the two clades, are not on the same chromosome as those genes found in Poales/BOP, suggesting a different evolutionary arc for these genes.