Molecular evidence and ecological niche modeling reveal an extensive hybrid zone among three Bursera species (section Bullockia)

Abstract

The genus Bursera, includes ~100 shrub and trees species in tropical dry forests with its center of diversification and endemism in Mexico. Morphologically intermediate individuals have commonly been observed in Mexican Bursera in areas where closely related species coexist. These individuals are assumed to result from interspecific hybridization, but no molecular evidence has supported their hybrid origins. This study aimed to investigate the existence of interspecific hybridization among three Mexican Bursera species (Bullockia section: B. cuneata, B. palmeri and B. bipinnata) from nine populations based on DNA sequences (three nuclear and four chloroplast regions) and ecological niche modeling for three past and two future scenario projections. Results from the only two polymorphic nuclear regions (PEPC, ETS) supported the hybrid origin of morphologically intermediate individuals and revealed that B. cuneata and B. bipinnata are the parental species that are genetically closer to the putative hybrids. Ecological niche modeling accurately predicted the occurrence of putative hybrid populations and showed a potential hybrid zone extending in a larger area (74,000 km2) than previously thought. Paleo-reconstructions showed a potential hybrid zone existing from the Last Glacial Maximum (~ 21 kya) that has increased since the late Holocene to the present. Future ecological niche projections show an increment of suitability of the potential hybrid zone for 2050 and 2070 relative to the present. Hybrid zone changes responded mostly to an increase in elevational ranges. Our study provides the first insight of an extensive hybrid zone among three Mexican Bursera species based on molecular data and ecological niche modeling.

About the Authors:

Eduardo Quintero Melecio

Contributed equally to this work with: Eduardo Quintero Melecio, Yessica Rico

Roles Data curation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing

Affiliation: Red de Diversidad Biológica del Occidente Mexicano, Centro Regional del Bajío, Instituto de Ecología, A.C., Pátzcuaro, Michoacán, Mexico

Yessica Rico

Contributed equally to this work with: Eduardo Quintero Melecio, Yessica Rico

Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Red de Diversidad Biológica del Occidente Mexicano, Centro Regional del Bajío, Instituto de Ecología, A.C., Pátzcuaro, Michoacán, Mexico

ORCID logo https://orcid.org/0000-0002-0468-8928

Andrés Lira Noriega

Roles Conceptualization, Resources, Supervision, Writing – review & editing

Affiliations Red de Estudios Moleculares Avanzados, Instituto de Ecología, A.C., Xalapa, Veracruz, Mexico, CONACyT, Ciudad de México, Mexico

Antonio González Rodríguez

Roles Resources, Supervision, Writing – review & editing

Affiliation: Laboratorio de Genética de la Conservación, Instituto de Investigaciones en Ecosistemas, Universidad Nacional Autónoma de México, Morelia, Michoacán, Mexico

Introduction

Interspecific gene flow resulting in the formation of intermediate individuals across hybrid zones is a common phenomenon in nature with important consequences for evolution and conservation [1–3]. Hybridization can reinforce the development of reproductive barriers, introduce potentially adaptive genetic variation into a population driving ecological divergence, and lead to the formation of new lineages [4–6]. Alternately, hybridization can result in maladaptive genetic introgression, displacement of pure parental species by novel phenotypes, formation of highly invasive lineages, and species extinction through genetic assimilation [7–9]. The study of hybrid zones has been a prolific area of research since it is key for understanding the factors that promote reproductive isolation, ecological adaptation, and speciation in natural and human-modified habitats [3, 10, 11]. Hybridization is common in plants [12] but there are still many plant taxa in which hybridization is poorly understood [13].

The genus Bursera encompasses a diverse group of nearly 100 species of deciduous and resinous shrub and trees occurring from the south-western USA to Peru, the Galapagos, Bahamas, and the Greater Antilles. The genus comprises two monophyletic sections, Bursera and Bullockia, which diversified in the early Miocene with the increase of aridification in Mesoamerica and expansion of the Tropical Dry Forest (TDF) [14]. The highest species diversity is found in the Pacific drainages of western Mexico, and nearly 90% of the species are endemic to the country [15, 16]. Bursera species have economic and cultural importance as they are used for the extraction of aromatic resins for religious and medicinal purposes [17], for elaboration of wood handcrafts, and for ex-situ propagation in restoration programs or as living fences in vast extensions throughout tropical landscapes [18].

Several taxonomic studies have documented the occurrence of morphologically intermediate individuals in areas where closely related Bursera species co-occur, likely as the result of secondary contact [16, 19–23]. Interspecific hybridization and introgression have been inferred as important drivers of speciation in several Bursera species [16], although only few species have been recognized to have a hybrid origin based on molecular and biochemical evidence [21, 24]. The most relevant study of interspecific hybridization was performed by Weeks and Simpson [21] who confirmed the hybrid origin of three endemic Bursera species (B. brunea, B. gracilipes, and B. ovata) from Hispaniola in the Caribbean based on nuclear and chloroplast DNA sequences. However, a subsequent study by Weeks and Tye [25] did not find evidence for a hybrid zone in the Galapagos Islands between B. malacophylla and B. graveolens, but instead suggested that B. malacophylla and the intermediate individuals shared ancestry of B. graveolens, a species of high morphological variation across the islands. Weeks and Tye [25] implied that the presence of morphologically intermediate individuals in areas of species co-occurrence cannot be taken as unambiguous evidence of interspecific hybridization as many Bursera species exhibit morphological plasticity across their ranges [16, 26].

Besides the use of molecular data in studies of interspecific hybridization, recent approaches have incorporated ecological niche modeling for inferring the geographic occurrence of potential hybrid zones and the role of ecological factors influencing the development of ecological divergence and local adaptation [5, 27–29]. In this context, past and future climatic projections of hybrid zones can provide a complementary perspective for understanding the origin, maintenance, and temporal dynamics of hybrid zones in association with global climatic changes [30, 31].

Bursera species are dominant or codominant woody elements of the TDF, thornscrub, and desertscrub of Mexico [32, 33], and some species can also occur in secondary habitats of oak and pine forests [19]. The TDF of the Bajío in Central-West Mexico was one of the most widespread ecosystems in the region, which is currently disappearing due to anthropogenic land use changes [34]. In this region, between the states of Michoacán and Guanajuato, Rzedowski and Guevara-Féfer [19] reported the ocurrence of putative Bursera hybrids with morphologically intermediate leaf characteristics between B. bipinnata, B. cuneata, and B. palmeri. These Bursera species are mostly allopatric with overlapping ranges in the northeastern Michoacán. Despite the observations that hybridization is considered frequent in Mexican Bursera [16], there is no molecular evidence confirming this phenomenon, while no information exists on the presence of hybrid zones between co-occurring Bursera species.

In this study we examined putative hybrid populations using nuclear and chloroplast DNA sequences coupled with ecological niche modeling to ask the following questions: (i) Are morphologically intermediate individuals of hybrid origin from B. cuneata, B. palmeri and B. bipinnata? (ii) To which parental Bursera species do putative hybrids show the closest genetic affinity? (iii) Do the putative hybrids occur in a similar ecological niche to the parental species? (iv) Can we spatially and environmentally predict the formation of the hybrid zone through ecological niche modeling? (v) What is the role of past and future climatic conditions on the formation and maintenance of this hybrid zone?

Material and methods

Target species and sampling

Bursera palmeri, B. cuneata and B. bipinnata are dioecious and deciduous trees in the Bullockia section which is characterized by having well-developed cataphylls, bilocular ovaries, fruits with two valves, and smooth, grey, not exfoliating barks [16]. There is no information on the chromosome number in these species, but they are expected to be diploid as in other Bursera [35].

Leaf shape can differentiate the putative hybrids from their potential parental species: B. bipinnata has small bipinnate leaves, B. cuneata has larger, oblong to lanceolate, once-pinnate leaves of leathery texture and rugose appearance in the upper side. Bursera palmeri and B. cuneata are morphologically more similar, but B. palmeri has leaves with membranous texture and leaflets not rough on the upper side. The putative hybrids exhibit intermediate leaf characteristics including partial bipinnate leaves with large variation in size, texture, and level of leaf divisions (Fig 1). The bipinnate characteristic suggests that B. bipinnata is one of the parental species. We have not identified intermediate individuals between B. cuneata and B. palmeri.

[Figure omitted. See PDF.]

Fig 1. Leaf characteristics of (A) Bursera bipinnata, (B) B. cuneata, (C) B. palmeri and (D) putative hybrids, and (E) map of sampling localities.

Numbers within circles represent the number of samples collected per each species and the putative hybrids per locality; color designations: B. bipinnata (green), B. cuneata (blue), B. palmeri (yellow), and putative hybrids (pink). Sampling localities: Aguascalientes AGS, Jalisco JAL, Churintzio CH, Zacapu ZA, Pátzcuaro PTZ, Tarímbaro TM, Acumbaro AC, Ciudad de Mexico CDMX, Oaxaca, OAX.

https://doi.org/10.1371/journal.pone.0260382.g001

According to Rzedowski and Guevara-Féfer [19] and herbaria specimens deposited in the IEB herbarium (Instituto de Ecología A.C., Mexico) we identified various putative hybrid populations. During May 2017 to November 2019, we randomly collected fresh leaf samples from adult or juvenile trees of the three Bursera species and their putative hybrids in nine localities: four representing hybrid populations as identified by the co-occurrence of putative hybrid trees and their potential parental species and four populations of the parental species where neither of the other two Bursera species were present, thus representing “pure” parental populations (Table 1; Fig 1). Leaves were preserved in sealable plastic bags containing silica gel until DNA extractions were performed. GPS coordinates were recorded for each locality.

[Figure omitted. See PDF.]

Table 1. Sampling localities for three Bursera species and their putative hybrids in Mexico.

https://doi.org/10.1371/journal.pone.0260382.t001

Laboratory procedures

Genomic DNA from 20 mg of dried tissue was extracted following the CTAB extraction protocol of Doyle and Doyle [36]. We amplified four chloroplast (cpDNA) intergenic spacers (psbA-trnH, rbcL, trnK-matK, rps16) and three nuclear regions including the external transcribed spacer ETS, fourth intron of the phosphoenolpyruvate carboxylase PEPC, and third intron of the NIA-i3. These regions have been used in phylogenetic and hybridization studies of Bursera (e.g., [14, 21]). PCR reactions were performed using the MyTaq™ DNA polymerase kit (BIOLINE, London, United Kingdom) as follows: 5 μL of PCR buffer, 0.4 μL (10 μM) of each forward and reverse primer, 0.2 μL (5 U/ μL) of Taq polymerase and 5–10 ng of genomic DNA in a final volume of 25 μL. PCR conditions were as follows: 5 min denaturation at 95 °C, 35 cycles of 60 s denaturation at 95 °C, 60 s annealing at 53–56 °C, 90 s extension at 72 °C, and a final extension of 10 min at 72 °C. PCR products were run on a 1% agarose gel, stained with GelRed™ (BIOTIUM, Fremont, California) and visualized under UV light. Unpurified, amplified PCR products were sent for sequencing to MAGROGEN Inc. (Seoul, South Korea). We included a positive and negative control to check for contamination. Quality of chromatogram files were checked and edited using CHROMAS v2.6.5 (Technelysium Pty Ltd); unambiguous nucleotide sites for the nuclear regions were coded using the IUPAC nomenclature. Two cpDNA regions (psbA-trnH, rbcL) and NIA-i3 consistently failed to amplify, while the other two cpDNA regions (trnK-matK, rps16) resulted non-polymorphic in a subset of DNA samples for each species (GENBANK accessions: trnK-matK: OK413940—OK413942, rps16: OK438213—OK438215). Thus, our results are based on the polymorphism of the ETS and PEPC genes. We successfully sequenced 101 individuals, and these sequences are available in NCBI GENBANK (www.ncbi.nlm.nih.gov) with accession numbers: (1) ETS: B. cuneata (OK413888—OK413898), B. palmeri (OK413914—OK413917), B. bipinnata (OK413900—OK413913); (2) PEPC: B. cuneata (OK413918—OK413925), B. palmeri (OK413927—OK413930), B. bipinnata (OK413931—OK413939).

Recombination, phylogenetic relationships, and haplotype networks

Sequences were aligned using MUSCLE v3.8.31 [37] with default parameters and manually adjusted in MEGA X [38]. The gametic phase of heterozygous individuals was determine using PHASE 2.1 [39] as implemented in DNASP v6.12.03 [40]. The algorithm was run with 1000 iterations, 20% burn-in, a thinning interval of 10, and 0.9 as the minimum posterior probability of haplotypes. We performed all subsequent analysis with the phased sequences.

Evidence of recombination, which would be expected from interspecific hybridization, was evaluated with the pairwise homoplasy index (PHI-test) as implemented in SplitsTree v4.17.1 [41]. Moreover, we ran seven different detection recombination methods (RDP, GENECONV, BootScan, MaxChi, Chimaera, SiScan, and 3Seq) in the software RDP v4.101 [42]. We used the default parameters and considered the recombination events that were supported for three out of the seven selected methods.

We reconstructed phylogenetic relationships with the concatenated matrix using a Bayesian (BI) method. We used as outgroup sequences of other Bursera species retrieved from NCBI GenBank: B. simaruba (accession number: AF445946, GQ377992), B. kerberi (accession number: JF919125, AF445937) and B. fagaroides (accession number: AY309365, AY964604). We also sequenced a specimen of B. penicillata from Chapala Jalisco (Section Bullockia; accession number: OK413926, OK413899). The best model of sequence evolution was selected using the Akaike information criterion (AIC) in JMODELTEST v2.1.10 [43]. BI analyses were carried out using BEAST v2.5 [44]. We used BEAUti v1.8.1 [45] to adjust a substitution rate GTR gamma with invariant sites, a log normal relaxed molecular clock model to estimate the substitution rate. The analysis was run for 100 million generations for three independent runs. Posterior values were sampled every 1,000 generations and a burn-in of 10%. The BEAST output file was evaluated using TRACER v1.7 [46]. The less likely trees (probability > 0.8) were discarded using TreeAnnotator v1.8.2 [47]. The consensus tree was visualized with FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

We then constructed an unrooted phylogenetic network with the neighbor-net algorithm implemented in SplitsTree v4.17.1 [41] using a concatenated matrix of unique haplotypes for each species and putative hybrids. The neighbor-net algorithm is optimal for visualizing reticulated relationships among sequences, as it is expected from interspecific hybridization [41]. Additionally, haplotype relationships were observed with the median-joining (MJ) algorithm, and excluding gaps, as implemented in POPART [48]. We constructed a MJ network for each marker separately to visualize if the most frequent haplotypes for each species were shared with the putative hybrids irrespective of the nuclear marker.

Genetic diversity and differentiation

For each species and the putative hybrids, we estimated the number of haplotyes (h), haplotype diversity (Hd), nucleotide diversity (π), and expected heterozygosity (Het) using ARLEQUIN v3.5.2.2 [49]. We used BAPS v5.3 [50] to identify the most likely number of genetic clusters and genetic assignments between the potential parental species and their putative hybrids. We used the clustering of linked molecular data option, which is appropriate for diploid sequence data. First, we performed few exploratory runs from K = 2 to 15 and few replicates for each K (n = 5) to determine the most likely genetic clusters. After these initial searches, we found a likely number between 5 and 7. Thus, we carried out three independent runs from K = 3 to 9 with 10 replicates each. The most likely number of K clusters was determined by its higher posterior probability and likelihood value. Then, we used the binary output file of the mixture clustering to run an admixture analysis. We used 1000 iterations, 5 reference individuals for each population and 100 iterations for the reference individuals. These parameters were defined by performing several exploratory runs and checking the consistency among runs. To assess patterns of genetic differentiation among the three Bursera species and the putative hybrids, we carried out a Discriminant Analysis of Principal Components (DAPC) in the R library adegenet [51] and used a binary coded matrix of presence (1) and absence (0) of haplotypes for each marker per individual. The DAPC is a multivariate method free of Hardy-Weinberg and linkage disequilibrium assumptions [51]. We used the three species and the putative hybrids a priori within the dapc function, while the optimal number of PCs to retain was optimized using the function xvalDapc [51]. Pairwise FST comparisons between the four groups were calculated with 1000 permutations to test for statistically significant differences in ARLEQUIN v3.5.2.2 [49].

Ecological niche modeling and niche overlap zone

We compiled herbarium records for the three Bursera species and the putative hybrids from different web databases, including the Global Biodiversity Information Facility (GBIF; http://www.gbif.org/), SEINet (https://swbiodiversity.org/seinet/), Missouri Botanical Garden (http://www.tropicos.org), National Herbarium of Mexico (MEXU; http://www.ib.unam.mx/botanica/herbario/), and the IEB Herbarium in Mexico. In addition, we requested Bursera herbarium records, including putative hybrids, from 33 Mexican universities and institutional herbaria (S1 Table). Specimen photography, collector and identifier were consulted to corroborate the veracity of specimens. Our search retrieved 1,529 records, which included the entire distribution range of the three parental Bursera species and 19 records of putative hybrids. To reduce sampling biases for the niche modeling, the occurrences were spatially thinned at 20 km for B. bipinnata and B. palmeri, and 15 km for B. cuneata. The variable radii were determined based on the number of available records and the spatial extent of each species distribution. After this filtering, the dataset comprised 313 occurrences for B. bipinnata, 105 for B. palmeri, and 93 for B. cuneata. Each species occurrences were equally split in two datasets, one for building the niche model and the second for data testing. No spatial thinning was applied for the putative hybrid records. These records were used only to carry out the similarity and equivalency test (see below).

We selected 11 bioclimatic layers in 30 arc-seconds from WorldClim (https://worldclim.org/data/worldclim21.html) [52]: Annual mean temperature (BIO1), mean diurnal range (BIO2), isothermality (BIO3), temperature seasonality (BIO4), temperature annual range (BIO7), annual precipitation (BIO12), precipitation of the wettest month (BIO13), precipitation of the driest month (BIO14), precipitation seasonality (BIO15), precipitation of the warmest quarter (BIO18), precipitation the coldest quarter (BIO19). This layer selection was based on other niche model and distributional studies in Bursera species [53, 54]. Also, Bursera occurrence at finer spatial scales is determined by slope aspect and roughness of the ground [55], thus we generated a topographic layer based on a SRTM digital elevation model (https://www2.jpl.nasa.gov/srtm/). The calibration areas were delimited using the WWF biogeographic regions to delimit the species accessible areas [56]. We discarded highly correlated variables (≥ 0.8) using Spearman´s rank correlation coefficients and removed variables with the lowest contribution in the final model measured by Jackknife analysis using the R library raster [57].

Niche models under current climatic conditions using the above-mentioned variables were constructed using the maximum entropy algorithm implemented in MAXENT v3.4.1 [58] using the R library kuenm v1.1.1 [59]. Levels of model complexity were evaluated, such as over-fitting by varying the regularization multiplier (RM) (0.1, 0.5, 1, 2, 3, 4) and feature classes linear (L), quadratic (Q), product (P), threshold (T), and hinge (H) in 5 combinations (i.e., L, LQ, LQP, LQPT, LQPTH), and 10000 background points for B. bipinnata and B. palmeri and 7000 points for B. cuneata due to its smaller range. This resulted in 150 candidate niche models for each parental species. The niche models were evaluated based on the statistical significance of the partial Receiver Operating Characteristic (pROC) [60], 10000 random points, 500 iterations, and 0.5% omission rates. We generated several metrics to rank the best models: AUC, omission rate, AICc, delta AICc, parameter number, and WAICc. The niche model selected with the best parameters was constructed and projected in the same polygon with ten bootstrap iterations. The best niche model for each Bursera species were set to thresholds with the 10 percentile training presence values to produce binary raster maps. These binary maps were overlapped and the areas of conjunction of pixels among the three parental species were defined as the niche overlap zone. This procedure was carried out using QGIS v2.12 [61]. The 19 putative hybrid records were mapped on the overlapping niche zone to determine their spatial correspondence.

Niche equivalency and similarity tests

To test for ecological niche divergence as evidence for reproductive isolation [62] among parental Bursera species and putative hybrids, we used the approach based on principal component analyses of climatic and topographic variables (PCA-env) to perform multivariate comparisons of niche overlap between pairs of species and the putative hybrids [63]. PCA-env uses a kernel density function to compute the density of occurrences in the multivariate PCA space, to account for potential bias from unequal sampling effort. This procedure consists in three steps: (i) calculation of the density of occurrences and environmental factors across the axes of the PCA-env, (ii) evaluation of the niche superposition along the gradient of multivariate analysis, and (iii) comparisons with repetitions of randomly generated simulated values of the Schoener’s D metric [64] and the modified Hellinger distance I [65]. Based on ecological niche modeling of Bursera species, we performed the PCA with the most common and important variables according to the permutation importance. We generated four environmental axes (PCA-env) and selected the two with the highest variance explained. With the selected two PCA-env axes, we extracted the environmental values of the occurrence points and the calibration areas for each species, the putative hybrids, and for the total occurrence points, i.e., the four groups, in the study area. The overlap niche zone was used as the hybrid calibration area. With this, we constructed the species object of each parental Bursera species and the putative hybrids, which represents the environmental tolerance threshold of each species. The species objects were used to compute the D and I metrics, which calculate the niche overlap ranging from 0 (no overlap) to 1 (complete overlap). In this case, the niche overlap was calculated between parental species and between a parental species and the putative hybrids. Then, using the two metrics, we performed niche equivalency and similarity tests that compare the niche overlap values to a null distribution. The equivalency test evaluates whether two niches are significantly different from each other (not environmentally identical), while the niche similarity test evaluates whether two niches are more similar (niche A can predict niche B and vice versa) than expected by chance [66]. We used the option “alternative” to test for niche conservatism, i.e., greater, in the R library ecospat [63] by performing 5,000 permutations and 100 replicates.

Paleodistributions and future climatic projections

To reconstruct the historical dynamics of the potential hybrid zone, i.e., overlap niche zone, we projected the ecological niche model under current climatic conditions into three past scenarios: Last Interglacial (LIG ~ 120–140,000 ya) at 30 seconds [67] downloaded from WorldClim (http://www.worldclim.com), Last Glacial Maximum (LGM ~ 21,000 ya) and late Holocene (LH ~ 2,200 ya) at 2.5 min resolution and using the global models MIROC-ESM (Model for Interdisciplinary Research on Climate) [68] and CCSM4 (The Community Climate System Model, https://www.cesm.ucar.edu/models/ccsm4.0/) [69] obtained from PaleoClim (http://www.paleoclim.org/). To project future climatic scenarios of the hybrid zone, we projected niche models at the years 2050 and 2070. In both cases we employed the global models CCSM 4.0 and the MIROC5 based on two scenarios of Representative Concentration Pathway (RCP) referred as RCP4.5 (optimistic) and RCP8.5 (pessimistic) [70]. The RCP 8.5 is the scenario with the highest predicted greenhouse gas emissions compared with RCP 4.5 [71, 72]. All projections were set to thresholds with 10 percentile training presence values to produce raster binary maps. The binary maps of each species were overlapped for each time scenario and the areas of conjunction of pixels among the three parental species were defined as the niche overlap zone for each time. This procedure was performed in QGIS v2.12 [61]. Lastly, to assess range shifts in the niche overlap zone in response to the climatic changes, we extracted elevation values from 500 random points from the binary maps of the different time niche scenarios for the CCSM4 and MIROC5. We constructed density plots with the 500 random points to observe the altitudinal range differences between scenarios for the overlap niche zone.

Results

Recombination, phylogenetic and haplotype relationships

The concatenated ETS (421 bp) and PEPC (299 bp) aligned matrix included 202 phased sequences and had a total of sites excluding gaps or missing data of 720 bp with 38 parsimony informative sites. The PHI-test of recombination was statistically significant (P < 0.027). RDP found evidence of two recombination events, but only one of them was supported for three different methods (MaxChi, SiScan, 3Seq, start-end position 25–328, P = 0.001). This recombination event involved an individual from B. cuneata, whose major parent was a putative hybrid, and a minor parent was an individual of B. cuneata. These three individuals were from Acumbaro, Michoacán (AC).

The best-fit model of sequence evolution was HKY gamma with invariable sites for the ETS and GTR gamma with invariant sites for the PEPC. The BI phylogenetic tree resulted in a high-supported Bullockia clade (posterior probability = 1), but the three parental species were not monophyletic. Putative hybrid sequences were intermixed among subclades of the three parental species in high proportion and as sister subclades of B. cuneata and B. bipinnata (S1 Fig). Results from the neighbor-net in SplitsTree showed a highly reticulated tree, where putative hybrid and parental haplotypes were intermixed among major clades. Like the BI tree, the three Bursera species were not monophyletic, of which Bursera cuneata haplotypes were the most intermixed with putative hybrids (Fig 2A).

[Figure omitted. See PDF.]

Fig 2. Neighbor-net tree (A) and median-joining haplotype networks for 202 phased sequences of the (B) ETS and (C) PEPC nuclear genes in Bursera bipinnata (green), B. cuneata (blue), B. palmeri (yellow), and putative hybrids (pink).

Size of nodes is proportional to the number of sequences per haplotype, dashed lines represent the number of mutational steps separating haplotypes, and black nodes represent unsampled haplotypes.

https://doi.org/10.1371/journal.pone.0260382.g002

The MJ haplotype network for the ETS gene showed 43 haplotypes, seven exclusive to B. cuneata, three to B. palmeri, nine to B. bipinnata and 16 to the putative hybrids. The three most frequent haplotypes in each of the three Bursera species were shared with the putative hybrids. For B. palmeri none of the haplotypes were shared with any of the other two parental species, while B. cuneata and B. bipinnata shared two of the most frequent haplotypes (Fig 2B). The MJ haplotype network for the PEPC gene showed 26 haplotypes, four exclusive to B. cuneata, seven to B. bipinnata, one to B. palmeri, and eight to the putative hybrids. A large proportion of the putative hybrid sequences were shared among the most common haplotypes of the three Bursera species. These common haplotypes were shared among the parental species, although in lower proportions relative to the putative hybrids (Fig 2C).

Genetic diversity and differentiation.

Estimates of genetic diversity showed moderate to high diversity values (Table 2). The putative hybrids exhibited the highest genetic diversity, including expected heterozygosity, closely followed by B. bipinnata. Of the three Bursera species, B. palmeri showed the lowest diversity values, which haplotype diversity was four-fold lower than the putative hybrids.

[Figure omitted. See PDF.]

Table 2. Genetic diversity indices of the concatenated ETS and PEPC nuclear genes in three Bursera species and their putative hybrids.

https://doi.org/10.1371/journal.pone.0260382.t002

Results from BAPS showed that the most likely number of genetic clusters was K = 7 (log marginal likelihood = -2633.06, PP = 0.99). Admixture analysis showed that the putative hybrids showed a high proportion of individuals assigned to each of the three parental Bursera species, but also individuals of each Bursera species were assigned to the putative hybrids, of which B. bipinnata showed the larger number of individuals assigned to the putative hybrids (Fig 3A). Few individuals showed admixed ancestry between two or three clusters, but all of them shared ancestry with the putative hybrids (Fig 3A). The DAPC plot showed the putative hybrids in-between individuals of the three parental species, but with closer affinity to B. cuneata and to a lesser extent to B. palmeri. The DA 1 axis showed the separation of B. palmeri from the rest, while DA 2 axis showed a larger separation of B. bipinnata, although there was a considerable overlap with the putative hybrids (Fig 3B). Pairwise FST genetic distances showed that B. cuneata and the putative hybrids were genetically closer relative to the other two species, while B. palmeri was the most distant to the putative hybrids (S2 Table).

[Figure omitted. See PDF.]

Fig 3. (A) BAPS barplot of the best-supported K = 7 clusters and (B) DAPC plot of haplotypes in Bursera bipinnata (Bp), B. cuneata (Bc), B. palmeri (Bm), and the putative hybrids (Bh).

The insets showed the proportion of variance explained by each of the two DA eigenvalues (DA 1 & 2 respectively) and the number of retained PCA eingenvalues.

https://doi.org/10.1371/journal.pone.0260382.g003

Ecological niche modeling and overlap niche analyses

According to the permutation values, the most important variables among the three Bursera species were temperature seasonality (BIO4), precipitation seasonality (BIO15), annual mean temperature (BIO1), and precipitation of the wettest quarter (BIO13). The AUC values for the niche models were high for the three species: B. bipinnata, 0.914 ± 0.03, B. palmeri, 0.938 ± 0.007 and B. cuneata, 0.900 ± 0.009, indicating that occurrence points were strongly differentiated from background locations, so model distributions were not random.

The B. bipinnata ecological niche model exhibited high suitability along the mountains of the Pacific coast from Nayarit up to Guatemala, and largely comprising the Bajío and the Balsas regions, and with altitudinal range between 900 to 1,780 m.a.s.l. For B. palmeri the ecological niche model extended in the southwestern Mexican Plateau and few areas of the Sierra Madre del Sur and with altitudinal range between 1,175 to 2,050 m.a.s.l. In contrast, the niche of B. cuneata was smaller and projected in the central portion of the Trans-Mexican Volcanic Belt and some areas of the Sierra Madre del Sur and with higher altitudinal range between 1,750 to 2,400 m.a.s.l. (S2 Fig). The predicted niche overlap zone among the three species showed a sympatric area in the Bajío region, largely between the states of Michoacán and Guanajuato, with few small and fragmented areas across the Sierra Madre del Sur. The total extent of the overlap zone was 74,023 km2, and accurately predicted the presence of the putative hybrids (Fig 4). The PCA-env was generated with the most important and common variables. The first axis explained 54% of the environmental variation (BIO1, BIO4 and BIO13) and the second axis explained the 25.5% of variation (BIO15) (S3 Fig)

[Figure omitted. See PDF.]

Fig 4. Ecological niche projections of the overlap zone among Bursera bipinnata, B. cuneata, and B. palmeri at the present and past scenarios.

(A) Hybrid leaf sample, (B) present time, (C) late Holocene (LH, CCSM4), (D) late Holocene (LH MIROC), (E) Last Glacial Maximum (LGM, CCSM), (F) Last Glacial Maximum (LGM, MIROC) and (G) Last Interglacial (LIG). Red areas on the map indicate the predicted overlap niche zones, purple dots represent the occurrence of hybrid populations, numbers in km represent the total extension of the hybrid zone.

https://doi.org/10.1371/journal.pone.0260382.g004

The similarity test showed that the B. palmeri niche significantly predicted the niche of B. cuneata (D = 0.45, I = 0.67, P < 0.05), while the other species pairwise comparisons were not significant. Moreover, the B. palmeri niche significantly predicted the niche of the putative hybrids (D = 0.13, I = 0.36 P < 0.05), but B. cuneata predicted the niche of the putative hybrids with higher overlap values (D = 0.31, I = 0.54, P < 0.05). Neither B. bipinnata nor the putative hybrids could predict the niche of each other (D = 0.04, I = 0.20, P > 0.05) (Table 3, S3 Fig). The equivalency test indicated that the niches of B. bipinnata and B. cuneata (D = 0.26, I = 0.45, P < 0.05), and B. cuneata and B. palmeri were not significantly equivalent (i.e., divergent) (D = 0.45, I = 0.67, P < 0.05), while the niches of B. bipinnata and B. palmeri were equivalent (D = 0.46, I = 0.62, P > 0.05). The comparisons with respect to the putative hybrids and the parental species showed that the putative hybrids and B. bipinnata had not equivalent niches (D = 0.04, I = 0.20, P < 0.05), but the putative hybrids had equivalent niches with B. palmeri (D = 0.13, I = 0.36, P > 0.05) and B. cuneata (D = 0.31, I = 0.54, P > 0.05).

[Figure omitted. See PDF.]

Table 3. Pairwise niche overlap values of Schoener’s D and Hellinger’s I between the three Bursera species and their putative hybrids.

https://doi.org/10.1371/journal.pone.0260382.t003

Temporal projections of the hybrid zone

The ecological niche model during the warm period of the LIG (~ 130 kya) did not predict a niche overlap zone among the three Bursera species, but during the cooling conditions of the LGM (~ 21 kya) the CCSM4 and MIROC projections showed a niche overlap towards some disconnected areas in the south coast of the Pacific and the Sierra Madre del Sur. The MIROC model predicted a larger area of niche overlap (102,471 km2) relative to the CCSM4 model (87,591 km2). These two overlap niche areas did not coincide with presence of the putative hybrid records at the present. With the increase of temperature during the LH (~ 2.2 kya) under the MIROC model, the overlap niche zone shifted north towards the Bajío region and Sierra Madre del Sur, whereas for the CCSM4 model it was limited to the east of the Sierra Madre del Sur. For both models, the niche overlap zone decreased in extension relative to the LGM. For the LH, the putative hybrid records completely match the overlap niche zone predicted under the MIROC, and partially match under the CCSM4 (Fig 4).

For the future scenarios, the overlap niche zone was located mainly in the Bajío region. For CCSM4 global models, the predictions of the overlap zones were larger than for the MIROC models. For 2050, the CCSM4 scenario showed a relatively smaller overlap niche zone under the optimistic scenario (93,371 km2, RCP 4.5), and a larger area under the pessimistic scenario (104,668 km2, RCP 8.5). For 2070, the pessimistic and optimistic CCSM4 scenarios predicted very similar locations and extensions of the niche overlap zones. All 2050 and 2070 scenarios coincided with the presence of the putative hybrid records. In contrast to the CCSM4, the MIROC models predicted a decrease of the overlap niche zone from 2050 to 2070, although the location coincided with the CCSM4 models (Fig 5). According to the altitudinal density plots, a range shift of the niche overlap zone through time was evident: in the past (since the LGM) the hybrid zone occurred in the lowlands, whereas at the present and the future, the hybrid zone is shifting to greater elevational ranges (S4 and S5 Figs).

[Figure omitted. See PDF.]

Fig 5. Ecological niche predictions of the overlap niche zone among Bursera bipinnata, B. cuneata, and B. palmeri for future climatic scenarios: (A) 2050 optimistic scenario (RCP 4.5, CCSM4), (B) 2050 optimistic scenario (RCP 4.5, MIROC), (C) 2050 pessimistic scenario (RCP 8.5, CCSM4), (D) 2050 pessimistic scenario (RCP 8.5, MIROC), (E) 2070 optimistic scenario (RCP 4.5, MIROC), (F) 2070 pessimistic scenario (RCP 8.5, CCSM4), (G) 2070 pessimistic scenario (RCP 8.5, CCSM4), (H) 2070 pessimistic scenario (RCP 8.5, MIROC).

Red areas on the map indicate the predicted overlap niche zones, purple dots represent the occurrence of hybrid populations, and numbers in km represent the total extension of the hybrid zone.

https://doi.org/10.1371/journal.pone.0260382.g005

Discussion

In this study, analysis of the ETS and PEPC nuclear sequences confirmed the hybrid origin of morphologically intermediate individuals in five populations with the co-occurrence of B. bipinnata, B. palmeri and B. cuneata in the Bajío region in Mexico. Results showed that the putative hybrids shared genetic ancestry from the three Bursera species, but with closer genetic affinity to B. cuneata and B. bipinnata. Ecological niche modeling accurately predicted the occurrence of known putative hybrid records and showed that the overlap zone among the three Bursera species, where potential hybridization may occur, extends in a larger area than previously thought. The putative hybrids did not show ecological divergence from B. palmeri and B. cuneata but did from B. bipinnata, while B. palmeri significantly predicted the niche of B. cuneata. Past and future ecological niche projections showed that the overlap zone likely occurred since the LGM (~ 21 kya) and have ascended altitudinally through time.

Genetic evidence of interspecific hybridization in Bursera

Presence of morphologically intermediate individuals between Bursera species has been reported extensively from Baja California [20], the Bajío region [16, 19], the western Balsas [22], and to the Tehuacán-Cuicatlán region in Mexico [16, 73]. In the Bajío, Rzedowski and Guevara-Féfer [19] identified putative hybrids derived from B. bipinnata with B. cuneata and B. palmeri. Our findings from phylogenetic reconstructions revealed that these putative hybrids were intermixed within subclades of the three Bursera species. A similar pattern was observed by Weeks and Simpson [74] among hybridizing Hispaniolan Bursera. Moreover, the shared origin of the putative hybrids from the three Bursera species was clearer from the ETS and PEPC MJ haplotype networks, which showed that most common species haplotypes were shared in large proportion with the putative hybrids.

Interspecific hybridization is expected to be common in closely related species that have not developed strong reproductive isolation. For the Bullockia section, two natural groups were recognized based on the fruit and flower characteristics [22] and further confirmed with molecular data [74, 75]. Taxonomic identity of the three Bursera species here has been confirmed with morphological and molecular data [14, 16, 19]. According to the pairwise FST estimates, genetic differentiation between B. bipinnata and B. cuneata (FST = 0.47) is lower than between B. bipinnata and B. palmeri (FST = 0.52). Likewise, the putative hybrids were more differentiated to B. palmeri relative to the other two species. This pattern was also evident from the DAPC plot, in which the putative hybrids were largely overlapped with B. cuneata and B. bipinnata. Our field observations and genetic data suggest that most of the putative hybrids were likely the result of interbreeding between B. bipinnata and B. cuneata, and to a lesser extent to B. bipinnata and B. palmeri, and then to B. cuneata and B. palmeri. Of the three species, B. bipinnata is the species that is known to hybridize with a larger number of Bursera species across its range [16].

Nothing is known about the development of reproductive barriers in Bursera which could help explain the propensity of some species to hybridize more frequently than others. Hybrid infertility is expected to limit the occurrence of backcrosses acting as a postzygotic reproductive barrier [76]. A recent study showed that pollen of these putative hybrids is fertile [77], which may facilitate the occurrence of backcrosses and F2 hybrids. Results from BAPS revealed genetic admixture between putative hybrids and the three Bursera species, but it was noticeably the higher genetic affinity of some B. bipinnata and B. cuneata individuals to the putative hybrids. This result implies that patterns of hybridization among these Bursera species are more complex than previously thought. This observed genetic pattern may reflect backcrossing events or F2 hybrid generations in which morphologically intermediacy is no longer recognized. Although the ETS and PEPC sequences have proven to be useful in phylogenetic and hybridization studies in Bursera including the present study [21, 25], with these molecular markers is not feasible to unambiguously characterize F2 hybrids and backcrosses scenarios. Moreover, we cannot completely rule out the possibility that shared haplotypes between parental species and putative hybrids were the result of incomplete lineage sorting or ancestral character retention [78].

Our estimates of genetic diversity showed that the putative hybrids had the highest levels of genetic diversity, including heterozygosity, expected from their mixed origins. Field observations have noted that these putative hybrids can form vigorous trees with abundant fruit and flower production [19], but there is no specific data on their local fitness. Hybridization is an important evolutionary force that can increase genetic diversity of the populations involved [12, 79]. If backcrosses are possible, hybridization can influence levels of genetic diversity and structure in parental species with consequences for local adaptation [80]. Further studies with a much higher number of unlinked molecular markers, i.e., single nucleotide polymorphism SNPs, in a large number of populations across the potential hybrid zone are needed to better characterize complex hybridization scenarios, including F1, F2 hybrids and of introgression events, and to understand the role of interspecific hybridization in population genetic composition and adaptive introgression.

The overlap niche zone and ecological divergence

Our results from ecological niche modeling showed a wide overlap niche zone among the three Bursera species that extends across the Bajío region comprising mostly the states of Jalisco, Michoacán, Guanajuato, and Querétaro. This ecological overlap zone accurately predicted the occurrence of known putative hybrid records, which provides evidence to suggest that it corresponds to a potential hybrid zone [10]. We found that the most important environmental variables associated with the potential hybrid zone were temperature and precipitation, the ecological factors that have been associated to the distribution of Bursera species in the TDF [53, 81]. This potential hybrid zone largely coincides with the habitat of B. palmeri and B. cuneata. It is important to emphasize that hybridization rates and hybridization success likely vary across species ranges where they co-occur [82, 83]. Our model predictions on the extension of the overlap niche zone should be interpreted cautiously and not as strong evidence of hybrid presence as few hybrid records exist. Hence, it would be key to verify and corroborate the presence of hybrids at other sites across the potential hybrid zone where Bursera hybrids have not been reported, but which are predicted to exist.

The similarity test showed that B. cuneata had an ecological niche similar to B. palmeri, but the equivalence test showed that these species had divergent niches, while B. bipinnata only showed ecological divergence with B. cuneata. Among all Bursera species, B. cuneata occurs at higher altitudinal ranges (1850 to 2500 m.a.s.l.), usually associated with transition zones between the TDF and pine-oak forests [19]. These results suggest that ecological divergence might be a mechanism promoting reproductive isolation between B. cuneata and B. bipinnata which are likely the species that hybridize more frequently, and the most genetically similar. Ecological divergence has been suggested to contribute to the establishment of reproductive barriers between hybridizing species [84, 85] for several taxa including insects [27], reptiles [86, 87], mammals [88], and invasive plants [89].

Ecological divergence of the hybrids with respect to the parental species can be a signature of their ecological adaptation across hybrid zones [5] and their potential to facilitate the development of novel adaptations to changing environments [90, 91]. We found that putative hybrids showed ecological divergence to B. bipinnata, but not to B. cuneata and B. palmeri, despite that B. palmeri showed the highest genetic differentiation to the putative hybrids. The putative hybrids exhibited the highest niche similarity with B. cuneata. The absence of ecological divergence of the putative hybrids with B. cuneata further confirms their high relationship and suggests that hybrids might not form and establish independent populations from this parental species (which we have observed in the field). However, signatures of ecological divergence might not be detected at this regional scale. Both endogenous and exogenous sources of selection on hybrids can vary substantially across a hybrid zone which can differentially influence hybridization and introgression rates and hybrid fitness at large and local scales [92–94]. For plants, it has been observed that hybridization rates can be influenced by pollinator behavior [95], reproductive phenology [96], ecological disturbance [97], and elevational ranges [92]. The patchy distribution of the TDF and the high landscape heterogeneity in which these Bursera species currently persist [34], may impose varying selective forces at short spatial scales with implications in local fitness. Future research is needed to explore the role of the ecological context on hybridization rates and hybrid fitness by increasing the number of populations analyzed across the potential hybrid zone.

Temporal dynamics of the hybrid zone

The temporal overlap niche projections using past bioclimatic variables predicted the occurrence of a potential hybrid zone during the Last Glacial Maximum (~ 21 kya), although its distribution was largely located in southern Mexico in the Sierra Madre del Sur. Later during the warm conditions of the late Holocene (~ 2.3 kya), this potential hybrid zone shifted northward, but it contracted in extension. In this period, both models predicted the presence of the current hybrid records. In Mexico, during the interglacial and glacial climatic cycles of the Pleistocene, forest species experienced repeated expansions and contractions in response to climatic changes [98]. There is scarce information on the historical temporal changes of TDF species during the Pleistocene in Mexico. Based on ecological niche modeling for areas of endemism in Bursera, Gámez et al. [54] found elevational changes of Bursera species from the LGM to the present in response to climate warming and that may be related to limited historical dispersal [14]. The potential hybrid zone location during the LGM for the MIROC model corresponds to an area known as the Western Balsas, while for the CCSM model, the potential hybrid zone was predicted in the Eastern Balsas. Both areas were suggested as refugia during the cooler and wetter conditions of the LGM [54]. From the late Holocene to the present there was an increase in the extension of the potential hybrid zone. Overall, future scenario projections under the CCSM model showed a potential increment of the suitability for the hybrid zone for 2050 and 2070 relative to the present. For both circulation models, the predictions indicated the increase in elevational ranges, but with a larger extent under the CCSM4 model relative to MIROC, which spatial extent decreases by 2070. The increase in elevational ranges coincided with the observed trends of the TDF to ascend altitudinally in response to climatic changes [54, 99] which likely will create the conditions for other Bursera species to come into contact. This may increase the opportunities for interspecific gene exchange as it has been predicted for several taxa [11, 100] and with important consequences for TDF biodiversity including species losses and gains at local scales [54].

Conclusion

This study provides the first insight into existence of a potential hybrid zone among three Mexican Bursera species (Bullockia section) based on molecular evidence and ecological niche modeling. The examination of the temporal dynamics of the hybrid zone in response to climatic changes provided useful information in understanding the role of environmental factors in shaping the formation and maintenance of potential hybrid zones in this plant taxa with a high number of co-occurring endemic species.

Supporting information

S1 Fig. BEAST phylogenetic tree of the concatenated ETS and PEPC nuclear genes.

The tree shows the relationships between the tree Bursera species and the putative hybrids. Outgroup species were B. simaruba, B. kerberi, and B. fagaroides. Values above nodes represent posterior probabilities.

https://doi.org/10.1371/journal.pone.0260382.s001

(PDF)

S2 Fig. Predicted suitable areas for each Bursera species through ecological niche modeling.

The model was constructed based on the best calibration based on Kuenm. Green dots represent the training set and blue dots the testing set. The black and gray shading represent an altitudinal gradient with lighter areas being altitudinally higher.

https://doi.org/10.1371/journal.pone.0260382.s002

(PDF)

S3 Fig. Results from the similarity test between the putative hybrid and each of the parental Bursera species.

The left panel shows the environmental space of Bursera palmeri (orange), B. cuneata (blue) and B. bipinnata (green) and the putative hybrids (purple). The gray area denotes the overlap area. The right panel shows the distribution of observed vs expected values of Schoener’s D, the blue lines denote the range of 95% of expected values and the red lines the observed values.

https://doi.org/10.1371/journal.pone.0260382.s003

(PDF)

S4 Fig. Density plots for the CCSM4 models showing the altitudinal range values predicted for the overlap niche zone (i.e., hybrid zone) for each climatic niche scenario.

Scenarios: LGM Last Glacial Maximu, LH Late Holocene. Values of the x-axis denotes the elevational ranges.

https://doi.org/10.1371/journal.pone.0260382.s004

(PDF)

S5 Fig. Density plots for the MIROC models showing the altitudinal range values predicted for the overlap niche zone (i.e., hybrid zone) for each climatic niche scenario.

Scenarios: LGM Last Glacial Maximum, LH Late Holocene. Values of the x-axis denotes the elevational ranges.

https://doi.org/10.1371/journal.pone.0260382.s005

(PDF)

S1 Table. Species occurrence used for the niche modelling.

https://doi.org/10.1371/journal.pone.0260382.s006

(XLSX)

S2 Table. Pairwise FST (below) and GST (above) genetic distances among the three Bursera species and the putative hybrids.

All comparisons were statistically significant P = 0.0001.

https://doi.org/10.1371/journal.pone.0260382.s007

(PDF)

Acknowledgments

For assistance in the field and sample preparation we thank Benjamín Castillo Ponce, Bruno A. Gutiérrez Becerril, Stephanie Aguilera López, Tania Andrade Ortiz, and Víctor Reyes Pino, and to several field assistants from local communities. A sampling permit was granted by SEMARNAT SGPA/DGGFS/712/1062/18. We are thankful to two external reviewers and the editor for comments that helped us significantly improve the quality of our manuscript.

Citation: Quintero Melecio E, Rico Y, Lira Noriega A, González Rodríguez A (2021) Molecular evidence and ecological niche modeling reveal an extensive hybrid zone among three Bursera species (section Bullockia). PLoS ONE 16(11): e0260382. https://doi.org/10.1371/journal.pone.0260382

References

1. Soltis PS, Soltis DE. The role of hybridization in plant speciation. Annu Rev Plant Biol. 2009;60: 561–588. pmid:19575590

2. Grant V. Plant Speciation. Columbia University Press; 1981. https://doi.org/10.7312/gran92318

3. Harrison RG. Hybrid Zones and the Evolutionary Process. Harrison R.G., editor. Oxford University Press, Inc; 1993. https://books.google.com.mx/books?hl=en&lr=&id=aFJFkVKskYIC&oi=fnd&pg=PA3&dq=Harrison+hybridization&ots=MFdXknOfPH&sig=lHy52RktdkRydeFKSp-HnZShLBs#v=onepage&q=Harrisonhybridization&f=false

4. Widmer A, Lexer C, Cozzolino S. Evolution of reproductive isolation in plants. Heredity. Nature Publishing Group; 2009. pp. 31–38. pmid:18648386

5. Ortego J, Gugger PF, Riordan EC, Sork VL. Influence of climatic niche suitability and geographical overlap on hybridization patterns among southern Californian oaks. J Biogeogr. 2014;41: 1895–1908.

6. Abbott R, Albach D, Ansell S, Arntzen JW, Baird SJE, Bierne N, et al. Hybridization and speciation. Journal of Evolutionary Biology. 2013. pp. 229–246. pmid:23323997

7. Rhymer JM, Simberloff D. Extinction by hybridization and introgression. Annu Rev Ecol Syst. 1996;27: 83–109.

8. Todesco M, Pascual MA, Owens GL, Ostevik KL, Moyers BT, Hübner S, et al. Hybridization and extinction. Evol Appl. 2016;9: 892–908. pmid:27468307

9. Schierenbeck KA, Ellstrand NC. Hybridization and the evolution of invasiveness in plants and other organisms. Biol Invasions. 2009;11: 1093–1105.

10. Taylor SA, Larson EL, Harrison RG. Hybrid zones: Windows on climate change. Trends Ecol Evol. 2015;30: 398–406. pmid:25982153

11. Grabenstein KC, Taylor SA. Breaking Barriers: Causes, Consequences, and Experimental Utility of Human-Mediated Hybridization. Trends in Ecology and Evolution. Elsevier Ltd; 2018. pp. 198–212. pmid:29306562

12. López-Caamal A, Tovar-Sánchez E. Genetic, morphological, and chemical patterns of plant hybridization. Revista Chilena de Historia Natural. Sociedad de Biologia de Chile; 2014. pp. 1–14.

13. Abbott RJ. Plant speciation across environmental gradients and the occurrence and nature of hybrid zones. Journal of Systematics and Evolution. Wiley-Liss Inc.; 2017. pp. 238–258.

14. de-Nova JA, Medina R, Montero JC, Weeks A, Rosell JA, Olson ME, et al. Insights into the historical construction of species-rich Mesoamerican seasonally dry tropical forests: The diversification of Bursera (Burseraceae, Sapindales). New Phytol. 2012;193: 276–287. pmid:21955031

15. Rzedowski J. Familia Burseraceae. In: Ramamoorthy TP, Bye R, Lot A, Fa J, editors. Biological diversity of Mexico: Origins and distribution. Oxford, UK: Oxford University Press; 1988. pp. 129–144.

16. McVaugh R, Rzedowski J. Synopsis of the Genus Bursera L. in Western Mexico, with Notes on the Material of Bursera Collected by Sesse & Mocino. Kew Bull. 1965;18: 381.

17. López AM. Identidad y simbolismo del copal prehispánico y reciente. Arqueología. 2004. https://www.revistas.inah.gob.mx/index.php/arqueologia/article/view/6278

18. Bonfil-Sanders C, Cajero-Lázaro I, Evans RY. Germinación de semillas de seis especies de Bursera del centro de México. Agrociencia. 2008;42: 827–834.

19. Rzedowski J, Guevara-Fefér F. FAMILIA BURSERACEAE. Fascículo 3. Flora del Bajío y de Regiones Adyacentes; 1994. https://libros.inecol.mx/index.php/FB/catalog/book/129

20. Pérez-Navarro JJ. El género Bursera Jacq. ex. L. (Burseraceae) en la Península de Baja California. Centro de Investigaciones Biológicas del Occidente, S.C. 2001.

21. Weeks A, Simpson BB. Molecular genetic evidence for interspecific hybridization among endemic Hispaniolan Bursera (Burseraceae). Am J Bot. 2004;91: 976–984. pmid:21653453

22. Toledo-Manzur CA. El género Bursera (Burseraceae) en el estado de Guerrero (México). Universidad Nacional Autónoma de México. 1982.

23. Rzedowski J, Kruse H. Algunas tendecias evolutivas en Bursera (Burseraceae). Taxon. 1988;28: 103–116.

24. Rzedowski J, Ortiz E. Estudios quimiotaxonómicos de Bursera (Burseraceae). II. Una especie nueva de origen hibrido de la Barranca de Tolantongo, estado de Hidalgo. Acta Bot Mex. 1988;1: 11–19.

25. Weeks A, Tye A. Phylogeography of palo santo trees (Bursera graveolens and Bursera malacophylla; Burseraceae) in the Galápagos archipelago. Bot J Linn Soc. 2009;161: 396–410.

26. Rosell JA, Olson ME, Weeks A, De-Nova JA, Lemos RM, Camacho JP, et al. Diversification in species complexes: Tests of species origin and delimitation in the Bursera simaruba clade of tropical trees (Burseraceae). Mol Phylogenet Evol. 2010;57: 798–811. pmid:20723609

27. Takami Y, Osawa T. Ecological differentiation and habitat unsuitability maintaining a ground beetle hybrid zone. Ecol Evol. 2016;6: 113–124. pmid:26811778

28. Sullivan AR, Owusu SA, Weber JA, Hipp AL, Gailing O. Hybridization and divergence in multi-species oak (Quercus) communities. Bot J Linn Soc. 2016;181: 99–114.

29. Cabaña I, Chiaraviglio M, Di Cola V, A G, Broennimann O, Gardenal CN, et al. Hybridization and hybrid zone stability between two lizards explained by population genetics and niche quantification. Zool J Linn Soc. 2020;190: 757–769.

30. Sánchez-Guillén RA, Muñoz J, Hafernik J, Tierney M, Rodriguez-Tapia G, Córdoba-Aguilar A. Hybridization rate and climate change: Are endangered species at risk? J Insect Conserv. 2014;18: 295–305.

31. Morales-Rozo A, Tenorio EA, Carling MD, Cadena CD. Origin and cross-century dynamics of an avian hybrid zone. BMC Evol Biol. 2017;17: 1–18.

32. Miranda F, Hernández-X E. Los tipos de vegetación de México y su clasificación. Boletín de la Sociedad Botánica de México; 1963. https://dx.doi.org/10.17129/botsci.1084

33. Rzedowski J. La vegetación de México. Second. Limusa; 1978.

34. Hernández-Oria J. Desaparición del Bosque Seco en El Bajío mexicano: Implicaciones del ensamblaje de especies y grupos funcionales en la dinámica de una vegetación amenazada. Zo Áridas. 2007;11: 13–31.

35. Guimarães R, Forni-Martins ER. Chromosome numbers and their evolutionary meaning in the Sapindales order: an overview. Brazilian J Bot 2021. 2021; 1–15.

36. Doyle JJ, Doyle JL. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem Bull. 1987;19: 11–15. Available: http://www.sciepub.com/reference/60627

37. Edgar RC. MUSCLE: Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32: 1792–1797. pmid:15034147

38. Kumar S, Stecher G, Li M, Knyaz C, Tamura K. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol Biol Evol. 2018;35: 1547–1549. pmid:29722887

39. S M, S NJ, D P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68: 978–989. pmid:11254454

40. Rozas J, Ferrer-Mata A, Sanchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, et al. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34: 3299–3302. pmid:29029172

41. Huson DH, Bryant D. Application of Phylogenetic Networks in Evolutionary Studies. Mol Biol Evol. 2006;23: 254–267. pmid:16221896

42. M D, R E. RDP: detection of recombination amongst aligned sequences. Bioinformatics. 2000;16: 562–563. pmid:10980155

43. Darriba D, Taboada GL, Doallo R, Posada D. JModelTest 2: More models, new heuristics and parallel computing. Nature Methods. Nature Publishing Group; 2012. p. 772. https://doi.org/10.1038/nmeth.2109 pmid:22847109

44. Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2019;15: e1006650. pmid:30958812

45. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29: 1969–1973. pmid:22367748

46. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67: 901–904. pmid:29718447

47. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7: 1–8.

48. Leigh JW, Bryant D. POPART: Full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6: 1110–1116.

49. Excoffier L, Lischer HEL. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10: 564–567. pmid:21565059

50. Cheng L, Connor TR, Sirén J, Aanensen DM, Corander J. Hierarchical and spatially explicit clustering of DNA sequences with BAPS software. Mol Biol Evol. 2013;30: 1224–1228. pmid:23408797

51. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet 2010;11: 1–15.

52. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005.

53. Hernández-Pérez E, González-Espinosa M, Trejo I, Bonfil C. Distribution of the genus Bursera in Morelos state (Mexico) and its relation to climate. Rev Mex Biodivers. 2011;82: 964–976.

54. Gámez N, Escalante T, Espinosa D, Eguiarte LE, Morrone JJ. Temporal dynamics of areas of endemism under climate change: A case study of Mexican Bursera (Burseraceae). J Biogeogr. 2014;41: 871–881.

55. Moisés-Méndez T. Factores ambientales gradiente bosquel toropical caducifolio. 2015.

56. Barve N, Barve V, Jiménez-Valverde A, Lira-Noriega A, Maher SP, Peterson AT, et al. The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecol Modell. 2011;222: 1810–1819.

57. Hijmans RJ. Raster Geographic Data Analysis and Modeling. 2015. https://cran.r-project.org/web/packages/raster/index.html

58. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Modell. 2006;190: 231–259.

59. Cobos ME, Townsend Peterson A, Barve N, Osorio-Olvera L. Kuenm: An R package for detailed development of ecological niche models using Maxent. PeerJ. 2019;2019: e6281. pmid:30755826

60. Peterson AT, Papeş M, Soberón J. Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol Modell. 2008;213: 63–72.

61. QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation. 2009. http://qgis.org

62. Funk DJ, Nosil P, Etges WJ. Ecological divergence exhibits consistently positive associations with reproductive isolation across disparate taxa. Proc Natl Acad Sci U S A. 2006;103: 3209–3213. pmid:16492742

63. Broennimann O, Fitzpatrick MC, Pearman PB, Petitpierre B, Pellissier L, Yoccoz NG, et al. Measuring ecological niche overlap from occurrence and spatial environmental data. Glob Ecol Biogeogr. 2012;21: 481–497.

64. Schoener TW. The Anolis Lizards of Bimini: Resource Partitioning in a Complex Fauna. Ecology. 1968.

65. Warren DL, Glor RE, Turelli M. Environmental niche equivalency versus conservatism: Quantitative approaches to niche evolution. Evolution (N Y). 2008;62: 2868–2883. pmid:18752605

66. Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37: 4302–4315.

67. Otto-Bliesner BL, Marshall SJ, Overpeck JT, Miller GH, Hu A. Simulating arctic climate warmth and icefield retreat in the last interglaciation. Science (80-). 2006;311: 1751–1753. pmid:16556838

68. Watanabe S, Hajima T, Sudo K, Nagashima T, Takemura T, Okajima H, et al. MIROC-ESM 2010: Model description and basic results of CMIP5-20c3m experiments. Geosci Model Dev. 2011;4: 845–872.

69. Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, et al. The community climate system model version 4. J Clim. 2011;24: 4973–4991.

70. Stocker TF, Qin D, Plattner G-K, Tignor MMB, Allen SK, Boschung J, et al. Climate Change 2013 The Physical Science Basis Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change Edited by. 2013. www.cambridge.org

71. Tocker TF, Qin D, Plattner G-K, Tignor MMB, Allen SK, Boschung J, et al. Climate Change 2013 The Physical Science Basis. Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. 2013.

72. van Vuuren DP, Edmonds J, Kainuma M, Riahi K, Thomson A, Hibbard K, et al. The representative concentration pathways: An overview. Clim Change. 2011;109: 5–31.

73. Medina-Lemos R. Burseraceae. Fascículo 66. Flora del Valle de Tehuacán-Cuicatlán; 2008. http://www.mobot.org/

74. Weeks A, Simpson BB. Molecular phylogenetic analysis of Commiphora (Burseraceae) yields insight on the evolution and historical biogeography of an “impossible” genus. Mol Phylogenet Evol. 2007;42: 62–79. pmid:16904915

75. Becerra JX. Evolution of Mexican Bursera (Burseraceae) inferred from ITS, ETS, and 5S nuclear ribosomal DNA sequences. Mol Phylogenet Evol. 2003;26: 300–309. pmid:12565038

76. Jerry A. Coyne, H. Allen Orr. Speciation -. In: Sinauer Associates [Internet]. 2004 [cited 14 Jun 2021] pp. 1–545. https://global.oup.com/academic/product/speciation-9780878930890?cc=mx&lang=en&

77. Rico Y, Reyes-Estanislao L. Pollen viability and germinability of putative Bursera hybrids (section Bullockia; Burseraceae) in Mexico. Acta Bot Mex. 2019;126: 1–10.

78. Jakob SS, Blattner FR. A Chloroplast Genealogy of Hordeum (Poaceae): Long-Term Persisting Haplotypes, Incomplete Lineage Sorting, Regional Extinction, and the Consequences for Phylogenetic Inference. Mol Biol Evol. 2006;23: 1602–1612. pmid:16754643

79. Arnold ML. Natural hybridization as an evolutionary process. Annu Rev Ecol Syst. 1992;23: 237–261.

80. Goad DM, Baxter I, Kellogg EA, Olsen KM. Hybridization, polyploidy and clonality influence geographic patterns of diversity and salt tolerance in the model halophyte seashore paspalum (Paspalum vaginatum). Mol Ecol. 2021;30: 148–161. pmid:33128807

81. Reséndiz-López MA. Analisis de la distribucion ecológica de Bursera Jacq. ex. L. (Burseraceae), Sección Bullockia en México. Univesidad Nacional Autónoma de México. 2009.

82. Klein EK, Lagache-Navarro L, Petit RJ. Demographic and spatial determinants of hybridization rate. J Ecol. 2017;105: 29–38.

83. Mandeville EG, Parchman TL, Thompson KG, Compton RI, Gelwicks KR, Song SJ, et al. Inconsistent reproductive isolation revealed by interactions between Catostomus fish species. Evol Lett. 2017;1: 255–268. pmid:30283654

84. Gross BL, Rieseberg LH. The ecological genetics of homoploid hybrid speciation. Journal of Heredity. Oxford Academic; 2005. pp. 241–252. pmid:15618301

85. Agrawal AF, Feder JL, Nosil P. Ecological divergence and the origins of intrinsic postmating isolation with gene flow. Int J Ecol. 2011.

86. Tarroso P, Pereira RJ, Martínez-Freiría F, Godinho R, Brito JC. Hybridization at an ecotone: Ecological and genetic barriers between three Iberian vipers. Mol Ecol. 2014;23: 1108–1123. pmid:24447270

87. Hekkala ER, Platt SG, Thorbjarnarson JB, Rainwater TR, Tessler M, Cunningham SW, et al. Integrating molecular, phenotypic and environmental data to elucidate patterns of crocodile hybridization in Belize. R Soc Open Sci. 2015;2: 150409. pmid:26473062

88. Otis JA, Thornton D, Rutledge L, Murray DL. Ecological niche differentiation across a wolf-coyote hybrid zone in eastern North America. Divers Distrib. 2017;23: 529–539.

89. Manzoor SA, Griffiths G, Obiakara MC, Esparza-Estrada CE, Lukac M. Evidence of ecological niche shift in Rhododendron ponticum (L.) in Britain: Hybridization as a possible cause of rapid niche expansion. Ecol Evol. 2020; 1–11. pmid:32128136

90. Barton NH. The role of hybridization in evolution. Mol Ecol. 2001;10: 551–568. pmid:11298968

91. Rieseberg LH, Raymond O, Rosenthal DM, Lai Z, Livingstone K, Nakazato T, et al. Major ecological transitions in wild sunflowers facilitated by hybridization. Science (80-). 2003;301: 1211–1216. pmid:12907807

92. Aldridge G. Variation in frequency of hybrids and spatial structure among Ipomopsis (Polemoniaceae) contact sites. New Phytol. 2005;167: 279–288. pmid:15948849

93. Carson EW, Tobler M, Minckley WL, Ainsworth RJ, Dowling TE. Relationships between spatio-temporal environmental and genetic variation reveal an important influence of exogenous selection in a pupfish hybrid zone. Mol Ecol. 2012;21: 1209–1222. pmid:22269008

94. Maxwell LM, Walsh J, Olsen BJ, Kovach AI. Patterns of introgression vary within an avian hybrid zone. BMC Ecol Evol. 2021;21: 14. pmid:33509089

95. Aldridge G, Campbell DR. Genetic and morphological patterns show variation in frequency of hybrids between Ipomopsis (Polemoniaceae) zones of sympatry. Heredity (Edinb). 2009;102: 257–265. pmid:18971956

96. Campbell LG, Shukla K, Sneck ME, Chaplin C, Mercer KL. The effect of altered soil moisture on hybridization rate in a crop-wild system (Raphanus spp.). PLoS One. 2016;11: e0166802. pmid:27936159

97. Beddows I, Rose LE. Factors determining hybridization rate in plants: A case study in Michigan. Perspect Plant Ecol Evol Syst. 2018;34: 51–60.

98. Mastretta-Yanes A, Moreno-Letelier A, Piñero D, Jorgensen TH, Emerson BC. Biodiversity in the Mexican highlands and the interaction of geology, geography and climate within the Trans-Mexican Volcanic Belt. J Biogeogr. 2015;42: 1586–1600.

99. Prieto-Torres DA, Navarro-Sigüenza AG, Santiago-Alarcon D, Rojas-Soto OR. Response of the endangered tropical dry forests to climate change and the role of Mexican Protected Areas for their conservation. Glob Chang Biol. 2016;22: 364–379. pmid:26367278

100. Christmas MJ, Breed MF, Lowe AJ. Constraints to and conservation implications for climate change adaptation in plants. Conserv Genet. 2016;17: 305–320.