ORIGINAL RESEARCH

published: 15 November 2017

doi: 10.3389/fpls.2017.01891

Domestication Genomics of the

Open-Pollinated Scarlet Runner
Bean (Phaseolus coccineus L.)
Azalea Guerra-García 1,2* , Marco Suárez-Atilano 3 , Alicia Mastretta-Yanes 4 ,
Alfonso Delgado-Salinas 5 and Daniel Piñero 2
1
Posgrado en Ciencias Biológicas, Universidad Nacional Autónoma de México, Ciudad de México, Mexico,
2
Departamento de Ecología Evolutiva, Instituto de Ecología, Universidad Nacional Autónoma de México, Ciudad de México,
Mexico, 3 Departamento de Ecología de la Biodiversidad, Instituto de Ecología, Universidad Nacional Autónoma de México,
Ciudad de México, Mexico, 4 CONACYT-CONABIO, Comisión Nacional para el Conocimiento y Uso de la Biodiversidad,
Ciudad de México, Mexico, 5 Departamento de Botánica, Instituto de Biología, Universidad Nacional Autónoma de México,
Ciudad de México, Mexico

The runner bean is a legume species from Mesoamerica closely related to common
bean (Phaseolus vulgaris). It is a perennial species, but it is usually cultivated in small-
scale agriculture as an annual crop for its dry seeds and edible immature pods. Unlike
the common bean, P. coccineus has received little attention from a genetic standpoint.
In this work we aim to (1) provide information about the domestication history and
Edited by:
domestication events of P. coccineus; (2) examine the distribution and level of genetic
Alejandro Casas, diversity in wild and cultivated Mexican populations of this species; and, (3) identify
Universidad Nacional Autónoma
candidate loci to natural and artificial selection. For this, we generated genotyping
de México, Mexico
by sequencing data (42,548 SNPs) from 242 individuals of P. coccineus and the
Reviewed by:
Peter J. Prentis, domesticated forms of the closely related species P. vulgaris (20) and P. dumosus (35).
Queensland University of Technology, Eight genetic clusters were detected, of which half corresponds to wild populations and
Australia
Gonzalo Gajardo,
the rest to domesticated plants. The cultivated populations conform a monophyletic
University of Los Lagos, Chile clade, suggesting that only one domestication event occurred in Mexico, and that it took
*Correspondence: place around populations of the Trans-Mexican Volcanic Belt. No difference between
Azalea Guerra-García
wild and domesticated levels of genetic diversity was detected and effective population
azalea.guerra@iecologia.unam.mx
sizes are relatively high, supporting a weak genetic bottleneck during domestication.
Specialty section: Most populations presented an excess of heterozygotes, probably due to inbreeding
This article was submitted to
depression. One population of P. coccineus subsp. striatus had the greatest excess
Evolutionary and Population Genetics,
a section of the journal and seems to be genetically isolated despite being geographically close to other wild
Frontiers in Plant Science populations. Contrasting with previous studies, we did not find evidence of recent gene
Received: 31 July 2017 flow between wild and cultivated populations. Based on outlier detection methods,
Accepted: 18 October 2017
Published: 15 November 2017
we identified 24 domestication-related SNPs, 13 related to cultivar diversification and
Citation:
eight under natural selection. Few of these SNPs fell within annotated loci, but the
Guerra-García A, Suárez-Atilano M, annotated domestication-related SNPs are highly expressed in flowers and pods. Our
Mastretta-Yanes A,
results contribute to the understanding of the domestication history of P. coccineus,
Delgado-Salinas A and Piñero D
(2017) Domestication Genomics and highlight how the genetic signatures of domestication can be substantially different
of the Open-Pollinated Scarlet Runner between closely related species.
Bean (Phaseolus coccineus L.).
Front. Plant Sci. 8:1891. Keywords: domestication, genotyping by sequencing, Phaseolus coccineus, adaptative variation, population
doi: 10.3389/fpls.2017.01891 genomics

Frontiers in Plant Science | www.frontiersin.org 1 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

INTRODUCTION included in these studies, and they focused on European

domesticated populations. Phylogenetic analyses including more
The scarlet runner bean (Phaseolus coccineus L.) is one of the five samples from the wide distribution of P. coccineus could bring
Phaseolus species that were domesticated in Mesoamerica, and it clues about the number of domestication events that took place
is the third-most economically important, after P. vulgaris L. and in this species. For example, if cultivars are grouped in one
P. lunatus L. The domestication process of this species continues monophyletic clade, it would suggest one domestication event.
today both in the Americas and Europe, where it was introduced Another interesting feature of P. coccineus domestication
by the Spaniards. One of its main characteristics is its ability to history is that similar levels of genetic variation have been
tolerate cooler climates than other Phaseolus and up to date it is reported in wild and cultivated populations (Escalante et al.,
an important food source for smallholders and indigenous groups 1994; Spataro et al., 2011; Rodriguez et al., 2013). This
in Mexico (Salinas, 1988). Despite the cultural value, economic contradicts the population genetics models that predicts a genetic
importance, and agronomic potential of P. coccineus, little is diversity reduction and increased divergence between wild and
known about its domestication history and the genetic variability domesticated forms due to demographic factors and selection
of its wild and cultivated forms. at target loci (Meyer and Purugganan, 2013). This pattern has
Wild P. coccineus are perennial climbing plants, occurring been described in crops like sunflower (∼30%; Renaut and
mostly at mid-high elevations (1,000–3,000 m.a.s.l.), from Rieseberg, 2015); soybean (∼30%; Li et al., 2013); maize (∼17%;
northern Mexico (Chihuahua) to Panama (Salinas, 1988). It Hufford et al., 2012); and in cultivated Agave species (from
has 11 pairs of chromosomes and an estimated genome size 21 to 66%; Eguiarte et al., 2013). Also, in the Mesoamerican
of 660 Mb (Plant DNA C-values database). Contrasting with common bean a ∼20% reduction in genetic variation has been
the autogamous common bean, the scarlet runner bean is an reported (Schmutz et al., 2014). However, the amount of genetic
open-pollinated species. The high morphological diversity of diversity that is lost along domestication depends on several
this species has been classified under two subspecies (Freytag factors, including the severity and the number of bottlenecks, the
and Debouck, 2002): P. coccineus subsp. coccineus (mostly with strength of selection and human management (Gepts, 2014). To
red flowers), including 11 wild varieties and the domesticated properly assess impact of domestication on the genetic diversity
form, and P. coccineus subsp. striatus (purple or mauve flowers), of P. coccineus, genomic data comparing wild and cultivated
conformed by eight wild varieties. No genetic evidence supports populations is necessary.
these subspecies and varieties, but given the environmental and The use of genomic tools also allows to characterize diversity
cultural heterogeneous landscape where P. coccineus occurs, it is and differentiation patterns across genomes. Regions or variants
expected that the species should be genetically structured. that departure from neutral predictions are probably influenced
As a cultivated species, P. coccineus is currently grown in by selective pressures and are tagged as candidates. Applying
Mexico, Guatemala, Honduras and Costa Rica, and in lesser this approach to crop species and their wild relatives allows
degree in South America. In Europe, it is mostly cultivated in the to distinguish loci affected during domestication, whereas
United Kingdom, Netherlands, Italy, and Spain (Rodiño et al., comparisons between landraces and/or improved cultivars
2006). In Mexico, the scarlet runner bean is cultivated both as measure the effect of subsequent selection (Tang et al.,
a self-sufficiency crop by smallholder farmers (<5 ha) and also 2010; Gepts, 2014). Furthermore, hypotheses about phenotypic
commercially for urban areas. Besides its native cultivars, in convergence in crops can be tested. In other words, if the same
Mexico there is one breeding line (Blanco Tlaxcala) developed genes or genomic regions were affected during the domestication
using a multi linear method (Vargas-Vázquez et al., 2012). Feral process of different species.
populations are common, but it is unknown if they originated Here, we aim to deal with the previous knowledge gaps
from hybridization between wild and domesticated populations, by using genomic data to (1) provide information about the
or if they escaped from cultivation. Wild, feral and domesticated domestication history of P. coccineus and its current evolutionary
distributions overlap in Mesoamerica, suggesting that there are dynamic in Mexico, in particular to analyze the occurrence of a
plenty of opportunities for gene flow to occur, making the single or multiple domestication events in Mexico; (2) examine
domestication history of P. coccineus difficult to disentangle the extent of the domestication bottleneck in this species by
without high resolution genetic markers. comparing the levels of genetic diversity and geographic patterns
The domestication history of the scarlet runner bean has of the wild, feral and domesticated Mexican populations; and (3)
been explored previously with low resolution molecular markers, identify candidate loci under natural and artificial selection in
and multiple domestication events were suggested. Specifically, P. coccineus genome.
chloroplast and nuclear SSRs of P. coccineus accessions including
European domesticated populations, Mesoamerican landraces
and wild samples from Mexico, Guatemala, and Honduras MATERIALS AND METHODS
(Angioi et al., 2009; Spataro et al., 2011; Rodriguez et al.,
2013) suggest that P. coccineus domestication took place in Plant Material and SNP Genotyping
the Guatemala-Honduras area, or that alternatively another Phaseolus coccineus individuals from 10 wild, three feral and
domestication event occurred in Mexico followed by extensive 11 cultivated Mexican populations and one cultivar from Spain
hybridization with the cultivated populations from Guatemala were analyzed, as well as plants from the breeding line Blanco
and Honduras. However, few Mesoamerican samples were Tlaxcala. Taxonomic and wild/feral/domesticated categories were

Frontiers in Plant Science | www.frontiersin.org 2 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

assigned based on morphology and habitat observations. Only genetic groups, both cultivated and wild, and if each cluster
one of the wild populations that were sampled corresponds to forms a monophyletic clade. This phylogenetic analysis was also
subsp. striatus, the rest belong to subsp. coccineus. A population used as a preliminary approach to identify the plausible number
was classified as feral if it was growing out of cultivation and of domestication events for the Mexican cultivated P. coccineus
presented intermediate traits between wild and domesticated (see below for other analyses). Specifically we examined if the
forms. The Mexican samples cover the species distribution and cultivated samples was recovered as a monophylogenetic group.
main cultivation areas at the national level. As outgroups, samples For the phylogenetic analysis, wild and cultivated samples of
from the closely related species P. vulgaris (three wild and one P. coccineus, P. vulgaris, and P. dumosus were analyzed under
cultivated) and P. dumosus (seven cultivated) were included three schemes:
(Supplementary Table S1). For the three species, the samples size First, a Maximum-Likelihood based approach was carried
of each population varied between three to 16 individuals. out with the FastTree software (Price et al., 2009). For
Sampling was performed during September–December of this, a mix of Nearest-Neighbor Interchange and Subtree-
2014 and 2015. In the case of the wild populations, tissue Prune and Regraft moves (NNI+SPR) was considered for
from young leaves was collected and stored in silica until topology and branch-length optimization and the General-Time
processed. Seeds from cultivars were collected and germinated Reversible with a single rate per site model (GTR+CAT) was
at the Instituto de Ecología, UNAM. DNA was extracted using included as nucleotide substitution model. Because FastTree
DNeasy Plant Mini Kit (Qiagen). DNA samples were genotyped only considers those SNPs identified as fixed within individuals
at the Institute for Genomic Diversity at Cornell University (i.e., homozygous), but polymorphic among individuals, only the
(Services | Institute of Biotechnology, 2017). Sequencing libraries 82% of the total VDC subset (41,223 SNPs) were considered
were constructed using enzymes PstI and BfaI following the in this analysis. Second, a phylogenetic network based on the
Genotype by Sequencing (GBS) protocol of Elshire et al. (2011). Neighbor-net algorithm and Patristic Distances with GTR+I+G
A total of 326 samples were processed in four plates of ninety correction was estimated with SplitsTree (Huson and Bryant,
six samples each, multiplexed and sequenced on four lanes of 2006) software. Lastly, we employed a Bayesian multispecies
Illumina HiSeq 2500 (100 bp, single-end reads). coalescent model (Rannala and Yang, 2003) to estimate the
Reads were aligned to P. vulgaris reference genome v1.0 phylogenetic relationships among well-supported clades within
(Phytozome) DOE-JGI and USDA-NIFA, http://phytozome.jgi. P. coccineus solely. We used the program SNAPP 1.3.0 (Bryant
doe.gov/ (Phytozome, 2017) using bwa v 0.7.8-r455; (Li and et al., 2012), included in the package BEAST 2.4.5 (Bouckaert
Durbin, 2009). Demultiplexing, initial quality control, assembly et al., 2014) to infer species trees directly from biallelic genetic
and SNP discovery were made with TASSEL pipeline v3.0.174 data. We used the eight main genetic clusters (see section Results)
(Glaubitz et al., 2014). Assembly and SNP discovery were inferred by Admixture as a priori designated species and the
performed independently for two sets of data, one containing Wild-TMVB cluster was partitioned in two, taking into account
samples from P. vulgaris, P. dumosus, and P. coccineus (VDC the ML topology of that cluster. Because SNAPP does not
group), which are the domesticated species of the Vulgaris incorporate missing data, we selected a subset of our taxonomic
clade (Delgado-Salinas et al., 2006); and the other data set only sampling that maximized the number of SNPs available. The final
including P. coccineus samples. SNPs were filtered in VCFtools analysis retained a total of 600 SNPs under linkage equilibrium;
0.1.15 (Danecek et al., 2011) using the following parameters without any missing data and considering a minimum of five
for the two data sets: (1) VDC group: maximum missingness individuals from each cluster of the designated species. We
threshold 20% per individual; minimum mean depth 10X; used SNAPP’s default settings and ran the analysis for 1,000,000
minimum allele frequency (MAF) 0.01; minimum allele count generations sampling every 1,000 generations. We evaluated the
90%; and only SNPs mapped in chromosomes. (2) Phaseolus convergence (i.e., short variation in -lnL scores, ESS > 100)
coccineus: maximum missingness threshold 30% per individual; from our runs by examining log files with the program Tracer
minimum mean depth 5X; MAF 0.02; minimum allele count 80%; 1.5 (Drummond and Rambaut, 2007). We analyzed the tree
and only SNPs mapped in chromosomes. files with SNAPP-TreeSetAnalyser 2.4.5, to identify species trees
Filtered SNP data, species occurrence data and scripts used for that were contained in the 95% highest posterior density (HPD)
the analyses are available at Dryad Repository under the identifier set and using 10% of topologies as burn-in. Resulted tree files
doi: 10.5061/dryad.q343c. (cloudgrams) were visualized using DensiTree (Bouckaert, 2010).

Inferring Population Structure and Population Genetics Statistics

Phylogenetic Relationships To evaluate the existence and degree of the domestication
We inferred the population structure of P. coccineus because bottleneck on P. coccineus we estimated genetic diversity and
different genetic clusters are expected to occur due to the isolation differentiation indices of the genetic groups inferred by the
and environmental and cultural heterogeneity in which this Admixture analysis (see section Results). Specifically, we used
species occurs. For this, the software Admixture v1.3 (Alexander the Hierfstat package (Goudet, 2005) in R (R Core Team, 2017)
et al., 2009) was used to infer population structure of P. coccineus. to estimate per site heterozygosity and F IS , as well as pairwise
Values of K ranging from one to twenty were tested, and the F ST among groups, performing a bootstrap (1,000) to obtain
value that exhibited the lowest cross-validation error was chosen. confidence intervals. To test the hypothesis that ni = nj (where
Then, we examined the phylogenetic relationships between the ni is the number of loci of the cluster i where H Ei > H Ej , and

Frontiers in Plant Science | www.frontiersin.org 3 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

n2 is the number of loci of the cluster j where H Ej > H Ei ) we principal components analysis (PCA), to ensure that at least
used a pairwise χ2 tests with Bonferroni correction to avoid false one (or more) scenarios would produce simulated datasets close
positive results (Sokal and Rohlf, 1995). Also, we estimated the enough to the empirical data. The PCA was based on a set of
heterozygosity and F IS at the sampling location (P. coccineus 5,000 simulated datasets, generated from the parameters’ prior
dataset) and at the species level (VDC dataset) applying the same distributions (Supplementary Figure S3).
test.
Identifying Candidate Loci
Multiple vs. Single Domestication Events We used the wild and cultivated samples of P. coccineus to
Test identify candidate loci related to domestication, to cultivar
In order to confirm the hypothesis of a single domestication diversification, and to natural selection. Before the candidate
event in Mexico suggested by our phylogenetic analyses SNPs analysis, an additional filter based on linkage disequilibrium
(see section Results) we applied the Approximate Bayesian (LD) was applied. To determine the threshold distance at which
computation (ABC; Beaumont et al., 2002) method implemented there is no LD, we estimate the inter-variant allele correlations
in DIYABC 2.04 (Cornuet et al., 2014). Preliminary tests included (r2 ) using PLINK 1.9 (Chang et al., 2015). To distinguish LD due
comparisons among three scenarios with 3 × 106 simulated to physical distance (bp), the r2 was estimated for SNPs located in
datasets (1 × 106 each scenario) in which the position of the same and in different chromosomes. The distance threshold
the Wild-Sierra Madre Occidental (Wild-SMOCC) clade was was established in 3,000 bp, so that SNPs closer than this distance
evaluated (see section Results, Supplementary Figure S1). Our were removed.
final estimation included 4 × 106 simulated datasets (2 × 106 This LD-filtered dataset was analyzed with two different
each scenario) considering the Wild-SMOCC population fixed approaches for outlier detection: the R package pcadapt (Luu
as sister clade of the Wild-Trans-Mexican Volcanic Belt (Wild- et al., 2017) and BayeScan 2.1 (Foll and Gaggiotti, 2008). Only
TMVB) populations (see section Results). The number of loci identified by the pcadapt and BayeScan methods were
domestication events was tested as follows: multiple events considered as candidate loci. Pcadapt detects candidate SNPs
(Scenario 1, Supplementary Figure S2) vs. a single one (Scenario assuming that these are outliers with respect to how they are
2, Supplementary Figure S2). The DIYABC approach was also related to population structure. By contrast to population-based
applied to estimate the time at which domestication occurred, approaches, pcadapt does not require grouping individuals into
as well as other demographic parameters such as effective populations and handles admixed individuals (Luu et al., 2017).
population size (Ne). A subsample from the SNAPP dataset BayeScan instead uses differences in allele frequencies of pre-
(279 SNPs) and the scheme of eight clusters were used to set defined populations, in this case the genetic clusters previously
populations in DIYABC (Figure 2B). Priors were set as follow: established by Admixture.
log-uniform distributions across all parameters, Ne ranging from In both approaches, three separate analyses were performed
100 to 100,000 individuals, mutation rate set to 10−8 –10−6 across with each method to detect signatures of different types of
SNPs, and divergence times among populations set to 10–100,000 selective pressure. First, to detect candidate domestication loci,
generations ago (Table 1). wild and cultivated samples of P. coccineus were included, and
We compared the fit of the single vs. multiple domestication feral individuals were removed. In this case, for the pcadapt
events scenarios by estimating their posterior probabilities: with analysis, only the first principal component was assessed because
the obtained reference tables from each scenario, we ranked the it explains the difference between wild and cultivated populations
simulated datasets in order of increasing distance to the observed (see section Results). Also, an additional SNPs filter was made
data considering direct and logistic approaches (Beaumont et al., and MAF were adjusted to consider SNPs present in at least five
2002; Cornuet et al., 2014). Distance between datasets was individuals. For this dataset, that is MAF = 0.023. For Bayescan
based on summary statistics, estimated from the empirical and no additional filter was made. Second, to identify loci related to
simulated sets. We performed a pre-evaluation step using a diversification in the context of domestication, only cultivated

TABLE 1 | Estimations of effective population sizes of the best-fit DIYABC model (single domestication) for Phaseolus coccineus in Mexico.

Genetic group Minimum prior value Maximum prior value Average posterior value 95%CI

Wild-SUR-CH 100 1 × 105 1.0 × 105 9.98 × 104 –1 × 105

Wild-SMOCC 100 1× 105 8.94 × 104 8.01 × 104 –9.5 × 104
Wild-TMVB 100 1 × 105 8.94 × 104 8.01 × 104 –9.5 × 104
Wild-striatus 100 1 × 105 9.68 × 104 8.9 × 104 –1 × 105
Cult-OV 100 1 × 105 8.83 × 104 7.5 × 104 –1 × 104
Cult-SUR-CH 100 1× 105 6.39 × 104 5.7 × 104 –9.38 × 104
Cult-TMVB 100 1× 105 9.36 × 103 8.39 × 103 –1.54 × 104
Cult-SMOCC 100 1× 105 8.57 × 103 6.01 × 103 –1.2 × 104

Please refer to the text to understand what the acronyms stand for.

Frontiers in Plant Science | www.frontiersin.org 4 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

FIGURE 1 | Distribution map of genotyped populations of Phaseolus coccineus in Mexico, assigned to eight genetic clusters with the program Admixture (colors).
Circles indicate wild populations and triangles show cultivated populations (Figure 2). Frames indicate the seven first-level ecoregions (9–15) present in Mexico (level
1) and boundaries represent 21 second-level nested ecoregions. Shaded area shown the potential geographical distribution of P. coccineus in Mexico (López-Soto
et al., 2005).

samples were analyzed. In the pcadapt analyses, the first six RESULTS
components were assessed because they explain the genetic
structure of populations, and MAF threshold was set to 0.038 to Sampling and SNP Genotyping
excluded alleles present in less than five individuals. Notice that A total of 296 individuals representing four ecoregions of Mexico
in this case, diversification refers to the phase that follows initial (as defined in Instituto Nacional de Estadística, Geografía e
domestication and involves the spread and adaptation to different Informática (INEGI), Comisión Nacional para el Conocimiento
agro-ecological and socio-cultural environments (Meyer and y Uso de la Biodiversidad (CONABIO), and Instituto Nacional de
Purugganan, 2013). Lastly, to detect natural selection signatures, Ecología (INE), 2008) were sampled and successfully genotyped
we focused both methods on wild samples. Again, for the pcadapt (Figure 1). After assembly and SNP discovery, the VDC
analyses the first six components were assessed and the MAF group dataset contains 241 individuals of P. coccineus, 20
threshold was set 0.055 to exclude SNPs present in less than of P. vulgaris and 35 of P. dumosus, 50 273 SNPs, 2.24%
five individuals. In all cases, no additional filter was made for mean missing data per individual, and a mean depth per site
BayeScan. of 58.63. The P. coccineus dataset includes 242 individuals
The false discovery rate threshold applied in pcadapt and (91 wild; 20 feral; 131 cultivated), 42,548 SNPs, 3.97% mean
BayeScan were 0.005 and 0.05, respectively. To compare how missing data per individual, and a mean depth per site of
genetic variance is explained by candidate SNPs and by data set 50.41.
LD filtered, PCAs were made using the SNPrelate package (Zheng
et al., 2012). Inferring Population Structure and
Using Phytozome’s JBrowser, the putative function and tissue
of expression of these loci was examined by looking for the Phylogenetic Relationships
annotation of the selected SNPs in P. vulgaris genome v 2.1 The K-value that presents the lower error rate in Admixture
(DOE-JGI and USDA-NIFA1 ). For each annotated loci we looked analysis was eight (Supplementary Figure S4). Half of the
for homologous proteins with the highest similarity in other genetic groups correspond to the cultivars from the Trans-
plants, and examined if the homolog genes in Glycine max Mexican Volcanic Belt (Cult-TMVB), Sierra Madre del Sur and
(soybean) were among the domestication-related loci associated Chiapas Highlands (Cult-SUR-CH), Sierra Madre Occidental
with flowering time and seed size in this species (Zhou et al., (Cult-SMOCC) and Oaxaca Valley (Cult-OV). The other half
2015). of the genetic clusters belong to wild populations from the
Trans-Mexican Volcanic Belt (Wild-TMVB), Sierra Madre del
1
http://phytozome.jgi.doe.gov/ Sur and Chiapas Highlands (Wild-SUR-CH), Sierra Madre

Frontiers in Plant Science | www.frontiersin.org 5 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

FIGURE 2 | Overall phylogenetic relationship among 242 individuals of P. coccineus from Mexico. Numbers and colors represents the eight main genetic clusters
solved by Admixture: (1) Wild-SUR-CH, (2) Wild-TMVB, (3) Wild-striatus, (4) Wild-SMOCC, (5) Cult-SUR-CH, (6) Cult-OV, (7) Cult-TMVB and (8) Cult-SMOCC. Please
refer to the text to understand what the acronyms stand for. (A) Maximum-Likelihood tree, main doted clades indicated support bootstrap values >75%. Black dot
next to cluster 8 indicates the Blanco Tlaxcala breeding line; dark triangle shows the Spanish population; feral samples are indicated with an asterisk. (B) Individual
assignment based on 42,548 SNP’s solved with Admixture. (C) Rooted Neighbor-Net topology achieved by SplitsTree.

Occidental (Wild-SMOCC) and subsp. striatus population, took place, although the Wild-SUR-CH genetic cluster can be
located in the TMVB (Wild-striatus; Figure 2). The genetic discarded.
clusters seem to be related to geographic distances (Figure 1), The ML and Neighbor-Net topologies in which P. dumosus
except the population Wild-striatus, which is geographically and P. vulgaris were included, positioned P. dumosus as a sister
close to populations of P. coccineus subsp. coccineus but group of P. coccineus (Figure 2A). However, the SplitsTree
seems genetically isolated. Samples from the Spanish population method indicated a basal reticulate pattern among P. dumosus,
(Figure 2B, triangle) were assigned to the Cult-TMVB genetic P. coccineus, and P. vulgaris (Figure 2C), suggesting ancestral
group, but unlike the individuals of this cluster, samples gene flow, but not recent. Furthermore, there is no evidence of
from Spain do not present a mixed ancestry. Regarding recent gene flow between wild and cultivated groups, but only
samples of the breeding line Blanco Tlaxcala (Figure 2B, within genetic clusters (Figure 2C).
circle), they are grouped with landraces from Cult-SMOCC Regarding SNAPP cloudgram (Figures 3B,C), 53 single
cluster. topologies summarize the 95% HPD consensus tree, indicating
The phylogenetic hypotheses constructed with FastTree and a different divergence pattern in which Wild-TMVB populations
SplitsTree (Figures 2A,C) are consistent with the Admixture are the closest clade to the domesticated group. Nevertheless,
genetic groups (Figure 2B). Nevertheless, both analysis suggested the complex assignment of individuals within Wild-TMVB
the Wild-TMVB group as a paraphyletic clade. ML topology and Wild-striatus are shown in a non-solved pattern within
revealed a finer-scale structure, identifying three paraphyletic the cloudgram as well as in low values of nodal support
clades within this genetic cluster, and Wild-striatus cluster is in the consensus topology (Figure 3C). Despite these main
a nested clade differentiated from the rest of the Wild-TMVB inconsistencies between ML and Neighbor-Net vs. SNAPP
group (Figure 2). Remarkably, the domesticated populations topologies, all hypotheses favor the occurrence of a single
integrate a monophyletic clade statistically well supported, domestication event.
suggesting a unique domestication event for the Mexican In regards of the ABC-based computations, the model
populations. Nevertheless, these phylogenetic hypotheses do not comparisons in preliminary trials indicated scenarios where the
allow to distinguish the genetic pool from which domestication Wild-SMOCC population that are paraphyletic to Wild-TMVB

Frontiers in Plant Science | www.frontiersin.org 6 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

FIGURE 3 | (A) Cloudgram depicting topologies of 9,999 species trees obtained from an analysis of 600 single nucleotide polymorphism loci from 124 P. coccineus
using SNAPP; (B) average assignment probability achieved by Admixture of selected individuals considered in species tree analyses based on nine groups; (C)
associated root canal depicting a consensus topology from SNAPP analysis. Nodal support values on the root canal are posterior probabilities that correspond to
strongly supported nodes designated a priori in the species tree analysis.

yielded a higher probability in both direct and logistic approaches P. coccineus showed the highest diversity and P. dumosus the
(Supplementary Figures S1, S2). A final test indicated that the lowest.
most likely scenario was a single domestication event, being Outstandingly, H O was greater than H E in all the genetic
the Wild-TMVB group the closest to the domesticated clade groups except in the Wild-SUR-CH cluster, resulting in
(Figure 4; Scenario 2, direct P = 0.786, logistic P = 1.0), which negative values of F IS . Within the groups with an excess of
is congruent with the results of SNAPP phylogenetic analyses. observed heterozygosity, Wild-striatus had the lowest inbreeding
Evaluation of the posterior predictions via PCA indicated that coefficient (Figure 5). On the contrary, at the species level
parameter values and summary statistics from the simulated P. vulgaris showed a deficit of heterozygotes, showing a high
datasets based on Scenario 1 closely matched the empirical data F IS . The inbreeding coefficient is positive when estimated taking
(Supplementary Figure S3). into account all P. coccineus samples. This is caused by the
Wahlund effect, which is the reduction of heterozygosity due
Wild and Domesticated Population to subpopulation structure. Regarding pairwise differentiation
Genetics Statistics index, F ST values ranged from 0.022 (Cult-TMVB vs. Cult-
High levels of genetic diversity were found in wild and cultivated SMOCC) to 0.178 (Cult-OV vs. Wild-striatus; Figure 6). As
populations (Figure 5). At the genetic cluster level, the Wild- expected, the pair F ST values are greater between wild genetic
TMVB group presented the highest diversity and the Cult- groups than between cultivated genetic clusters (Figure 6).
OV group the lowest. No clear pattern in the amount of Cultivated populations of P. coccineus show smaller effective
diversity was observed between wild and cultivated clusters. population sizes than wild populations. In some cases, like in
There were cultivated groups with high genetic variance (Cult- Cult-TMVB and Cult-SMOCC, Ne was one order of magnitude
SUR-CH and Cult-TMVB), and wild clusters that presented lower smaller than in the rest of the populations. On the contrary, the
diversity than cultivated populations (Wild-SMOCC). At the genetic cluster Wild-SUR-CH had the biggest Ne (Table 1). The
location level (Supplementary Table S2), the samples from Spain most recent split was estimated to happen 3.9 × 103 generations
(H E = 0.134) and Oaxaca Valley (H E = 0.148) presented the ago, and occurred between the Cult-SMOCC and the Cult-TMVB
lowest diversity, and the highest was found in wild population clusters. On the contrary, the oldest split event was dated in
located in Tlalpan, Mexico City (H E = 0.208). Regarding species, 4.95 × 105 generations ago between the Wild-SUR-CH and the

Frontiers in Plant Science | www.frontiersin.org 7 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

FIGURE 4 | Best-fitted domestication scenario of P. coccineus achieved with DIYABC. Split times in generations (tn) indicated the average posterior value estimated
after Bayesian Computations (95% CI).

FIGURE 5 | HE in green (variance interval), HO in red (variance interval) and inbreeding coefficient in blue (95% IC) estimated for genetic clusters of P. coccineus and
for species of VDC group. Letters show groups that are statistically different and are decreasingly ordered according HE .

rest of P. coccineus clade. The split event that separates wild and Identifying Candidate Loci
domesticated samples was dated about 2.1 × 104 generations ago Before LD filtering, the mean r2 value among SNPs located in
(Figure 4). Since P. coccineus is usually treated as an annual when the same chromosome separated by a maximum distance of
cultivated, that represents 21,000 years ago. In the case of wild, 10,000 bp was 0.151. After eliminating SNPs closer than 3,000 bp,
perennial plants, one generation could be more than a year. the mean r2 was 0.063 (Supplementary Figure S5). In the case of

Frontiers in Plant Science | www.frontiersin.org 8 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

FIGURE 6 | Heatmap representing the pair F ST values between genetic clusters.

SNPs from different chromosomes, the mean r2 was 0.022. This genetic and geographic structure of wild and cultivated groups
low LD is not due to the closeness, but rather by factors like can be recovered by these few candidate SNPs (Figures 7B,C)
populations structure. Interestingly, the pattern in the decay of and a clear separation of wild and domesticated populations is
LD differed between genetic groups, with the fastest decay and observed (Figure 7A).
lowest r2 in cultivated and wild populations from the TMVB. Four SNPs of the candidate domestication loci were found
Meanwhile, Wild-striatus, Wild-SURCH and Cult-OV had the to be annotated in P. vulgaris genome, one of the candidate
slowest LD decay and highest r2 values (Supplementary Figure loci under natural selection and none of the candidate loci
S5). After filtering, the data set for candidate loci contained for cultivar diversification (Supplementary Table S3). Three of
11,693 SNPs distributed across the 11 chromosomes. In the the annotated candidate domestication loci (Phvul.001G232200,
central region of most of the chromosomes, there is a reduction Phvul.007G256000, Phvul.009G156400) are highly expressed in
in SNP density, probably due to centromeres (Supplementary flowers, flower buds or young pods, and the remaining locus
Figure S6). (Phvul.002G145600) is highly expressed in green mature pods.
Using the pcadapt package, 47 SNPs were identified All these loci have their highest similarity homologs in G. max
as candidate domestication loci; 342 involved in cultivar genome v2.0 (Schmutz et al., 2010), but none of these correspond
diversification; and 1,030 potentially under natural selection. to the domestication-related loci previously identified by Zhou
Despite the great number of candidate SNPs that were identified, et al. (2015). The annotated candidate locus for natural selection
few are shared among selection types (Supplementary Figure (Phvul.003G197500) is highly expressed in roots and steam and
S7). In the case of the BayeScan analyses, 469 candidate corresponds to a calmodulin binding protein-like, which also has
SNPs for domestication were identified; 16 related to cultivar an homolog in G. max.
diversification; and 12 candidates associated with natural
selection. None of these SNPs were shared among the three
BayeScan analysis. DISCUSSION
Twenty four SNPs related to domestication, 13 to cultivar
diversification and eight to natural selection were detected A Single Domestication Event for
by both approaches and considered as candidate loci for Mexican P. coccineus in the TMVB
further analyses (Supplementary Table S3). The genetic variance Spataro et al. (2011) and Rodriguez et al. (2013), using SSR
explained by the candidate SNPs compared to the 11,693 SNPs data, suggested two domestications events of P. coccineus, one
used previously changed dramatically (Figure 7). Notably, the in Mexico and the other in Guatemala-Honduras. The genomic

Frontiers in Plant Science | www.frontiersin.org 9 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

FIGURE 7 | Principal components analysis plot for the first two principal components using LD dataset (left) and candidate SNPs identify in both in PCAdapt and
BayeScan analysis (right). (A) Analysis including cultivated and wild samples (no feral) to detect domestication-related SNPs. (B) Analysis of cultivated populations to
distinguish cultivar diversification-related loci. (C) Analysis of wild samples to detect signatures of natural selection.

data generated in this work indicates a unique domestication one generation per year in cultivated populations, this represents
event for the cultivated populations from Mexico (Figures 2, 3). 3,950 years. But divergence between the cultivated and wild
This includes Chiapas populations (Cult-SUR-CH), which are clades was dated in 21,000 generations (Figure 4, t5). This date
geographically and culturally closer to Guatemala than to Central is out of range of any plant domestication event and it seems
and Northern Mexico. However, no samples from Guatemala unlikely. There are evolutionary processes that may affect these
and Honduras were included, therefore a second domestication estimations. Processes like selection, population subdivision and
event in this area cannot be discarded with the present data. incomplete lineage sorting may result in an overestimations of
Nonetheless, based on the results from SNAPP and DIYABC divergence times because increase the time to coalescence, that
analyses, we were able to identify Wild-TMVB as the genetic pool is, the time it takes for the two sequences to find their common
from which domestication started in Mexico (Figure 3). ancestor (Albrechtsen et al., 2010; Angelis and Dos Reis, 2015). In
The most recent divergence time, that corresponds to the P. coccineus, the selection made by humans during domestication
separation between cultivated groups of SMOCC and TMVB, and the high population structure in wild and domesticated
was dated in 3,950 generations ago (Figure 4, t1). Assuming groups probably has resulted in overestimated divergence times.

Frontiers in Plant Science | www.frontiersin.org 10 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

Also, it has to be considered that wild populations are perennial The breeding line Blanco Tlaxcala grouped with SMOCC
and thus generation times may be longer than a year. landraces. Probably, breeding practices have acted over specific
The genetic findings suggest that P. coccineus domestication regions rather than over all the genome. The individuals of
likely occurred from TMVB’s material, pinpointing the this breeding line did not present mixed ancestry, despite
domestication of this species to a particular region within that Blanco Tlaxcala was developed using a multi linear
the large Mexican territory where it is cultivated nowadays. method (Vargas-Vázquez et al., 2012). This suggests that
Other sources of information could be incorporated to confirm all lines used to generate Blanco Tlaxcala belonged to the
this, using our findings as a geographic reference. If confirmed, same genetic cluster (Cult-SMOCC), and they were submitted
identifying the TMVB as the area where domestication started for to several rounds of strong selection, decreasing genetic
this species is interesting and important from an evolutionary, variation.
cultural and conservation perspective. The TMVB is the most Contrary to what was reported by Spataro et al. (2011)
recent mountainous region of Mexico, a biodiversity hotspot and and Rodriguez et al. (2013), samples from Spain clustered
it has a complex bio- and phylogeographic history characterized within the TMVB landraces, indicating that this European
by following a sky-island dynamic during the last 2 Myr population was originated by the introduction of individuals of
(Mastretta-Yanes et al., 2015). Culturally it became prominent the Cult-TMVB group into Spain. Nevertheless, because just one
during the Mexica Empire, and has been the most populated part European population was analyzed, no general pattern can yet
of Mexico since little before the Spanish conquest (Bataillon, be inferred. Notably, Spanish samples did not present mixed
1972). This has derived in other important cases of domestication ancestry, meanwhile the rest of the individuals of this genetic
to occur in this region. For instance, this central region was group did (Figure 2B). Probably the genetic bottleneck that
where the introgression of Zea mays ssp. parviglumis and Z. mays originated European populations and the isolation from wild
ssp. mexicana occurred during the domestication of maize (van relatives and American landraces, have decreased the amount of
Heerwaarden et al., 2011). However, human occupation in this shared ancestral polymorphisms between cultivars from TMVB
area is also a concern for conservation, because the growth of and Spain.
urbanization and high-input agriculture in this area threat both It has been suggested that hybridization and introgression
P. coccineus landraces and wild populations (CONABIO and have played a major role in P. coccineus evolution, both
IUCN, 2016). in cultivated and wild populations (Escalante et al., 1994;
Besides genetic data, a Mexican domestication origin of Angioi et al., 2009; Spataro et al., 2011; Rodriguez et al.,
P. coccineus is also supported by the several names that this bean 2013). Our results showed mixed ancestry both in wild and
has among different cultures. For instance, it is called tekómari cultivated clusters. However, little evidence of introgression and
in Chihuahua (Tarahumara indigenous language); tasukhu in hybridization was detected, and mixed ancestry can also be
Hidalgo and Puebla (Otomi); ayocote in central states of Mexico due to shared ancestral polymorphisms. Nevertheless, wild and
(Nahuatl); shaushana or xaxana in Veracruz (Totonaco); ma-má- cultivated populations frequently coexist, therefore hybridization
ja (Mazateco) in Oaxaca; and botil or shbotil chenec in Chiapas cannot be discarded and a formal test considering the number
(Tzeltal) (Salinas, 1988). Associated to these groups, there is also and size of introgressed regions and the direction of gene flow
considerable traditional knowledge regarding the cultivation and must be done.
use of P. coccineus species (e.g., Monroy and Quezada-Martínez,
2010). Phaseolus coccineus Is Highly Diverse
and Structured
Phaseolus coccineus wild populations are divided in four genetic
Historic and Recent Gene Flow among clusters that show considerable population differentiation.
Wild, Feral and Domesticated Similar levels of differentiation have been observed in several
Populations other highland species, which has been related to the high
The individuals identified as feral clustered in the domesticated environmental variability and the complex geologic and climatic
clade (Figure 2A), suggesting that they are escaped cultivars. history of Mexico (Mastretta-Yanes et al., 2015). The extent of this
This questions the hypothesis of an hybrid origin between differentiation in crop wild relative species has been mostly done
wild and cultivated populations (Salinas, 1988) and contrasts with low resolution neutral makers (Bellon et al., 2009; Piñero
with previous studies of feral P. vulgaris populations in Mexico et al., 2009) so it still needs to be further explored with genomic
(Papa and Gepts, 2003), where weedy populations appear to data. However, the present study and analyses in teosinte (van
be genetically intermediate between domesticated and wild Heerwaarden et al., 2011; Aguirre-Liguori et al., 2017), highlight
populations, and not cultivar escapees. Interestingly, the three that there is high diversity contained in the genetic pools of crop
collected feral populations belonged to the same genetic cluster wild relatives from Mexico.
(Cult-SUR-CH) and presented high levels of mixed ancestry, Besides the diversity contained in wild relatives, one of the
of which only a small proportion corresponds to wild clusters most important determinants in crop evolution is the level of
(Figure 2B). Since little evidence of gene flow was found in genetic diversity contained in the domesticated populations,
SplitsTree (Figure 2C), probably the mixed ancestry is due to especially with reference to the wild ancestral gene pool.
shared polymorphisms or ancestral gene flow, rather than recent Genetic diversity reduction has been widely described in crop
introgression events. domestication (Hufford et al., 2013; Li et al., 2013; Schmutz

Frontiers in Plant Science | www.frontiersin.org 11 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

et al., 2014; Renaut and Rieseberg, 2015). This reduction of Adaptative Variation in Wild and
genetic diversity is caused by genetic drift resulting from Domesticated Populations
population bottlenecks, and by artificial selection (Gepts, 2014).
Mexico is an environmentally and culturally heterogeneous
This phenomenon was also described in P. vulgaris (Schmutz
country, which favored crop genetic diversity. The distribution
et al., 2014) but in P. coccineus no clear pattern of genetic
of Phaseolus, both cultivated and wild, involves an interaction
reduction was found between the wild or cultivated genetic
with a wide range of different cultures, and isolated populations
groups (Figure 5). The Wild-TMVB cluster presented the highest
are exposed to diverse environmental conditions. For example,
genetic variation, followed by the Cult-TMVB and Cult-SUR-
compared to P. vulgaris, P. coccineus grows in more humid
CH groups. On the contrary, the Cult-OV and Wild-SMOCC
environments, at cooler temperatures and at higher altitudes.
clusters showed the lowest H E . Regarding effective population
Nevertheless, there are few studies that aim to elucidate the
sizes, these were greater in wild than in cultivated genetic
genetic basis of adaptation, especially for the wild populations
clusters, which is expected due to the genetic bottlenecks
of Phaseolus crop species (Bitocchi et al., 2017). Our outlier
associated to domestication process. Nevertheless, Ne estimations
analyses listed some candidate SNPs that could be under artificial
of domesticated groups are in the order of 103 –104 . Taking
selection during the domestication and diversification stages, and
together all results, these suggest that the genetic bottleneck
others that could be under natural selection. Although most of
during domestication was not severe. Other factors that may
these outliers are still not annotated, they could serve as a base
favor the maintenance of genetic diversity in P. coccineus
for identifying population differentiation in adaptive variation,
are its high outcrossing rate (Escalante et al., 1994) and
which is a needed step for genetic resources and crop wild
the fact that the genetic cluster from which domestication
relatives conservation (Maxted et al., 2012). Our study is based
started (Wild-TMVB) presents the highest diversity. Little
on GBS data, so P. coccineus genome is not fully saturated, and
evidence of recent gene flow was detected, but early gene
likely there are loci under selection that we did not sample.
flow could also favor the amount of genetic diversity in
Nevertheless, this set of outliers are a first approximation to
cultivars.
identify candidate loci to domestication and natural selection in
Analyzing the genetic variance at the location level, Spanish
runner bean.
samples presented the lowest diversity (Supplementary Table S2),
The fact that no loci overlapped between domestication,
which may be due to the recent demographic bottleneck that
diversification and natural selection categories shows that
occurred during its introduction to Europe. Nevertheless, Oaxaca
different selective processes were detected. This is to be
Valley also showed low genetic variation (Supplementary Table
expected because, in general, loci under natural selection and
S2) and the ancestry analysis (Figure 2B) suggests that it has been
artificial selection related to domestication and diversification are
genetically isolated from the other genetic clusters.
expected to differ across the genome (Meyer and Purugganan,
Regarding the inbreeding coefficient, the wild and cultivated
2013).
genetic clusters presented negative F IS values, indicating an
The loci involved in domestication are expected to be specially
excess of heterozygotes, except in the Wild-SUR-CH group.
related to the phenotypic changes of the domestication syndrome
A possible explanation for this pattern is inbreeding depression,
(Koinange et al., 1996), that is modifications in morphological
which effect in progeny has been studied in cultivars from
and physiological traits like seed dispersal, seed dormancy,
Spain, finding that selfing affected germination, survival rate
gigantism, increased harvest index and flowering time (Hammer,
and seed weight (González et al., 2014). Also, a negative
1984). Most of the domestication-related loci identified here
correlation was found between outcrossing rate and seed
are still of unknown function, but the four that are annotated
abortion in wild populations studied by Escalante et al.
are highly expressed in flowers or pods (Supplementary Table
(1994). In the case of domesticated populations, the bottlenecks
S3). This is interesting because in the soybean, another legume,
that they suffered during domestication may promote the
several domestication-related loci associated with flowering time
accumulation of deleterious alleles and the increase of inbreeding
have been identified (Zhou et al., 2015). However, no overlap
depression, resulting in lower values of the inbreeding coefficient
among those loci and the ones identified here was found.
(Morrell et al., 2011). Opposite to what was expected, in
P. coccineus the population with the lowest F IS was a wild
cluster (Wild-striatus). This population was previously studied CONCLUSION
by Búrquez and Sarukhán (1984), who found evidence of
self-incompatibility, which is congruent with our results. The SNPs generated in this work provided high resolution
A possible explanation for this pattern is the accumulation data to understand the domestication of P. coccineus. Results
of deleterious alleles in the Wild-striatus cluster. Notably, no suggest one domestication event for Mexico, which started from
mixed ancestry was detected in this genetic group, indicating the wild genetic pool from TMVB. Furthermore, wild and
that it is genetically isolated from other populations despite domesticated populations are highly diverse and presented high
being geographically close to other wild and cultivated TMVB values of Ne, suggesting that the demographic bottleneck due
populations. It is necessary to evaluate other populations of to domestication was not severe. These genomic analyses allow
P. coccineus subsp. striatus to know if this is a common to highlight how the genetic signatures of domestication can be
pattern and to explore the ecological and genetic causes and substantially different even between species of the same genus
consequences of it. domesticated in the same geographic area. Common bean and

Frontiers in Plant Science | www.frontiersin.org 12 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

scarlet runner bean are closely related species, nevertheless their the analyses. All authors revised the results and wrote the
reproductive strategies and domestication histories seem to be manuscript.
different: P. vulgaris tends to self-crossing, which theoretically
facilities the domestication process, and it also suffered a severe
domestication bottleneck. On the contrary, P. coccineus is an FUNDING
open pollinated species that presents high levels of genetic
diversity and population structure, and its domestication did not This work was supported by Consejo Nacional de Ciencia y
result in a strong demographic bottleneck. Tecnología through the Ph.D. scholarship number 440709 to
Our findings also show that both wild and domesticated AG-G and CONACYT Grant 247730 to DP.
populations of P. coccineus are highly structured. Most of the
genetic clusters presented an heterozygotes excess, showing
evidence of inbreeding depression. Interestingly, the population ACKNOWLEDGMENTS
identified as P. coccineus subsp. striatus shows the greatest
excess of heterozygotes and seems to be genetically isolated from We thank Idalia Rojas, Myriam Campos, Erick García,
other wild and cultivated populations. Contrasting with previous Verónica González, Alfredo Villarruel, Nancy Gálvez, and
studies, our data shows that gene flow within and between wild Rocío González for fieldwork assistance, Tania Garrido for
and cultivated populations is not a common process. Fully testing laboratory technical assistance and Ernesto Campos Murillo
this represents an area where further research is needed. for bioinformatic assistance to execute analyses in a cluster
The levels of diversity and population differentiation found environment. We acknowledge funding from the CONACYT
here support that the runner bean is a potential source of grant number 247730 and IEUNAM to DP. Statistical
variability for several traits for plant breeding (Schwember et al., analyses were carried out in the CONABIO’s computing
2017). The data presented here highlights that for a better cluster, which was partially funded by Secretaría de Medio
characterization of P. coccineus wild and cultivated forms there Ambiente y Recursos Naturales (SEMARNAT) through the grant
is still a need of more sampling, specially including Central “Contribución de la Biodiversidad para el Cambio Climático”
American populations. Complete and annotated genomes of to CONABIO. This work constitutes a partial fulfillment
Phaseolus and other legume crops will facilitate not only of the Posgrado en Ciencias Biológicas at the Universidad
comparative genomics, but will give a better knowledge of the Nacional Autónoma de México (UNAM) for AG-G. Finally,
evolution and domestication of this group of plants that has been we thank to all farmers that share with us their seeds and
independently domesticated by several human groups across its knowledge.
distribution.
SUPPLEMENTARY MATERIAL
AUTHOR CONTRIBUTIONS
The Supplementary Material for this article can be found online
AG-G, DP, and AD-S designed the study. AG-G made the at: https://www.frontiersin.org/articles/10.3389/fpls.2017.01891/
molecular procedures. AG-G, AM-Y, and MS-A conducted full#supplementary-material

REFERENCES Beaumont, M. A., Zhang, W., and Balding, D. J. (2002). Approximate Bayesian
computation in population genetics. Genetics 162, 2025–2035.
Aguirre-Liguori, J. A., Tenaillon, M. I., Vázquez-Lobo, A., Gaut, B. S., Jaramillo- Bellon, M. R., Barrientos-Priego, A. F., Colunga-García, M. P., Perales, H., Reyes-
Correa, J. P., Montes-Hernandez, S., et al. (2017). Connecting genomic patterns Agüero, J. A., Rosales-Serna, R., et al. (2009). Diversidad y conservación
of local adaptation and niche suitability in teosintes. Mol. Ecol. 26, 4226–4240. de recursos genéticos en plantas cultivadas. Capital Nat. México 2,
doi: 10.1111/mec.14203 355–382.
Albrechtsen, A., Nielsen, F. C., and Nielsen, R. (2010). Ascertainment biases in SNP Bitocchi, E., Rau, D., Bellucci, E., Rodriguez, M., Murgia, M. L., Gioia, T., et al.
chips affect measures of population divergence. Mol. Biol. Evol. 27, 2534–2547. (2017). Beans (Phaseolus ssp.) as a model for understanding crop evolution.
doi: 10.1093/molbev/msq148 Front. Plant Sci. 8:722. doi: 10.3389/fpls.2017.00722
Alexander, D. H., Novembre, J., and Lange, K. (2009). Fast model-based estimation Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.-H., Xie, D., et al. (2014).
of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. doi: 10.1101/ BEAST 2: a software platform for bayesian evolutionary analysis. PLOS Comp.
gr.094052.109 Biol. 10:e1003537. doi: 10.1371/journal.pcbi.1003537
Angelis, K., and Dos Reis, M. (2015). The impact of ancestral population size and Bouckaert, R. R. (2010). DensiTree: making sense of sets of phylogenetic
incomplete lineage sorting on Bayesian estimation of species divergence times. trees. Bioinformatics 26, 1372–1373. doi: 10.1093/bioinformatics/
Curr. Zool. 61, 874–885. doi: 10.1093/czoolo/61.5.874 btq110
Angioi, S. A., Desiderio, F., Rau, D., Bitocchi, E., Attene, G., and Papa, R. Bryant, D., Bouckaert, R., Felsenstein, J., Rosenberg, N. A., and RoyChoudhury, A.
(2009). Development and use of chloroplast microsatellites in Phaseolus spp. (2012). Inferring species trees directly from biallelic genetic markers: bypassing
and other legumes. Plant Biol. 11, 598–612. doi: 10.1111/j.1438-8677.2008. gene trees in a full coalescent analysis. Mol. Biol. Evol. 29, 1917–1932.
00143.x doi: 10.1093/molbev/mss086
Bataillon, C. (1972). La Ciudad y el Campo en el México Central. Available Búrquez, A., and Sarukhán, J. (1984). Biología floral de poblaciones silvestres de
at: https://books.google.co.in/books/about/La_ciudad_y_el_campo_en_el_M% Phaseolus coccineus L. II. Sistemas reproductivos. Bol. Soc. Bot. México 46,
C3%A9xico_Centr.html?hl=&id=h_grAAAAMAAJ&redir_esc=y 3–12.

Frontiers in Plant Science | www.frontiersin.org 13 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., and Lee, J. J. Instituto Nacional de Ecología (INE) (2008). Ecorregiones Terrestres de México.
(2015). Second-generation PLINK: rising to the challenge of larger and richer Available at: http://www.conabio.gob.mx/informacion/gis/
datasets. Gigascience 4:7. doi: 10.1186/s13742-015-0047-8 Koinange, E. M. K., Singh, S. P., and Gepts, P. (1996). Genetic control of
CONABIO and IUCN (2016). Resultados del Segundo Taller del Proyecto the domestication syndrome in common bean. Crop Sci. 36, 1037–1045.
Salvaguardando los Parientes Silvestres de Plantas Cultivadas. Gland: Unión doi: 10.2135/cropsci1996.0011183x003600040037x
Internacional para la Conservación de la Naturaleza. Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with Burrows-
Cornuet, J.-M., Pudlo, P., Veyssier, J., Dehne-Garcia, A., Gautier, M., Leblois, R., Wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/bioinformatics/
et al. (2014). DIYABC v2.0: a software to make approximate Bayesian btp324
computation inferences about population history using single nucleotide Li, Y.-H., Zhao, S.-C., Ma, J.-X., Li, D., Yan, L., Li, J., et al. (2013). Molecular
polymorphism, DNA sequence and microsatellite data. Bioinformatics 30, footprints of domestication and improvement in soybean revealed by whole
1187–1189. doi: 10.1093/bioinformatics/btt763 genome re-sequencing. BMC Genomics 14:579. doi: 10.1186/1471-2164-
Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., et al. 14-579
(2011). The variant call format and VCFtools. Bioinformatics 27, 2156–2158. López-Soto, J. L., Ruiz-Corral, J. A., Sánchez-González, J. J., and Lépiz-
doi: 10.1093/bioinformatics/btr330 Ildelfonso, R. (2005). Adaptación climática de 25 especies de frijol
Delgado-Salinas, A., Bibler, R., and Lavin, M. (2006). Phylogeny of silvestre (Phaseolus spp) en la República Mexicana. Rev. Fitotec. Mex. 28,
the genus phaseolus (Leguminosae): a recent diversification in an 221–230.
ancient landscape. Syst. Bot. 31, 779–791. doi: 10.1600/036364406779 Luu, K., Bazin, E., and Blum, M. G. B. (2017). pcadapt: an R package to perform
695960 genome scans for selection based on principal component analysis. Mol. Ecol.
Drummond, A. J., and Rambaut, A. (2007). BEAST: bayesian evolutionary Resour. 17, 67–77. doi: 10.1111/1755-0998.12592
analysis by sampling trees. BMC Evol. Biol. 7:214. doi: 10.1186/1471-2148- Mastretta-Yanes, A., Moreno-Letelier, A., Piñero, D., Jorgensen, T. H., and
7-214 Emerson, B. C. (2015). Biodiversity in the Mexican highlands and
Eguiarte, L. E., Aguirre-Planter, E., Aguirre, X., Colín, R., González, A., Rocha, M., the interaction of geology, geography and climate within the Trans-
et al. (2013). From isozymes to genomics: population genetics and conservation Mexican Volcanic Belt. J. Biogeogr. 42, 1586–1600. doi: 10.1111/jbi.
of Agave in México. Bot. Rev. 79, 483–506. doi: 10.1007/s12229-013-9123-x 12546
Elshire, R. J., Glaubitz, J. C., Sun, Q., Poland, J. A., Kawamoto, K., Buckler, Maxted, N., Kell, S., Ford-Lloyd, B., Dulloo, E., and Toledo, Á. (2012). Toward
E. S., et al. (2011). A robust, simple genotyping-by-sequencing (GBS) approach the systematic conservation of global crop wild relative diversity. Crop Sci. 52,
for high diversity species. PLOS ONE 6:e19379. doi: 10.1371/journal.pone.001 774–785. doi: 10.2135/cropsci2011.08.0415
9379 Meyer, R. S., and Purugganan, M. D. (2013). Evolution of crop species: genetics of
Escalante, A. M., Coello, G., Eguiarte, L. E., and Pinero, D. (1994). Genetic domestication and diversification. Nat. Rev. Genet. 14, 840–852. doi: 10.1038/
structure and mating systems in wild and cultivated populations of Phaseolus nrg3605
coccineus and P. vulgaris (Fabaceae). Am. J. Bot. 81:1096. doi: 10.2307/24 Monroy, R., and Quezada-Martínez, A. (2010). Estudio Etnobotánico del frijol
45471 Yepatlaxtle (Phaseolus coccineus L.), en el Área Natural Protegida Corredor
Foll, M., and Gaggiotti, O. (2008). A genome-scan method to identify Biológico Chichinautzin, Morelos, México. Cuernavaca: Universidad Autónoma
selected loci appropriate for both dominant and codominant markers: del Estado de Morelos.
a Bayesian perspective. Genetics 180, 977–993. doi: 10.1534/genetics.108. Morrell, P. L., Buckler, E. S., and Ross-Ibarra, J. (2011). Crop genomics: advances
092221 and applications. Nat. Rev. Genet. 13, 85–96. doi: 10.1038/nrg3097
Freytag, G. F., and Debouck, D. G. (2002). Taxonomy, Distribution, and Ecology Papa, R., and Gepts, P. (2003). Asymmetry of gene flow and differential
of the Genus Phaseolus (Leguminosae-Papilionoideae) in North America, Mexico geographical structure of molecular diversity in wild and domesticated
and Central America. Forth Worth, TX: Botanical Research Institute of Texas common bean (Phaseolus vulgaris L.) from Mesoamerica. Theor. Appl. Genet.
(BRIT). 106, 239–250. doi: 10.1007/s00122-002-1085-z
Gepts, P. (2014). The contribution of genetic and genomic approaches to plant Phytozome (2017). Available at: http://phytozome.jgi.doe.gov/ [accessed June 26,
domestication studies. Curr. Opin. Plant Biol. 18, 51–59. doi: 10.1016/j.pbi.2014. 2017].
02.001 Piñero, D., Caballero-Mellado, J., Cabrera-Toledo, D., et al. (2009). “La diversidad
Glaubitz, J. C., Casstevens, T. M., Lu, F., Harriman, J., Elshire, R. J., Sun, Q., genética como instrumento para la conservación y el aprovechamiento de la
et al. (2014). TASSEL-GBS: a high capacity genotyping by sequencing analysis biodiversidad: estudios en especies mexicanas,” in Capital Natural de México,
pipeline. PLOS ONE 9:e90346. doi: 10.1371/journal.pone.0090346 p. 619. Comisión Nacional para el Conocimiento y Uso de la Biodiversidad, eds J.
González, A. M., De Ron, A. M., Lores, M., and Santalla, M. (2014). Effect Sarukhán, J. Soberón, G. Halffter and J. Llorente Bousquets (Mexico: Comisión
of the inbreeding depression in progeny fitness of runner bean (Phaseolus Nacional para el Conocimiento y Uso de la Biodiversidad).
coccineus L.) and it is implications for breeding. Euphytica 200, 413–428. Price, M. N., Dehal, P. S., and Arkin, A. P. (2009). FastTree: computing large
doi: 10.1007/s10681-014-1177-2 minimum evolution trees with profiles instead of a distance matrix. Mol. Biol.
Goudet, J. (2005). hierfstat, a package for R to compute and test hierarchical Evol. 26, 1641–1650. doi: 10.1093/molbev/msp077
F-statistics. Mol. Ecol. Notes 5, 184–186. doi: 10.1111/j.1471-8286.2004. R Core Team (2017). R: A Language and Environment for Statistical Computing.
00828.x Vienna: R Foundation for Statistical Computing.
Hammer, K. (1984). Das Domestikationssyndrom. Kulturpflanze 32, 11–34. doi: Rannala, B., and Yang, Z. (2003). Bayes estimation of species divergence times and
10.1007/BF02098682 ancestral population sizes using DNA sequences from multiple loci. Genetics
Hufford, M. B., Lubinksy, P., Pyhäjärvi, T., Devengenzo, M. T., Ellstrand, 164, 1645–1656.
N. C., and Ross-Ibarra, J. (2013). The genomic signature of crop-wild Renaut, S., and Rieseberg, L. H. (2015). The accumulation of deleterious mutations
introgression in maize. PLOS Genet. 9:e1003477. doi: 10.1371/journal.pgen.10 as a consequence of domestication and improvement in sunflowers and
03477 other compositae crops. Mol. Biol. Evol. 32, 2273–2283. doi: 10.1093/molbev/
Hufford, M. B., Xu, X., van Heerwaarden, J., Pyhäjärvi, T., Chia, J.-M., msv106
Cartwright, R. A., et al. (2012). Comparative population genomics of maize Rodiño, A. P., Paula Rodiño, A., Lema, M., Pérez-Barbeito, M., Santalla, M., and
domestication and improvement. Nat. Genet. 44, 808–811. doi: 10.1038/ De Ron, A. M. (2006). Assessment of runner bean (Phaseolus coccineus L.)
ng.2309 germplasm for tolerance to low temperature during early seedling growth.
Huson, D. H., and Bryant, D. (2006). Application of phylogenetic networks Euphytica 155, 63–70. doi: 10.1007/s10681-006-9301-6
in evolutionary studies. Mol. Biol. Evol. 23, 254–267. doi: 10.1093/molbev/ Rodriguez, M., Rau, D., Angioi, S. A., Bellucci, E., Bitocchi, E., Nanni, L.,
msj030 et al. (2013). European Phaseolus coccineus L. landraces: population structure
Instituto Nacional de Estadística, Geografía e Informática (INEGI), Comisión and adaptation, as revealed by cpSSRs and phenotypic analyses. PLOS ONE
Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO), and 8:e57337. doi: 10.1371/journal.pone.0057337

Frontiers in Plant Science | www.frontiersin.org 14 November 2017 | Volume 8 | Article 1891

Guerra-García et al. Domestication Genomics of Phaseolus coccineus

Salinas, A. D. (1988). “Variation, taxonomy, domestication, and germplasm Vargas-Vázquez, L. P., Muruaga-Martínez, J. S., Lépiz-Ildefonso, R., and Pérez-
potentialities in phaseoluscoccineus,” in Genetic Resources of Phaseolus Beans. Guerrero, A. (2012). La colección INIFAP de frijol ayocote (Phaseolus coccineus
Current Plant Science and Biotechnology in Agriculture, Vol. 6, ed. P. Gepts L.) I. Distribución geográfica de sitios de colecta. Rev. Mex. Cien. Agríc. 3,
(Dordrecht: Springer), 441–463. doi: 10.1007/978-94-009-2786-5_18 1247–1259.
Schmutz, J., Cannon, S. B., Schlueter, J., Ma, J., Mitros, T., Nelson, W., et al. Zheng, X., Levine, D., Shen, J., Gogarten, S. M., Laurie, C., and Weir, B. S.
(2010). Genome sequence of the palaeopolyploid soybean. Nature 463, 178–183. (2012). A high-performance computing toolset for relatedness and principal
doi: 10.1038/nature08670 component analysis of SNP data. Bioinformatics 28, 3326–3328. doi: 10.1093/
Schmutz, J., McClean, P. E., Mamidi, S., Wu, G. A., Cannon, S. B., Grimwood, J., bioinformatics/bts606
et al. (2014). A reference genome for common bean and genome-wide analysis Zhou, Z., Jiang, Y., Wang, Z., Gou, Z., Lyu, J., Li, W., et al. (2015). Resequencing
of dual domestications. Nat. Genet. 46, 707–713. doi: 10.1038/ng.3008 302 wild and cultivated accessions identifies genes related to domestication
Schwember, A. R., Carrasco, B., and Gepts, P. (2017). Unraveling agronomic and and improvement in soybean. Nat. Biotechnol. 33, 408–414. doi: 10.1038/nbt.
genetic aspects of runner bean (Phaseolus coccineus L.). Field Crops Res. 206, 3096
86–94. doi: 10.1016/j.fcr.2017.02.020
Services | Institute of Biotechnology (2017). Available at: http://www.biotech. Conflict of Interest Statement: The authors declare that the research was
cornell.edu/brc/genomic-diversity-facility/services [accessed June 26, 2017]. conducted in the absence of any commercial or financial relationships that could
Sokal, R. R., and Rohlf, F. J. (1995). Biometry: The Principles and Practice of Statistics be construed as a potential conflict of interest.
in Biological Research, 3rd Edn. New York: W. H. Freeman and Co.
Spataro, G., Tiranti, B., Arcaleni, P., Bellucci, E., Attene, G., Papa, R., et al. The handling Editor declared a shared affiliation and past co-authorship, though
(2011). Genetic diversity and structure of a worldwide collection of Phaseolus no other collaboration, with the authors and states that the process nevertheless
coccineus L. Theor. Appl. Genet. 122, 1281–1291. doi: 10.1007/s00122-011- met the standards of a fair and objective review.
1530-y
Tang, H., Sezen, U., and Paterson, A. H. (2010). Domestication and plant genomes. Copyright © 2017 Guerra-García, Suárez-Atilano, Mastretta-Yanes, Delgado-Salinas
Curr. Opin. Plant Biol. 13, 160–166. doi: 10.1016/j.pbi.2009.10.008 and Piñero. This is an open-access article distributed under the terms of the Creative
van Heerwaarden, J., Doebley, J., Briggs, W. H., Glaubitz, J. C., Goodman, Commons Attribution License (CC BY). The use, distribution or reproduction in
M. M., de Jesus Sanchez Gonzalez, J., et al. (2011). Genetic signals of other forums is permitted, provided the original author(s) or licensor are credited
origin, spread, and introgression in a large sample of maize landraces. and that the original publication in this journal is cited, in accordance with accepted
Proc. Natl. Acad. Sci. U.S.A. 108, 1088–1092. doi: 10.1073/pnas.10130 academic practice. No use, distribution or reproduction is permitted which does not
11108 comply with these terms.

Frontiers in Plant Science | www.frontiersin.org 15 November 2017 | Volume 8 | Article 1891