Abstract

• Plant domestication has strongly modified crop morphology and development. Nevertheless, many crops continue to display some atavistic characteristics that were advantageous to their wild ancestors, such as pod dehiscence (PD). Here, we provide the first comprehensive assessment of the inheritance of PD in a domesticated grain legume, common bean (Phaseolus vulgaris).
• We compared three methods to evaluate the PD phenotype and identified multiple, unlinked genetic regions controlling PD using a biparental population and two diversity panels, corresponding to the two major domestications of common bean. We subsequently assessed patterns of orthology between these loci and those controlling the trait in other species.
• Our results show that different genes were selected in each domestication and ecogeographic race. Foremost among these is a chromosome Pv03 dirigent-like gene involved in lignin biosynthesis, which shows a single base-pair substitution that is strongly associated with decreased PD and expansion of the Mesoamerican domesticates into northern Mexico, where the arid conditions promote PD.
• The  base-pair mutated in certain PD-resistant lineages has been conserved across species during millions of years of plant evolution. This contrast is a consequence of the markedly different fitness landscape imposed by domestication. Environmental dependency and genetic redundancy explain the maintenance of atavistic traits under domestication. Even core domestication traits in well-domesticated species show considerable genetic variation, which can ultimately affect food security.

Materials and Methods
 

Germplasm
A recombinant inbred (RI) population (n = 238), developed from a cross between ICA Bunsi (domesticated, PD-susceptible, Middle American) and SXB 405 (domesticated, PD-resistant, Middle American), was used for QTL mapping (Assefa et al. 2013; Berny Mier y Teran et al. 2019). For association mapping, different panels were used, the individual members of which are listed in Tables S2. Two-hundred eight members of the Andean Diversity Panel (ADP, Cichy et al. 2015) and 278 members of the Middle American Diversity (MDP, Moghaddam et al. 2016) were grown and phenotyped. Sequencing was performed in a diverse panel of 90 varieties representing six species were acquired from the National Plant Germplasm System. Eighteen varieties commonly grown at UC Davis with known PD phenotypes were also genotyped. Stringless snap bean varieties were specifically excluded from the analysis to avoid the epistatic effect of the Stringless (St) locus on PD.

Microscopy
Pods of G12873 (wild, high dehiscence), ICA Bunsi (domesticated dry bean, dehiscence-susceptible) SXB 405 (domesticated dry bean, dehiscence-resistant), and Midas (domesticated snap bean, dehiscence-susceptible) were Vibratome-sectioned to identify morphological differences that might be associated with PD. All sectioned pods were greenhouse-grown and harvested when pods were at full size with seeds filled, at the onset of pod color change. All sections were 100 micrometers thick and made in a transverse plane perpendicular to the fibers of interest. All sections were treated with Auramine O for at least 20 minutes to stain lignified tissue (Ursache et al. 2018). Fluorescence was visualized using an Olympus microscope. 

RI population cultivation and PD phenotyping
The ICA Bunsi/SXB 405 (IxS) RI population of 238 RILs was field-grown during the spring and summer of 2014. The spring planting was an un-replicated trial conducted at Coachella, California. At maturity, plots were visually evaluated for the presence or absence of PD, and the data were used as a phenotype for QTL mapping. During the summer of 2014, the RI population was grown in a replicated field trial in Davis, California. At maturity, dried non-dehiscing pods from 191 RILs were harvested from each plot; these were evaluated for susceptibility to PD by two methods. First, all pods were desiccated at 65°C for seven days, and then returned to room temperature for a minimum of seven additional days. The proportion of dehiscing pods after this process was recorded for each plot. Second, the amount of force required to induce pod fracture was measured using an Imada force measurement gauge (method modified from Dong et al., 2014). Force measurements were taken on pods that had not dehisced during the desiccation treatment. A bit mounted to the gauge was used to press the ventral side of each pod at the most apical seed, and the peak force required to cause fracture at the apical end of the pod beak was recorded. Force required for PD was normalized to account for small but significant differences between note-takers, and the standardized score was used for QTL mapping. Pods that failed to produce seeds were excluded from all phenotyping analyses.

Genotyping
Genomic DNA was extracted from parents and RILs of the IxS population using a modified CTAB protocol. DNA quality was confirmed using a NanoDrop spectrophotometer. The IxS population was genotyped using the Illumina Infinium II BARCBean6K_3 BeadChip (Song et al. 2015); 382 segregating SNPs were identified in the population. Primers spanning the transcribed sequence of Phvul.003G252100, also known as Phaseolus vulgaris Pod Dehiscence 1 (PvPdh1), a candidate gene underlying a major QTL identified in this study, were developed using the NCBI Primer-BLAST tool. Several differences in the genomic sequence exist between the Middle American and Andean gene pools, so a mixture of two forward primers was introduced into each PCR with a common reverse: PvPDH1 ALL Middle American Forward: CATCTCCCCCATTTTCCCCC; PvPDH1 ALL Andean Forward: CATCTCTCCCATTTTCTCCT; PvPDH1 ALL common Reverse: AACACGTGGAAGAGGAGGATT.  PCR conditions for this amplification included an initial denaturation at 95°C for 180s, 38 cycles of 95°C for 30s, 51°C for 30s, and 68°C for 60s, and a final elongation step of 68°C for 300s. Another set of primers was developed to specifically improve the amplification and sequencing of Andean common beans, with the sequences: PvPDH1 Andes Forward: TTTTTCTTGTGAGCAAAATTGAGTT; PvPDH1 Andes Reverse: GCAGAGGAAAAACACGTGGA. This primer set was amplified with an initial denaturation at 95°C for 300s, 34 cycles of 95°C for 30s, 46°C for 30s, and 72°C for 70s, and a final elongation step of 72°C for 300s. PCR products were cleaned using a GeneJET PCR Purification Kit and sequenced at the UC DNA Sequencing Facility by Sanger sequencing.

QTL mapping
Composite interval mapping was conducted using the R package R/qtl (Broman et al. 2003). Field dehiscence score, proportion dehiscing in a desiccator, and force measurements were separately used to identify PD QTLs marked by SNPs. The maximum LOD score of 1000 randomized permutations of the data was used as a significance threshold. Single QTL scans were performed using the scanone function. Multiple QTL mapping was conducted using the scantwo function in R/qtl and by running the analysis with RILs subset by genotype at the most significant marker near PvPdh1 on Pv03. QTL mapping results were based on maximum likelihood via the EM algorithm (Lander and Botstein 1989).

Validation of QTL mapping results using association mapping 
Two hundred and eight accessions of the ADP were grown in Davis, CA during summer 2016. PD in the field, proportion dehiscing in a desiccator, and force required for fracture were recorded. Principal component analysis was conducted on SNP data for the population, and the results were used as covariates to account for population structure. Two hundred seventy-eight members of the MDP were phenotyped for PD by desiccation in 2017. Association mapping was conducted using GLMs in TASSEL via SNiPlay (Bradbury et al., 2007; Dereeper et al., 2011). A minor allele frequency of 0.1 was used as a threshold for SNPs, and these SNPs were evaluated for significance based on a Bonferroni-corrected alpha of 0.05. QTL regions of significance were determined as the area between the first and last significant SNP on a chromosome arm. Individual significant SNPs without significant neighbors in the same population or others were not given further consideration, as these are likely All results were visualized using the qqman R package (Turner, 2018), including the Bonferroni-corrected significance thresholds at alpha=0.05 and 0.01 were shown, along with the positions of major candidate genes.

Expression and synteny mapping
Gene expression information from a variety of tissues and developmental stages was extracted from published data (O’Rourke et al. 2014) and visualized independently using R base graphics (R Core Team, 2013). Candidate genes related to PD were identified in significant QTL intervals based on definition line terms for gene families related to PD, which were downloaded with the PhytoMine interface of Phytozome 12 (Goodstein et al. 2012). Subsequent comparisons were made using the Basic Local Alignment Search Tool (BLAST) function with known amino acid sequences from related species. Synteny comparisons between common bean and soybean (Glycine max) were made using the Legume Information System 2.0 (Rice et al. 2015); these were verified using available literature (McClean et al. 2010, Schmutz et al. 2014). The CoGe SynMap (Lyons et al. 2008) and LegumeIP 2.0 (Li et al., 2016) synteny tools were used to compare syntenic regions between Arabidopsis (Col-0, TAIR10), common bean (G19833, Pvulgaris_V1.0_218; Schmutz et al. 2014), and soybean (Williams 82, Release 1.1; Schmutz et al., 2010). A neighbor-joining tree was produced to determine the pattern of homology between a common bean candidate gene (PvPdh1), a related soybean gene (GmPDH1), and other members of the dirigent gene family in these two species. The amino acid sequence of these proteins was BLASTed against the G. max and P. vulgaris proteomes identify closely related genes. These were then compared using a multiple BLASTP to develop a distance tree based on a Grishin protein distance matrix. A fast-minimum evolution tree (Desper & Gascuel 2004) was generated based on a maximum sequence difference of 0.85.

Amino acid conservation analyses
The complete amino acid sequence of PvPdh1 from accession G19833 was compared via BLASTP against the NCBI proteome database, using a BLOSUM62 matrix for comparison and existence and extension costs of 11 and 1, respectively (Altschul et al. 2005). The COnstraint-Based multiple ALignment Tool (COBALT; Papadopoulos & Agarwala, 2007, https://www.ncbi.nlm.nih.gov/tools/cobalt/re_cobalt.cgi) was used to align the most similar proteins known among several plant taxa and identify conserved residues based on the BLASTP results. The Protein Variation Effect Analyzer (PROVEAN; Choi & Chan, 2015) v1.1.3 software tool was used to estimate the effect of mutations of interest using default settings, including a cutoff threshold of -2.5 for identifying deleterious alleles.

Validation of the role of PvPdh1 in a wider population
Genomic DNA was extracted using a modified CTAB method; amplification and Sanger sequencing of PvPdh1 were conducted as described previously. Genotypes were separated into the Andean or Middle American gene pool based on an indel in the 3’ UTR of PvPdh1. This indel consistently predicted the gene pool in varieties of known ancestry. After sequencing, Middle American varieties were divided into groups based on amino acid at position 162 of PvPdh1. The degree of dehiscence between these groups was evaluated by Student’s t-test. Pod shatter phenotype data from the Germplasm Resource Information Network (GRIN) was compared with our sequencing data for varieties acquired from NPGS.

Landrace ecogeography 
Precipitation across the native range of Middle American beans was mapped in QGIS 2.18.19 using data from worldclim2 (Fick & Hijmans, 2017). National boundaries and coastlines were added using shapefiles available through Natural Earth (Kelso & Patterson, 2010). USGS topographical global raster data grids were also used to improve the visualization of coastlines (https://topotools.cr.usgs.gov/gmted_viewer/gmted2010_global_grids.php). Landraces genotyped by Kwak and Gepts (2009) were filtered by their ecogeographic race, those with values of 0.5 in STRUCTURE groups K6 (race Mesoamerica) and K9 (race Durango/Jalisco) were filtered for subsequent analysis. Delimited text layers were added in QGIS for varieties with latitude and longitude data that belonged to one of the ecogeographic races of interest. The average annual precipitation and elevation of the region where each landrace was collected using the “add raster values to points” function in QGIS, and the values between ecogeographic races were compared by student’s t-test.