Consequences of introgression and gene flow on the genetic structure and diversity of Lima bean (Phaseolus lunatus L.) in its Mesoamerican diversity area

Study sites and plant material

The study was conducted in the Yucatan Peninsula, Mexico, at a local and regional scale:

(1) Local scale. A specific search for wild-weedy-crop complexes of Lima bean was carried out in 2019; two complexes were found, which were named as Itzinté and Dzitnup based on local place names (from this moment on, the use of these names will only refer to the existence of the complexes, unless otherwise indicated). As species, P. lunatus is not in any risk category, all material plant used in this study was collected under scientific permits (SGPA/DGVS/008421/18 and SGPA/DGGFS/712/2913/17) issued by the Secretaría del Medio Ambiente y Recursos Naturales (SEMARNAT-Mexico) to researchers at Centro de Investigación Científica de Yucatán, A. C (CICY).

Itzinté. Located in the abandoned Mayan ruins (N 20.018611 and W −89.724444; Fig. 1, Table S1), near de town of Bolonchén, located at the municipality of Hopelchén in the Campeche state. In this complex, wild and domesticated populations of Lima bean have coexisted for at least 20 years (Martínez-Castillo et al., 2004), but probably for longer, since this is an important cultivation area for farmers of Bolonchén, founded in 1957. Itzinté comprises an area of flat land interrupted only by mounds of the Mayan ruins and some hills, it has deep soils where the so-called continuous milpa is practiced (i.e., plantation fields are used for many years continuously without a rest period). Within Itzinté, seeds were collected in three sites separated by a distance of 300–400 m: (a) site 1, the slope of a hill where only some wild plants were found; (b) site 2, a milpa where domesticated and wild plants were observed growing in close proximity (<1 m); and (c) site 3, a milpa with only one production cycle and whose seeds were brought from the town of Escárcega, located in the south of Campeche. A total of 27 seeds were sequenced from this complex, all from different plants: (a) three wild seeds from site 1, without morphological evidence of introgression; (b) 22 seeds from site 2, which spanned the entire range of morphological variation found and with possible introgression since they were collected from wild and domesticated plants growing near each other; and (c) two domesticated seeds from site 3, without morphological evidence of introgression (Fig. 2).

Dzitnup. Located near the town of Dzitnup, located in the municipality of Hecelchakán in Campeche state (Fig. 1, Table S1). Within this complex, seeds were collected in three sites separated by a distance of 1.5–2 km: (a) site 1, Xcalumkin archaeological zone (N 20.171806 and W −90.009722), where a wild population was found; (b) site 2, milpa located on the road from Dzitnup to Xcalumkin (N 20.163611 and W −90.05333), where domesticated and wild plants were growing intermingled or in close proximity (<1 m); the characteristics of the terrain (shallow soil with high stoniness) indicate that this is a milpa with 1–2 years of use; and (c) site 3, home garden at the exit of the town of Dzitnup (N 20.20011 and W −90.11611), heading for Xcalumkin. A total of 45 seeds from different plants were sequenced from this complex: (a) nine wild seeds collected at site 1, without morphological evidence of introgression; (b) 33 seeds collected at site 2, which spanned the entire range of morphological variation found and with possible introgression since they were collected from wild and domesticated plants growing near each other; and (c) three domesticated seeds collected at site 3, without morphological evidence of introgression (Fig. 2). In short, the SNP data set used for the analyses at the local scale considered 72 accessions (27 from Itzinté and 45 from Dzitnup).

(2) Regional scale

A total of 65 accessions in the broader Yucatan Peninsula were sequenced: 39 domesticated and 25 wild accessions; part of these accessions were collected in 2019 and others were obtained from the collection of Dr. Jaime Martínez-Castillo. In addition, we included data reported by Chacón-Sánchez & Martínez-Castillo (2017) comprising 13 domesticated and six wild accessions from this region, for a total of 83 accessions: 52 domesticated and 31 wild (Fig. 1, Table S2). Moreover, we incorporated a control group composed of samples from Chacón-Sánchez & Martínez-Castillo (2017) that did not include the Yucatan Peninsula: nine wild MI accessions (from Mexico), eight domesticated MI accessions (five from Mexico, two from Guatemala, and one from El Salvador), two wild MII accessions (one from Mexico and one from Guatemala), one domesticated MII accession from Costa Rica, and eight accessions from the Andean gene pool (five domesticated and three wild accessions from Ecuador and Peru) (Table S2). The purpose of incorporating the control group was to help in the assignment of individuals to one or several gene pools and in the prediction of chromosomal segments with possible introgression events. In short, the SNP data set used for the analyses at the regional scale considered 66 domesticated and 45 wild accessions.

DNA extraction

Seedlings from the selected accessions were used to obtain DNA. Total genomic DNA was extracted from leaf tissue using a silica extraction protocol (Echevarría-Machado et al., 2005). DNA quality was verified by electrophoresis and visualization on 1% agarose gels, and DNA quantification was carried out in a Quantus ™ fluorometer. Genotyping-by-sequencing (GBS) library construction and 150 bp paired-end sequencing in an Illumina Hi-seq 2000 was performed at the Australian Genome Research Facility (Melbourne, Australia). Total genomic DNA extracted from each sample was digested with the Ape KI restriction enzyme which resulted in a population of fragments, most of them in the range of 200–300 bp. For library preparation, 100 ng of digested DNA and 3.6 ng of each adaptor were used for ligation. After ligation, 5 µl from each one of 96 libraries were mixed, amplified by PCR and sequenced in an Illumina sequencing lane. Average sequencing depth per accession ranged from 5X to 10X.

SNP detection

The paired-end sequences obtained from the Illumina platform were de-multiplexed by their bar code using the Demultiplex option of the program NGSEP ver. 4.02 (graphical interface) (Duitama et al., 2014). Subsequently, data were cleaned using Trimmomatic v.0.36 (Bolger, Lohse & Usadel, 2014) to remove adapters and low-quality sequences (>Q30). Afterwards, reads were aligned to the Lima bean reference genome (available at: https://phytozome-next.jgi.doe.gov/info/Plunatus_V1) with the Bowtie2 tool and with default parameters in NGSEP (ver. 4.02); we then used the MultiSampleVariantsDetector option of NGSEP (ver. 4.02) to identify and genotype the genetic variants with the following parameters: 100 as the maximum number of alignments per start position, SNP heterozygosity rate of 0.0001 and minimum genotyping quality of Q40. The SNP matrix was filtered with the FilterVCF option of NGSEP (ver. 4.02) with the following options: MAF of 5%, missing data at a maximum of 20% (for loci and genotypes), minimum genotyping quality of Q40, maximum heterozygosity of 0.2 and removing SNPs from the repetitive regions of the genome. This filtering process resulted in 15,168 SNPs and 184 accessions. SNP density per chromosome varied from 15.2 SNPs/Mbp for chromosome PL10 to 41.8 SNPs/Mbp for chromosome PL03, with an average SNP density of 29 SNPs/Mbp. The VCF file was converted to other formats using TASSEL (Bradbury et al., 2007) and PGDSpider (Lischer & Excoffier, 2012).

Raw sequence data generated and used in this study have been deposited in the NBCI Sequence Read Archive (SRA) under the Bioproject accessions numbers PRJNA596114 and PRJNA823361. Details about the SRX number for each accession can be found in the Table S3. Also, data from both complexes analyzed and those from the Yucatan peninsula were included as Supplemental Information as a *.vcf file.

Data analysis

Introgression at local scale

In Itzinté, optimal K was 2; STRUCTURE recognized genetically pure wild (green color) and domesticated (red color) individuals and identified seven individuals collected at site 2 that showed different grades of ancestry of both gene pools, which can be considered as admixed individuals (Fig. 3A). In Dzinup, K optimal was also 2; STRUCTURE recognized pure wild (green color) and domesticated (red color) individuals and identified 19 individuals collected at site 2 that showed ancestry of both gene pools, which can be considered as admixed individuals (Fig. 3B). Figure 3C showed that, in both complexes, for domesticated individuals the percentage of ancestry from the wild gene pool is, in most cases, low; on the contrary, in wild individuals, the percentage of ancestry from the domesticated gene pool is more variable and includes some individuals that show high values of ancestry from the domesticated gene pool. For Itzinté, the PCoA showed two compact groups of accessions (wild and domesticated), and some accessions with a scattered distribution that, according with STRUCTURE, correspond to admixed accessions; the first principal coordinate explained 54.13% of variation, followed by 7.27% for the second one, both coordinates explained 61.4% of the total variation (Fig. S1). Similarly, for Dzitnup, the PCoA showed two compact groups (wild and domesticated), with the rest of the accessions showing a scattered distribution (but minor compared with Itzinté) that, according with STRUCTURE, are also admixed accessions; the first principal coordinate explained 28.51% of variation, followed by 10.46% for the second one, both coordinates explained only 38.97% of the total variation (Fig. S2).

For Itzinté, NGSEP identified 54 introgression events in four domesticated and seven wild individuals, all of them (except two) collected at site 2 (Figs. 4C and 4D; Table S4); these individuals corresponded with the seven admixed individuals observed with STRUCTURE. 33 of these events corresponded to introgressions from the domesticated to the wild gene pool. For Dzitnup, NGSEP identified 118 introgression events in 5 domesticated and 27 wild individuals, 26 of them collected at site 2 (Figs. 4A and 4B; Table S4); these individuals corresponded to 15 of the 19 admixed individuals observed with STRUCTURE. Of the 118 introgression events, 94 were in the domesticated to wild direction, and 24 events were in the opposite direction. It is worth mentioning that 11 wild individuals from Dzitnup share an introgression block from the domesticated gene pool at chromosome 3 of about 2 Mbp (from base 34,238,474 to base 36,221,738).

Although both STRUCTURE and NGSEP analyses indicated gene flow mainly from crop-to-wild in both complexes, these analyses also showed some evidence of gene flow in both directions. In order to delve into this, we conducted ABBA-BABA tests to assess whether the signals of introgression are stronger between wild and domesticated populations within the same complex than with other wild and domesticated populations of the Yucatan Peninsula. In effect, as it can be seen in Fig. S3 and Table S5, D values and f4-ratios are higher between domesticated and wild populations within the same complex (models m1, m3, m5 and m7) than with other wild/domesticated populations in the Peninsula of Yucatan (models m2, m4 and m6, except for model 8). The proportion of introgression, as measured by the f4-ratio, was higher in the older complex (for Itzinté was around 24%) than in the younger complex (for Dzitnup was around 15%).

To further investigate the introgression of domesticated alleles into wild populations (d-w), and vice versa (w-d), within both complexes at a fine genomic scale, we chose to investigate four models: m1, m3, m5 and m7. Fig. S4 and Table S6 show the distribution of D and df values in all eleven chromosomes in Itzinté and Dzitnup and the proportion of introgression in the direction wild to domesticated and vice versa. It can be observed that median df values for all four models were positive, indicating a higher contribution of domesticated/wild alleles in P2 (wild/domesticated populations within the same complex) than in P1 (wild MII or MI populations outside the complexes). For the older complex (Itzinté), median df values are very similar for models m1 and m3, indicating that the proportion of introgression is similar (df is around 12–13%) in both directions (d-w, w-d). In contrast, for the younger complex (Dzitnup) the proportion of introgression in the direction domesticated to wild (median df = 18%) is almost twice than in the opposite direction (median df = 10%). When comparing introgression proportions (median df) in both directions (d-w, w-d) at the chromosome level, df values are strikingly different in chromosome seven for Itzinté and in chromosome nine for Dzitnup, with df values larger in the direction domesticated to wild in both chromosomes (for details between the chromosomal regions found with evidence of introgression and its relation to domestication, see Table S6.1).

Figure 5 shows the distribution of introgressed blocks in both complexes. It should be noted that introgression blocks in Fig. 5 are smaller than in Fig. 4 because Dsuite calculates introgression at the population level (not at the individual level as NGSEP does) and at a finer genomic scale in windows of 10 SNPs/1 SNP. In Itzinté (the older complex) the number of introgressed blocks was larger (189 blocks) than in Dzitnup (the younger complex) (129 blocks), and the size of the genome occupied by all the introgressed blocks was larger in Itzinté than in Dzitnup (195 Mb versus 158 Mb) (Tables S7–S8). In contrast, the average size of the individual blocks was larger in Dzitnup (1.2 Mb/block) than in Itzinté (1.0 Mb/block). All this indicates that in the older complex, the introgressed blocks are more in number but smaller in length compared to the younger complex. Also, for Dzitnup the size of the genome occupied by introgressed blocks was larger in the direction d-w (88 Mb) than in the w-d direction (70 Mb), suggesting gene flow may be predominant from landraces to wild populations in this complex. The opposite pattern occurred in Itzinté (76 Mb and 119 Mb, respectively), suggesting gene flow may be predominant from wild to landrace populations.

In relation to the genetic diversity found in both complexes, H_E values were slightly higher in the domesticated accessions (Dzitnup = 0.36; Itzinté = 0.34) than in the wild ones (Dzitnup = 0.28, Itzinté = 0.31); but the other estimators were greater in the wild accessions than in the domesticated ones (Table 1), although the statistical significance for these analyses was not evaluated.

Structure, genetic diversity and wild-crop gene flow in Lima bean at regional scale

When the Lima bean accessions of the Yucatan Peninsula were analyzed together with the control group accessions, the Neighbor-Joining (N-J) analysis (Fig. 6) clearly separated all the accessions in two groups, according to the main gene pools of Lima bean: MI and MII. Into the MI group, two subgroups were observed: the green one integrated for all domesticated accessions from the Yucatan Peninsula and those collected in other regions of Mexico and central America, and the blue subgroup integrated for all wild accessions collected out of the Yucatan Peninsula. Into the MII group (indicated in red), we also observed two subgroups: one composed of 12 accessions from the Yucatan Peninsula plus one accession from central Mexico and another from Costa Rica; and another composed of 19 wild accessions collected in Tabasco and Campeche, one collected in Guatemala, plus a domesticated accession collected in Campeche. The Principal Coordinate Analysis (PCoA) (Fig. S5) generated a grouping pattern similar to that obtained with the N-J analysis. The first principal coordinate explained 45.69% of variation, followed by 17.02% for the second one; together, both coordinates explained 62.71% of the total variation. The optimal K obtained was 2; STRUCTURE analysis with K = 2 (Fig. S6A) grouped the wild MII accessions collected in the Yucatan Peninsula with the Andean accessions (green group), while all wild and domesticated accessions of the MI group were grouped together (red group). STRUCTURE with K = 3 (Fig. S6B) grouped the accessions according to the three large gene pools of Lima bean: MII (red group), MI (blue group) and Andean (green group). Finally, when STRUCTURE was run with K = 4 (Fig. S6C), the clustering pattern generated resulted very similar to the obtained with the N-J and the PCoA analyses.

Considering only the Lima bean accessions of the Yucatan Peninsula, the optimal K obtained was 3; STRUCTURE separated the domesticated accessions (all integrated in one group) from the wild ones (integrated in two groups) (Fig. 7): a pure wild group was integrated by 15 wild accessions collected in Tabasco and western Campeche, and an admixed wild group was integrated by 17 wild accessions collected in Yucatan, Quintana Roo and central-eastern Campeche. PCoA showed the existence of three groups (Fig. S7), supporting the STRUCTURE results; the first principal coordinate explained 61.51% of variation, followed by 9.97% for the second one, both coordinates explained 71.48% of the total variation. Likewise, the AMOVA supported the existence of three groups, showing that the strongest genetic differentiation is found between groups (81.74%) (Table S8).

The pure wild group from the Yucatan Peninsula showed a greater value of H_E (0.36) than the admixed wild group (H_E = 0.26), and the latter showed a slightly higher value than the domesticated group (H_E = 0.23); the admixed wild group showed higher values for the other estimators (P% and π), and the domesticated group showed the lowest values for all estimators (Table 1). Statistical differences were not analyzed.

In relation to the recent wild-crop gene flow in the Lima bean of the Yucatan Peninsula, STRUCTURE (K = 3) separated the wild accessions from the domesticated ones, with most accessions showing high ancestry coefficients (>90%). With K = 3, a group of wild accessions considered as genetically pure were collected in Tabasco and western Campeche where Lima bean is not planted and clustered apart from another group of wild accessions that showed some ancestry from the domesticated gene pool, and that were collected in areas of traditional agriculture of Yucatan, Quintana Roo and central-eastern Campeche, where wild and domesticated Lima bean populations grow in sympatry. Moreover, only three domesticated accessions showed some degree of ancestry from the wild gene pool, although it was low (Fig. 7).

Taking in account the existence of three genetic groups into the Yucatan peninsula, MIGRATE showed a higher historical mutation scaled immigration rate (M = m/µ) from domesticated group (157.17) to admixed wild group (77. 07), slightly higher from the pure wild group (130.61) to admixed wild group (123.14), and higher from the domesticated group (114.07) to the pure wild group (55.93).

Introgression and gene flow in Lima bean of the Yucatan Peninsula

STRUCTURE and NGSEP analyses showed predominantly a crop-to-wild introgression at the local scale, but also some evidence of a bidirectional gene flow in both studied complexes; ABBA-BABA tests showed bidirectional gene flow in both complexes, but predominantly in the direction wild-to-crop in Itzinté and in the opposite direction in Dzitnup. These tests also showed a higher introgression proportion in Itzinté than in Dzitnup.

Itzinté is the oldest of the two complexes, it has existed for at least 20 years (Martínez-Castillo et al., 2004), suggesting that the wild-crop gene flow could have been occurring continuously during this period. There, the practice of continuous milpa has generated an environment in permanent disturbance that has favored the presence of wild populations due to two characteristics of these ones: (1) a metapopulation dynamic (Barrantes et al., 2008), which has been identified as one of the factors that affect the dynamics of wild-crop introgression (Meirmans, Bousquet & Isabel, 2008; Ellstrand et al., 2013); and (2) a seed bank with an up to three years viability (Degreef et al., 2002). In this regard, some models predict that the dynamics of the seed bank and the establishment of seedlings under competition represent the strongest determinants of the persistence of weedy populations that occupy disturbed areas (Linder & Schmitt, 1994; Linder & Schmitt, 1995). In Itzinté, the fate of the resulting introgressed seeds can follow two possible scenarios: (1) if human selection is exerted, both directly on the plant material or indirectly by significant modification of the environment, the complex will progressively include only domesticated-type individuals; if farmers benefit from the generated variants, the presence of this type of complexes could be indicative of an on-going domestication process (Beebe et al., 1997); and (2) plants within the complex could evolve into wild forms under natural selection; it is known that selection within hybrid zones can act in two ways: (a) endogenous selection (not dependent on the environment), which assumes that hybrid generations may have inherent and consistent fitness advantages or disadvantages in relation to their parents and each other, and (b) exogenous selection (dependent on the environment), which considers that selection operating on hybrids may differ according to environmental conditions, implying the existence of genotype-environment interactions (Campbell & Waser, 2007; Kimball, Campbell & Lessin, 2008).

In Dzitnup, the wild-crop contact site was located in a slash-and-burn milpa. This type of milpa is commonly used only for 1–3 years and then left for a rest period (fallow) that can last up to 15 years (Hernández-Xolocotzi, 1959). It means that this complex will not remain for long, perhaps 2–3 years while the milpa is in use and one more year after production ends, allowing not harvested seeds to fall, germinate, reach the reproductive stage and breed with wild plants; it is the most common scenario for wild-crop introgression events in Lima bean of the Yucatan Peninsula. Characteristics of this complex suggests two possible scenarios for the resulting introgressed seeds: (1) if the seeds are harvested, they will be subject to human selection; current knowledge indicates that the genetic basis of domestication traits involve generally few genomic regions and therefore can be purged quite rapidly (Papa et al., 2005; Papa et al., 2007); and (2) if the introgressed seeds remain in the field, natural selection will act on the resulting hybrids and their maintenance will depend on their fitness (Mercer et al., 2014).

In both complexes, the identified introgressed individuals were mostly collected at the sites where wild and domesticated plants grew very close to each other (sites # 2). In hybrid zones, the relative frequencies of hybrid versus non-hybrid genotypes are expected to vary according to some characteristics, including distance between hybridizing individuals, relative size of each population, overlap in the flowering period, and pollen and seed dispersal (Nagy, 1997; Ellstrand, 2003; Jenczewski, Ronfort & Chévre, 2003). These characteristics were not evaluated in our study, but we observed that the wild and domesticated individuals in Dzitnup grew intermingled, while in Itzinté both types of individuals grew, mainly, in adjacent places. In relation to the flowering period of wild and domesticated populations of Lima bean, they overlap throughout the Yucatan Peninsula (Martínez-Castillo et al., 2004) and it is estimated that pollen and seed dispersal does not exceed 6 m, with some exceptions (Hardy et al., 1997; Baudoin et al., 1998); furthermore, both complexes were located in the region known as Los Chenes, Campeche, separated by approximately 40 km, and under very similar ecological-environmental conditions.

The asymmetry in the wild-crop introgression found at local scale (mainly crop-to-wild in Dzitnup, where more wild accessions were analyzed than domesticated ones, compared with Itzinté), is expected when the relative size of wild or domesticated populations is different, or when populations are fragmented (Ellstrand, 1992; Ellstrand, Prentice & Hancock, 1999); unfortunately, we did not evaluate these factors. Dzul-Tejero, Coello-Coello & Martínez-Castillo (2014), evaluated introgression in three Lima bean complexes in the Yucatan Peninsula using 11 SSR loci and ancestry coefficients obtained with STRUCTURE; these authors reported significant differences in the degree of introgression among the three complexes studied, the complex with the least introgression was a milpa where the farmer recognized and eliminated introgressed seeds, while the complex with the highest introgression was a milpa where the farmer recognized and positively selected the introgressed seeds; but did not find an asymmetry within any of the complexes. However, in this study we used a wider genomic sampling (15,168 SNPs) and the detection of introgression events at a finer level (chromosome); furthermore, the two studies were conducted in different areas of the Yucatan Peninsula and did not involve the same wild populations or the same landraces, so different factors (e.g., agricultural, genetic or environmental management) may explain the differences observed in both studies.

At the regional scale, gene flow showed a crop-to-wild asymmetry, according to STRUCTURE; however, this asymmetry was detected only in wild accessions that were collected in areas of traditional agriculture where wild and domesticated populations coexist, since the wild accessions collected in Tabasco and western Campeche (where Lima bean is not planted) did not showed genetic infiltration of the domesticated gene pool. This crop-to-wild asymmetry was also supported with the results from MIGRATE, which showed more historical gene flow from the domesticated gene pool to the wild one; the admixed wild group appears to be acting as a genetic bridge between the domesticated populations and the pure wild populations. Martínez-Castillo et al. (2007), using nine loci of SSRs, also reported an asymmetric historical gene flow in Lima bean from the Yucatan Peninsula, three times greater from the domesticated towards the wild gene pool. The crop-to-wild asymmetry observed in Lima bean has also been reported in common bean (Phaseolus vulgaris L.) from other regions of Mexico and the Americas (Beebe et al., 1997; Papa & Gepts, 2003; Papa et al., 2005; Zizumbo-Villarreal et al., 2005; Hoc et al., 2006).

Crop-to-wild gene flow and introgression asymmetry, as observed at the regional scale (and into Dzitnup), may depend on the fitness of several hybrid generations in relation to their wild counterpart (Ellstrand et al., 2013), which can also be strongly influenced by the environment (Mercer, Wyse & Shaw, 2006). In species such as Lima bean, different aspects related to agricultural management could also help explain the asymmetry: (1) In cultivated environments, farmers who grow common bean easily recognize and select against F1 hybrids, since its seeds generally have an intermediate size compared to the parentals, and a different color compared to the seeds of the domesticated maternal plant; however, in wild environments, the domesticated traits tend to be hidden by the wild alleles and natural selection will only operate in later generations (Papa & Gepts, 2003; Papa et al., 2005). In Lima bean, the recognition of hybrids is also based on the flavor of the harvested seeds. Wild seeds contain high concentrations of linamarin, a cyanogenic glycoside that makes them inedible (Maquet, 1991) and whose presence allows Mayan farmers to detect unpleasant flavors in the introgressed seeds of their landraces (Martínez-Castillo et al., 2007). (2) The intensification of traditional agriculture in the Yucatan Peninsula has generated an increase in the use of herbicides, a decrease in the rest period, the presence of maize monocultures and the use of agricultural machinery, among other changes (Martínez-Castillo et al., 2004). These changes can generate a reduction in density of wild populations present in domesticated areas, thereby favoring the existence of a higher pollen production of domesticated plants and an increased gene flow to wild plants that remain (Papa & Gepts, 2003; Papa et al., 2005). In the Yucatan Peninsula, some wild individuals have been observed growing near or among hundreds of domesticated Lima bean plants (Martínez-Castillo et al., 2004). (3) In some cities of the Yucatan Peninsula, landraces of Lima bean with white seeds and sweet flavors are preferred, which favors their cultivation; this has led producers in southern Yucatan to eliminate wild populations to avoid introgression with their landraces and thus obtain better sales (Martínez-Castillo et al., 2004); this could help explain the lower degree of genetic infiltration in the domesticated populations of southern Yucatan (Martínez-Castillo et al., 2007).

Several studies have reported the negative aspects of a crop-to-wild asymmetry in relation to the conservation of plant genetic resources, among these, the evolution of new weeds or invasive species (Schierenbeck & Ellstrand, 2009) contributing to the risk of local extinction (Ellstrand & Elam, 1993; Levin, Francis-Ortega & Jansen, 1996), as well as to the risk associated to transgene escape (Gepts & Papa, 2003; Jenczewski, Ronfort & Chévre, 2003; Pilson & Prendeville, 2004; Andow & Zwahlen, 2006). Considering the characteristics of the Mayan milpa, the risk of local extinction of wild Lima bean populations can be considered low, since the dynamics of the milpa allows their existence, which is also favored by their metapopulation behavior. Martínez-Castillo et al. (2006) reported the existence of a positive correlation between agricultural intensification and increase in the genetic diversity of the Lima bean populations, suggesting that these populations are favored by the intensification of disturbance in situations involving at least three years of fallow. Even in complexes like Itzinté (where there is a permanent wild-crop contact), where wild populations could show a greater risk of local extinction due to swamping (i.e., which occurs when a small population loses its genetic integrity and assimilates into a larger population, as a result of repeated episodes of hybridization and introgression; (Ellstrand, Prentice & Hancock, 1999), our results showed high levels of genetic diversity in the wild population, higher than those showed at the regional level. Regarding the negative effects of transgenic crops on wild relatives, it is known that the existence of gene flow and a crop-to-wild asymmetry, as we showed happens in the Lima bean from de Yucatan Peninsula, is necessary for transgene escape (Gepts & Papa, 2003). Currently there are no transgenic varieties of Lima bean, but our study and those of Martínez-Castillo et al. (2007) and Dzul-Tejero, Coello-Coello & Martínez-Castillo (2014) warn of the risks of their possible release in the Yucatan Peninsula.

Contrary, wild-crop introgression may have positive effects on the evolution of domesticated species (Harlan, 1965; Jarvis & Hodgkin, 1999), and wild-weedy-crop complexes are of great importance in the generation of genetic variability in landraces (Pernès & Lourd, 1984; Harlan, 1992; Beebe et al., 1997). Introgression represents one of the main sources of variation for genetic improvement of crops, including that carried out by farmers in traditional agroecosystems (Hajjar & Hodgkin, 2007; Worthington et al., 2012). For Lima bean from the Yucatan Peninsula, high levels of genetic diversity have been reported in individuals grown from milpas with high levels of introgression, compared to milpas with low introgression, as well as the selection and consumption of weedy variants by Mayan farmers (Dzul-Tejero, Coello-Coello & Martínez-Castillo, 2014). This evidence could help explain why the Yucatan Peninsula is the region with the highest number of landraces in Mexico and a center of genetic diversity for this species in Mesoamerica (Ballesteros, 1999; Martínez-Castillo et al., 2004).

Genetic structure and diversity of Lima bean in the Yucatan Peninsula

Using control groups for some analyses helped in the assignment of the individuals collected in the Yucatan Peninsula to the main gene pools of Lima bean (MI or MII); however, these analyses masked the existing genetic structure within this region (in particular, in the wild gene pool where two groups were found). Although we detected a group of admixed wild individuals, Lima bean in the Yucatan Peninsula showed a strong genetic structure based on the existence of wild MII and domesticated MI gene pools. This is an interesting result considering that gene flow and introgression processes between wild and domesticated Lima bean may have been occurring for many centuries in this region as revealed by archaeological evidence (Kaplan & Kaplan, 1988). Even though the self-pollinated tendency of Lima bean limits wild-crop gene flow and could thus explain the observed genetic structure, our study and that of Dzul-Tejero, Coello-Coello & Martínez-Castillo (2014) clearly show the existence of introgression events at a local scale, suggesting that selection is playing an important role in maintaining the genetic identity of wild and domesticated populations that grow in sympatry, as has been reported for common bean (Papa & Gepts, 2003).

In the Yucatan Peninsula we found higher H_E values in the wild gene pool (especially in the Pure wild group) than in the domesticated one. In contrast, although π values were higher for the wild gene pool compared with the domesticated one, for this estimator the Admixed wild group showed a higher nucleotide diversity than the Pure wild group, probably due to the higher percentage of polymorphic sites in the Admixed wild group. Although H_E is a good estimator because it is directly related to genetic diversity and is not influenced by differences in sample size (such as P%), the differences found in H_E values should be taken with caution considering that the statistical significance of this estimator was not evaluated among groups. Differences in genetic diversity between crops and their wild relatives can be explained by a founder effect due to domestication (Ladizinsky, 1985) (and subsequent introduction of the crop into the Peninsula), which has already been reported for this species using SNPs (Chacón-Sánchez & Martínez-Castillo, 2017), ribosomal DNA (Serrano-Serrano et al., 2012) and chloroplast DNA (Andueza-Noh et al., 2013). Chacón-Sánchez & Martínez-Castillo (2017), using 4,593 SNP markers and individuals from the entire Lima bean distribution area in the Americas, reported H_E values that are lower (wild MII = 0.138, domesticated MI = 0.079) than those found by us; although both studies considered different geographical scales, our results suggested that the Yucatan Peninsula is an important area of genetic diversity for Lima bean. Relatively high values of genetic diversity for wild and domesticated populations of Lima bean (considering its endogamous nature) from the Yucatan Peninsula have also been reported using SSR (Martínez-Castillo et al., 2006) and ISSRs (Martínez-Castillo, Colunga-GarcíaMarín & Zizumbo-Villarrreal, 2008; Camacho-Pérez et al., 2017) markers. Even though our study showed a crop-to-wild asymmetry in Lima bean at regional scale, which was also reported by Martínez-Castillo et al. (2007), we did not find evidence of a negative impact on the genetic diversity of the wild gene pool, as expected in crop-to-wild gene flow (Ellstrand, Prentice & Hancock, 1999).