A total of 182 dry bean lines from the Durango Diversity Panel (DDP), described in a previous study by Soler-Garzón et al. (2021), were used in this study to map and target candidate genes for bc-2/bc-2² resistance alleles and the unknown gene that interacts with bc-2. The Othello/VAX 3 (OV3) F_6:7 recombinant inbred line (RIL) population (n = 100) (Viteri et al., 2014) was used to narrow the bc-2² genomic interval. Othello pinto bean cultivar (NW-410/2/Victor/Aurora; Burke et al., 1995) possesses the bc-1/bc-2²/bc-u^d resistance genes (Soler-Garzón et al., 2021), which places it in HG-6. Othello is only susceptible to BCMV PG-VII strains such as US-6 and NL-4. VAX 3 is a germplasm release derived from a multiple-parent interspecific cross between common bean and tepary bean (Phaseolus acutifolius). VAX 3 has the I gene for resistance to BCMV (Singh et al., 2001).

Seven F₂ populations were generated from crosses among DDP lines selected for presence vs. absence of the bc-2 allele based on phenotypic response to different PGs, as reported in the literature, or from inoculations conducted by us. Furthermore, two HG-4 navy bean cultivars, namely, Sanilac and Michelite-62 with bc-2, were crossed with Poncho (DDP-041), Beryl (DDP055), and Matterhorn (DDP-033), which lack bc-2 to generate six additional F₂ populations segregating for the bc-2 gene originating from race Mesoamerica (Supplementary Table 1).

The DDP lines were evaluated for reaction to US-6 (PG-VII) BCMV strain (Silbernagel, 1969). DDP reactions to NL-8 (PG-III) and NL-3 (PG-VI) BCMNV strains were reported by Soler-Garzón et al. (2021). The OV3 RILs were evaluated for reaction to US-6, NL-8, and NL-3 strains and the BCMV US-4 (PG-IV) strain (Skotland and Burke, 1961) was used to inoculate a subset of OV3 RILs. The subsets of individual plants from F₂ populations and F_2:3 families (i.e., progeny lines) were evaluated for reaction to one or more strains, which included US-4 (PG-IV) and US-6 (PG-V11) BCMV strains as well as NL-8 (PG-III) and NL-3 (PG-VI) BCMNV strains.

All virus strains were maintained in the universal susceptible HG-1 differential cultivar “Sutter Pink,” with the youngest leaves with mosaic symptoms used as inoculum. All inoculations were conducted in the USDA-ARS greenhouses at Prosser, WA, United States, under controlled conditions (22–28°C temperature range and 14-h photoperiod using artificial lights as necessary). Test plants were sown in 9 cm³ pots (three plants maximum per pot) containing a commercial potting mix (Sun Gro Horticulture, Bellevue, WA, United States). About 10 days after planting, the two primary leaves of each plant were mechanically inoculated with the same strain using the procedure described by Drijfhout (1978).

Disease reactions were observed every week until 5 weeks post-inoculation (wpi). The symptoms for each plant, as reported by Soler-Garzón et al. (2021), were recorded as follows: NS = no apparent symptoms on inoculated leaves and no systemic symptoms; M = leaf curling and plant stunting with severe systemic chlorotic mosaic symptoms; mM = mild systemic chlorotic mosaic symptoms; dM = delayed severe systemic chlorotic mosaic symptoms observed 2–4 wpi; VN = restricted vein necrosis on inoculated leaves, no systemic symptoms; VN⁺ = restricted vein necrosis on inoculated leaves, with some small patches (10 mm²) of systemic restricted vein necrosis on upper trifoliolate leaves observed from 3 to 5 wpi; NLL = local necrotic pinpoint lesions on inoculated leaves, no systemic symptoms; TN = lethal systemic top necrosis by 7–10 days post-inoculation (dpi), resulting in plant death; and dTN = delayed top necrosis beginning >11 dpi, most often resulting in plant death. Host differentials with specific reactions to BCMV and BCMNV strains in parenthesis (US-4/US-6/NL-8/NL-3): HG-1, Sutter pink (M/M/M/M); HG-4, UI-34 (NS/NS/M/M); HG-6, Othello (NS/M/NS/NS); HG-9, Topcrop (NS/NS/VN/TN); HG-10, Beryl (NS/NS/VN/VN); and HG-9b, Jubila (TN/NS/VN/VN⁺) were included as controls.

Genomic DNA was extracted from 20 mg of leaf tissue from an individual plant for each DDP accession, check line, and HG differential cultivar using a QIAGEN DNeasy 96 Plant Kit (Hilden, Germany). For F₂ and F₃ plants and single plants of OV3 RILs, DNA was extracted from four-leaf disks (approximately 30 mm²) according to the alkaline extraction method described by Xin et al. (2003) with modifications (Soler-Garzón et al., 2021). A 1:7 dilution of the extracted DNA was placed in a 96-well plate with a final volume of 100 μl. Finally, 5 μl of diluted DNA template was used for PCR.

The DDP was genotyped, as reported by Soler-Garzón et al. (2021), with a filtered set of 1,269,044 biallelic single-nucleotide polymorphisms (SNPs) obtained by whole-genome resequencing. The linkage map for OV3, which was generated with polymorphic SNPs from the BARCBean6K_3 BeadChip assay (Song et al., 2015), was obtained from Viteri et al. (2014). The DDP and OV3 RILs were assayed for the I-, bc-1-, bc-3-, and bc-u^d-linked SNP markers described by Soler-Garzón et al. (2021). The individual F₂ plants from each population were similarly assayed with the same resistance gene-linked markers. Eventually, the abovementioned materials were assayed with the new putative markers for bc-2 and unknown gene (described below). Notably, 35 select F_2:3 progenies from individual F₂ plants, either fixed or heterozygous for one or more of the resistance gene-linked markers, were genotyped using this same set of SNP markers. Eight true-breeding F_2:3 families were used as controls for each gene in a homozygous state, and the remaining 27 families had at least 1 gene in a heterozygous state (Supplementary Table 1). A total of 217 F₂ plants and 1,345 F_2:3 plants were assayed for the I, bc-1, bc-3, bc-u^d, and new bc-2 and unknown gene, i.e., gene-linked SNP markers (Supplementary Table 2).

The phenotypic reactions to US-6, NL-8, and NL-3 strains and the genotypic SNP data for the DDP were integrated for the genome-wide association study using a multi-locus random-SNP-effect mixed linear model (mrMLM) described by Wang et al. (2016) and implemented in the “mrMLM” R package (Wen et al., 2018). A kinship matrix was generated using the efficient mixed-model association (EMMA) algorithm implemented in the Genome Association and Prediction Integrated Tool (GAPIT) R package (Lipka et al., 2012) with corrections for kinship and population structure. Five principal components (PCs) generated from GAPIT were included as covariates. The Bonferroni test was implemented to control the experiment-wise type I error rate at 0.05. The GWAS results were plotted using CMPlot v3.62 (Yin, 2016), and IntAssoPlot v3 (He et al., 2020) was used to represent regional-based marker-trait associations graphically.

Candidate genes located within the identified genomic intervals were obtained by perusing available scientific literature and aligning their sequences to the Andean G19833 v2.1 and Middle-American UI 111 v1.1 P. vulgaris reference genomes.¹ Subsequently, exon sequences of each candidate gene were amplified across select genotypes using PCR primers designed with Primer3 software (Koressaar and Remm, 2007; Untergasser et al., 2012; Supplementary Table 3).

Two standard PCRs, for each sample, were replicated in a volume of 25 μl containing 1.8 mM MgCl₂, 0.4 mM of deoxynucleoside triphosphate (dNTP) Mix (Promega™, Madison, WI, United States), 0.25 μM of each primer (forward and reverse), 25 ng of genomic DNA, and 1 unit Taq DNA polymerase (Promega), in 1× PCR buffer (Promega), under the following amplification conditions: 2 min at 95°C, followed by 38 cycles at 94°C for 20 s, specific annealing temperature for each primer set for 30 s and 72°C for 90 s, and a final extension at 72°C for 5 min. All amplifications were performed in a PCR Eppendorf Mastercycler (Eppendorf AG, Hamburg, Germany). PCR fragments were visualized by gel electrophoresis on 2% (w + v) agarose and then purified and Sanger sequenced by Eurofins MWG Operon (Louisville, KY, United States). Sequence trimming, alignment, and polymorphism discovery were performed with Geneious 9.1.2 software (Kearse et al., 2012).

For each polymorphic variant identified in coding regions of the sequenced candidate genes, a set of allele-specific primers, with GC tail of unequal length attached to their 5′ end, were designed using the Primer3 software (Koressaar and Remm, 2007; Untergasser et al., 2012), according to melting temperature (Tm)-shift SNP genotyping method developed by Wang et al. (2005). Fragments were amplified by PCR on an Eppendorf Mastercycler using a volume of 20-μl master-mix containing 1.5 mM MgCl₂, 0.2 mM of dNTP Mix (Promega), 0.15 μM of each primer, 1X EvaGreen™ (Biotium, Fremont, CA, United States), 20 ng of genomic DNA, 1X Taq buffer, and 0.1 μl Taq1 polymerase (Promega) under the following profile: 94°C for 2 min, then 38 cycles of denaturation at 92°C for 20 s, annealing for 20 s (i.e., the temperature was specific to each primer trio), extension at 72°C for 20 s, and a final extension at 72°C for 5 min. The melting point analysis for the allele determination of the template DNA was performed with a fluorescence-detecting thermocycler (LightCycler™ 4890 Instrument II, Roche, Basel, Switzerland) with EvaGreen™ Fluorescent Dye (Biotium). The fluorescent detection profile was for 1 min at 95°C, and the melting curve step was ramping from 65 to 95°C in increments of 1°C for every 20 s.

The GWAS was performed with 23 DDP lines with bc-2 or bc-2² alleles [i.e., 9 lines with bc-2 (NS to US-6 but M against NL-3) and 14 lines with bc-2² (12 lines with M against US-6 but NS to NL-3, and 2 lines with NLL to NL-3)] compared to 146 DDP lines absent the bc-2 and bc-2² alleles as indicated by M to both US-6 and NL-3, or mM, TN, or VN to NL-3. One significant peak interval (α = 0.05; p-value = 8.53E-26 to 4.33E-31) was detected between 9,117,985 and 10,468,271 based on Pv11 in G19833 v2.1 reference genome (Figure 1A and Supplementary Table 4a).

The interval for bc-2² was narrowed from 9.16 to 9.67 Mb in the OV3 RIL population, based on three RILs (i.e., OV3-92 with VN, OV3-107 with TN, and OV3-95 with NLL reactions to NL-3) with recombinant SNP haplotypes (Figure 1B and Supplementary Table 5). The narrowed region for bc-2², when applied to the whole-genome resequencing data for the DDP, revealed a low level of linkage disequilibrium (LD) between 74 single-nucleotide variants (SNVs), 31 insertions/deletions (InDel), and 11 short tandem repeats (STRs) polymorphic variants, distributed across 32 gene models (based on the G19833 v2.1 reference genome) (Figure 1C and Supplementary Table 6).

Additional alignment of the bc-2 region between G19833 v2.1 and UI 111 v1.1 reference genomes revealed a gap of 10,084 bases in “UI 111,” which just so happens to possess the bc-2 allele based on NS reaction to US-6 strain herein and from published reports (Singh et al., 2007, pp. 987). This 10-kb deletion in “UI 111” spans the Phvul.011G092600 and Phvul.011G092700 gene models of G19833 v2.1 reference genome (Figure 2A). The Phvul.011G092600 model, a subtilisin-like serine protein with endopeptidase activity, is completely missing, and the 3′ end of Phvul.011G092700, a vacuolar protein-sorting 4 (Vps4) AAA+ ATPase endosomal sorting complexes required for transport (ESCRT) protein involved in multivesicular endosome function, is truncated (Figures 2B,C). Amano et al. (2013) identified a Vps4 AAA+ ATPase ESCRT protein as a candidate gene (Cucsa.385040 in Phytozome v13) for the recessive zym^{A192– 18} gene conferring resistance to Zucchini yellow mosaic virus (ZYMV) in cucumber (Cucumis sativus L.). Using BLASTp, the Vps4 AAA+ ATPase ESCRT protein in cucumber had 92.2% identity with the Phvul.011G092700 gene model. Based on functional similarity to zym^{A192– 18}, Phvul.011G092700 was chosen as the most likely candidate gene for bc-2.

To identify the unknown gene that interacts with bc-2, GWAS was conducted with 9 DDP lines with bc-2 as detected by NS to US-6 and M to NL-3 strain compared to 160 DDP lines without bc-2 as determined by M, mM, NS, TN, and VN reactions against NL-3 strain. Notably, DDP lines without bc-2 included 12 lines with NS and 2 lines with NLL to NL-3 strain due to the presence of bc-2². The Manhattan plot revealed one significant peak (α = 0.05; p-value = 6.74E-43 to 2.62E-32) on Pv05, with an interval between 35,471,758 and 36,225,550 bases (Figure 3A and Supplementary Table 4b), which is putatively associated with the “unknown” gene that assists the expression of bc-2.

The Vps4 AAA+ ATPase ESCRT protein of the Phvul.011G092700 candidate gene for the bc-2 resistance gene was blasted against the G19833 v2.1 reference genome, and the most significant hit (i.e., 94.9% identity) was Phvul.005G125100, a Vps4 AAA+ ATPase ESCRT homolog. Phvul.005G125100 was within the 35,471,758- to 36,225,550-bp GWAS interval. A significant SNP S05_36225550 (p-value = 2.62E-32), a missense variant, from the GWAS analysis was located within Phvul.005G125100. Furthermore, SNP S05_36225550 was present in all nine DDP lines with the unknown gene and absent in all DDP lines without the unknown gene. Interestingly, among DDP lines, the SNP S05_36225550 was found to be 111,034 bases upstream and in high LD with the bc-u^d gene that interacts with bc-1 to confer resistance to BCMNV (Soler-Garzón et al., 2021; Figure 3B). Given these preliminary results, this new gene was assigned the preliminary gene symbol bc-4.

A total of 16 lines were chosen according to BCMV and BCMNV reactions, gene combinations, and origin, for the exon sequencing of the bc-2 candidate gene Phvul.011G092700 and the bc-4 candidate gene Phvul.005G125100.

The Phvul.011G092700 sequences were blasted to the G19833 v2.1 reference genome (which has the Bc-2 wild-type genotype), revealing eight exon regions (Table 1 and Figure 2A). The Phvul.011G092700 sequences in the race Durango genotypes with either bc-2 (i.e., UI-111 and UI129) or bc-2² (i.e., Othello and 92US-1006) recessive resistance alleles all possessed the same 10,084-bp deletion that truncates three exons from the 3′ end (Figure 2B), providing initial evidence of a single resistance allele. For the race Mesoamerican genotypes Robust, Michelite-62, and Sanilac, also with the bc-2 resistance allele, there was a single deletion-frameshift mutation (Pv11: 9,278,764 bases) in the first exon of Phvul.011G092700. To track these different bc-2 mutations, they were assigned a superscript in brackets: bc-2^{[UI 111]} for the large 10-kb deletion of three exons found in race Durango lines, and bc-2^[Robust] for the deletion-frameshift found in navy beans.

The G19833 Phvul.011G092700 protein is 446 amino acids (aa) in length and contains three domains, namely, (i) microtubule interacting and transport (MIT) domain, (ii) AAA-type ATPase (AAA+) domain, and (iii) Vps4 domain (Figure 2C). The same protein homolog PvUI111.11G095000 (but with a different code) in reference genome “UI 111” (which has bc-2 genotype) has only 319 aa and 2 domains, i.e., MIT and AAA+, due to the 10,084-bp deletion. The Sanilac (or Robust) homolog has only 128 aa due to 1 deletion-frameshift mutation generating a single aa replacement in aa position 53 (p.E53R) and a premature stop codon in position 77 (fsX77), resulting in an incomplete protein with only the MIT domain.

The in silico-translated Vps4 AAA+ ATPase ESCRT protein of sequenced open reading frame (ORF) for bc-4 candidate gene Phvul.005G125100 did not reveal disruption among the MIT, AAA+, and Vps4 domains. The SNP S05_36225550 missense variant was found to be a transversion T > G, which causes the substitution of a non-polar, neutral aa, Methionine (M), with a basic polar, positively charged aa, Arginine (R), at position 33 of the MIT domain (i.e., p.M33R) (Table 1 and Figure 4).

Tm-shift assays were developed to detect the polymorphisms in the bc-2 and bc-4 candidate genes (Table 2). Two distinct assays were designed for the different bc-2 candidate gene mutations. The bc-2^{[UI 111]} assay was developed using the first forward primer aligned to the “UI 111” gene (PvUI111.11G095000) sequence and the second forward primer aligned to the G19833 gene (Phvul.011G092700) sequence and a common reverse primer (Supplementary Material 1). The technical name for this marker assay is “Pvvps4_del.” For better amplification efficiency, no GC tails were added to “Pvvps4_del” forward primers. The assay for bc-2^[Robust] flanking the S11_9278764 in-frame-deletion in the navy beans Sanilac and Michelite-62 is named “Pvmit-2_C_del.”

The assay for bc-4 was developed using the S05_36225550, i.e., missense SNP, in the Phvul.005G125100 candidate gene. The technical name for this marker is “Pvmit-1_T_G.” These bc-2 and bc-4 markers were added to I, bc-1, bc-3, and bc-u^d markers developed by Soler-Garzón et al. (2021) for genotyping the lines and populations described below.

The DDP was assayed for both “Pvmit-2_C_del” and “Pvvps4_del” markers for bc-2, but only the latter bc-2^{[UI 111]}-linked marker was detected. All DDP lines with bc-2 (NS to US-6 and M to NL-3) or bc-2² resistance alleles (M to US-6 and NS to NL-3) possessed the “Pvvps4_del” marker for bc-2^{[UI 11]}, and all other DDP lines possessed the wild-type Bc-2 allele. The DDP lines with bc-2 possessed the bc-4 marker, and those with bc-2² possessed the bc-u^d marker, including those DDP lines in HG-11 (92US-1006 and Quincy) with NLL reaction to NL-3 conditioned by I, bc-2², and bc-u^d genes (Supplementary Table 7). These results further support that only one resistance allele, i.e., bc-2, exists at the Bc-2 locus. This bc-2 allele exhibits a differential effect based on the presence of bc-4 or bc-u^d.

Four genes (i.e., I, bc-1, bc-2, and bc-u^d) were segregated in the OV3 RIL population. Of the four genes, only the I gene had distorted (df = 1.0, p ≤ 0.05: X² = 5.99) segregation, 3:1, in favor of the dominant I gene, whereas 1:1 segregation was expected. Altogether, there were six distinct phenotypes to NL-3, seven phenotypes to NL-8, two phenotypes to US-6, and 16 distinct genotypes (Supplementary Table 8). The phenotypic reactions of all RILs to US-6, NL-8, and NL-3 strains matched the genotypes predicted by the resistance gene-linked markers. The most relevant results for this study were that RILs with bc-2 and bc-u^d expressed M to US-6 and NS to NL-3 when I gene was absent (∼HG-6) and NLL when I gene was present (∼HG-11), whereas the bc-2 allele had no discernible effect against any of the strains in the absence of bc-u^d. Additionally, OV3 RILs with unique genotypes were phenotyped with US-4 (PG-IV) to further validate their use as additional differential cultivars in HGs 1, 2, 6, 8, and 9b (Table 3).

Differential genotypes were distributed across 12 HGs for BCMV/BCMNV, and one unassigned genotype “UI 129” was assayed for the I, bc-1, bc-3, bc-u^d, bc-2^{[UI 111]}, bc-2^[Robust], and bc-4 markers (Table 3). Only genotypes in HG-4 and HG-5, previously reported to possess bc-2, possessed the bc-4 gene. In addition, HG-4 (NS to US-6 and M to NL-3) had a mix of genotypes with the two distinct bc-2 mutations, namely, Michelite-62 and Sanilac with bc-2^[Robust], and “UI 111” and “UI 34” with bc-2^{[UI 111]}. Moreover, all HG-4 genotypes exhibited NS to US-6 and M to NL-3 strains, regardless of the mutation.

The HG-5 genotype, “UI 114,” possesses bc-1, bc-2^{[UI 111]}, and bc-4. “UI 114,” similar to HG-4 genotypes, exhibits NS to US-6 and M to NL-3; however, it also possesses resistance to NL-8 strain, whereas HG-4 genotypes are susceptible. Soler-Garzón et al. (2021) showed that bc-1 alone provides protection to NL-8, and the addition of this gene in HG-5 lines explains the differential reaction between HG-4 and HG-5 genotypes to NL-8. Unassigned genotype “UI 129” also carries bc-1, bc-2^{[UI 111]}, and bc-4 like HG-5, but unlike HG-5 lines, it exhibits resistance to PG-V strains. Therefore, the hypothesis that “UI 129” has another unknown gene that contributes to its resistance against PG-V strains should be tested.

The HG-6 genotypes (bc-2 and bc-u^d) exhibit resistance to all BCMV and BCMNV PGs except PG-VII with BCMV strains such as NL-4 and US-6 (Drijfhout, 1978). However, most HG-6 genotypes possess bc-1, in addition to bc-2^{[UI 111]} and bc-u^d (e.g., Othello, UI-31, and UI-35). Furthermore, Monroe with bc-2^[Robust], bc-1, and bc-u^d is also classified as an HG-6 differential cultivar (Miklas et al., 2000), which indicates that the bc-2^[Robust] and bc-2^{[UI 111]} mutations react similarly with bc-u^d. Finally, the OV3 RILs with bc-2 and bc-u^d but without bc-1 were resistant to PG-III, IV, VI, and VII strains in our screenings, indicating that the combination of bc-2 and bc-u^d conditions resistance to PG-III and PG-V strains in the absence of bc-1.

The HG-7 cultivar IVT-7214 possessed the bc-3 and bc-4 genes, based on the marker assays, indicating that bc-4, not bc-u^d as previously considered, may be interacting with bc-3 to confer broad resistance against all PGs, except for the newly discovered PG-VIII (Feng et al., 2015). Therefore, we are in the process of developing a differential cultivar with bc-3 and bc-u^d to test whether bc-u^d may help to extend the range of effectiveness for bc-3 to condition resistance to this new PG-VIII as well. However, we did not observe our bc-2 gene-linked markers in IVT-7214, whereas its presence was reported by Drijfhout (1978). Perhaps IVT-7214 is heterogeneous, and the line we happened to test was absent bc-2. This is the best explanation because the results for the “Michelite-62 (bc-2/bc-4) x IVT-7214 (bc-2/bc-3/bc-4)” F₂ population from the study by Drijfhout (1978) support the segregation of a single gene, i.e., bc-3, with bc-2 and bc-4 genes present in both parents.

For the HGs 8–12 with I gene, bc-4 was absent in all, and bc-2 was only detected in HG-11. The HG-11 lines all possessed I, bc-2^{[UI 111]}, and bc-u^d genes. As was reported by Drijfhout (1978) and validated by Soler-Garzón et al. (2021), the bc-u^d and bc-2^{[UI 111]} gene combination is needed to condition NLL to BCMNV (PG-III and PG-VI) in the presence of the I gene. OV3 RILs with I, bc-1, bc-2^{[UI 111]}, and bc-u^d also expressed NLL to BCMNV, indicating no additional phenotypic effect from the presence of bc-1.