Highly divergent sympatric lineages of Leptotila verreauxi (Aves: Columbidae) suggest a secondary contact area in the Isthmus of Tehuantepec, Mexico

1.   Introduction

Mesoamerica is characterized by a complex geological history, a rugged topography (see Ferrusquia-Villafranca, 1993), and complex paleoclimatic dynamics (Shackleton, 2000; Wogau et al., 2019). These characteristics have been highlighted as major drivers that have shaped the biotic evolutionary history in the region (Morrone, 2010; Ornelas et al., 2013; Stull, 2023). Some of the most important effects of geological and climatic dynamics on Mesoamerican biota are in range shifts in the species distribution, demographic history changes, and genetic divergence, all of which has motivated the increase of phylogeographic studies throughout the region (García-Moreno et al., 2004; Navarro-Sigüenza et al., 2008; Vázquez-Miranda et al., 2009; Ortiz-Ramírez et al., 2016; Rodríguez-Correa et al., 2017).

Located in southeastern Mexico, the Isthmus of Tehuantepec (IT) has been long recognized as one of the most important geographic features in the diversification of a large number of taxa in northern Mesoamerica (Binford, 1989; Peterson et al., 1999; Barber and Klicka, 2010; Ornelas et al., 2013). Studies in montane-forest birds distributed east and west of the IT have suggested its crucial role as a geographic barrier acting in different times or “pulses” (Barber and Klicka, 2010; Ornelas et al., 2013; Tsai et al., 2019). However, the role that IT has played in the diversification of lowland birds is still poorly understood (see Sosa-López et al., 2013), probably because it has been long assumed that lowland bird species can disperse freely throughout the region (Edwards and Lea, 1955; Binford, 1989), thus preventing complete divergence in these taxa (e.g., Howell and Webb, 1995).

The IT represents the narrowest zone in Mesoamerica, allowing contact between the dry tropical lowland forests of the Pacific slope, and the humid tropical lowland forests of the Gulf slope. This gradient has been defined as an ecological barrier that may have limited gene flow between populations in a north to south axis, promoting phenotypic differentiation of lowland populations, and even leading to the recognition of subspecies in each habitat (Binford, 1989). In addition, a few lowland species show an east-west differentiation pattern (e.g., Sánchez-González et al., 2023). Hence, the geographical and ecological history of the IT implies intricate biogeographical patterns and processes within the region (Friedmann et al., 1950; Binford, 1989; Hogan, 1999; Zarza et al., 2008; Ornelas et al., 2013; Arcangeli et al., 2018; Gray et al., 2019; Castillo-Chora et al., 2021; Sánchez-González et al., 2021; Llanes-Quevedo et al., 2022).

The White-tipped Dove (Leptotila verreauxi) is a widely distributed polytypic species mainly distributed in tropical lowlands from extreme southern USA (the Lower Rio Grande Valley in southern Texas) through Mexico and Central America to northern Argentina (Friedmann et al., 1950). Currently, 13 subspecies are recognized (Gill et al., 2021), four of which are distributed in Mexico (Fig. 1): capitalis, endemic to the Tres Marias Islands; angelica, from southern Texas in the Gulf slope, and southern Sonora in the Pacific slope south to the IT; fulviventris, from the IT on the Gulf slope to the Yucatan Peninsula; and bangsi, from the southern IT and Chiapas (Friedmann et al., 1950) to Honduras and western Nicaragua (Hogan, 1999; Gill et al., 2021). The ranges of the three Mexican mainland L. verreauxi subspecies abut at the IT (Friedmann et al., 1950), although precise boundaries between these are unclear, suggesting a complex biogeographic history in the region. As currently described, the main differentiation occurs in a west-east axis (angelica – bangsi, fulviventris). Such west to east differentiation, apparently follows a geographic pattern opposite to what may be expected given the environmental gradient. A previous study on the biogeography of New World doves supported a similar pattern of genetic differentiation, however details on geographical/ecological limits for this pattern are unclear (Johnson and Weckstein, 2011).

Figure  1.  Geographical representation of the samples used in this study. Color of the circles correspond to detected clades. Green, NA clade; Blue, NCA clade; Yellow, YUC clade. Red and pink: South American clades. (A) Inset 1 shows the distribution of the haplogroups and the number of individuals in each locality in the Isthmus of Tehuantepec (IT). The black star (subdivided in the two-colored pie charts) at the southern IT depicts the secondary contact zone between the NA and NCA clades. Colors depict vegetation types at the IT. Dashed line shows putative boundaries between currently recognized subspecies: A = L. v. angelica, B = L. v. bangsi, C = L. v. fulviventris. (B) Inset 2 depicts the Tres Marías Islands, where the endemic L. v. capitalis is found. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

Here, we described genetic differentiation patterns at the IT using L. verreauxi as a model. Based on the distribution of three subspecies, which abut at the region (Fig. 1), as a starting framework, we expected: a) that genetic differentiation corresponds to the ecological structure, in which main differentiation would occur in a north to south axis across the environmental gradient. A west-east pattern of differentiation would not be expected due to the ecological continuity of the lowland forests in the region which has promoted the recognition of mainly phenotypically-based bird subspecies across a north to south ecological gradient (e.g., Howell and Webb, 1995); b) ecological niche differentiation, which would support the recognized differentiation in current taxonomic treatments; and c) that contact zones between differentiated taxa will be maintained by local environmental differences across the ecological gradient. We performed phylogenetic and population genetics analyses using two mitochondrial and one nuclear markers, and estimated divergence times between lineages to explore genetic variation and historical population dynamics, focusing particularly on the populations at the IT. In addition, we used correlative ecological niche models (ENM) to estimate paleodistributions of lineages and test hypotheses of niche conservatism, as well as evaluate the role of ecological differences as a biogeographic limit in these populations. All of these tests allowed us to establish if these distribution patterns may be the result of secondary contact, which would suggest a vicariant event at the IT. We also evaluated if there are differences in the ecological niche that could be contributing to the maintenance of differentiation.

2.   Materials and methods

2.1   Taxon sampling and laboratory protocols

We obtained 48 tissue samples from the collection of the Museum of Zoology “Alfonso L. Herrera” (MZFC), representing all Mexican subspecies (Appendix Table S1). We supplemented our samples with 20 sequences available in GenBank (Benson et al., 2013; http://www.ncbi.nlm.nih.gov; Appendix Table S1). Outgroups were selected based on a previous phylogenetic hypothesis (Johnson and Weckstein, 2011) and included six species: L. jamaicensis, L. cassini, L. megalura, L. rufaxilla, L. plumbeiceps and Leptotrygon veraguensis. In total, we included 68 individual genetic sequences, considering samples deposited in GenBank.

DNA was extracted using both a Qiagen DNeasy Blood and Tissue Kit following the manufacturer’s protocol, as well as a phenol-chloroform extraction protocol (Sambrook and Russell, 2006). We amplified and sequenced two mitochondrial DNA markers [NADH, dehydrogenase subunit 2 (ND2) and the cytochrome C oxidase subunit 1 (COI)]; and one nuclear marker [the β-fibrinogen intron 7 (β-fibint 7)] via the polymerase chain reaction (PCR) following standardized protocols (Lorenz, 2012). We used primers L5215-H5766 and 487L-H6313 (Sorenson et al., 1999) for ND2, L6626-H7005 (Hafner et al., 1994) for COI, and B17L-B17U (Prychitko and Moore, 1997) for β-fibint 7. These molecular markers provide good resolution when analyzing phylogenetic relationships and biogeographic patterns in Columbiformes, in addition to show high phylogenetic consistency between mtDNA and nDNA in some studies at the supraspecific level (Johnson and Clayton, 2000; Pereira et al., 2007; Johnson and Weckstein, 2011; Sweet and Johnson, 2015; Sweet et al., 2017; Sangster et al., 2018). PCR products were visualized in agarose gels (1.3%) stained with ethidium bromide (0.5 μL/mL). Purification and sequencing processes were performed at Macrogen Inc., South Korea. Chromatograms were edited using Geneious Prime 2019.0.4 (Kearse et al., 2012) and multiple alignments were performed using MAFFT (Katoh and Standley, 2013) with default parameters.

2.2   Model selection and phylogenetic analyses

Partition schemes and sequence evolution models for Bayesian Inference were selected using PartitionFinder v2.1.1 (Lanfear et al., 2017). Molecular evolution models selected for each of the partition schemes are shown in Appendix Table S2. Phylogenies were generated using Maximum Likelihood (ML) and Bayesian Inference (BI). For ML, we used RAxMLGUI (Silvestro and Michalak, 2012), with the selected partition scheme, the GTR substitution model and the algorithm “ML + rapid bootstrap” (Stamatakis et al., 2008), and 10,000 iterations. We used MrBayes for estimating BI trees, which were obtained using the partition schemes and selected models as previously described (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003). We ran 10,000,000 generations in four parallel runs and eight Markov Chain Monte Carlo (MCMC) sampling every 1000 generations. We used Tracer 1.7 (Rambaut et al., 2018) to determine a stationary phase and the burn-in, and used FigTree v1.4.3 (Rambaut, 2009) to visualize the consensus tree.

2.3   Divergence time estimation

Divergence times were estimated in BEAST v.1.10.4 (Suchard et al., 2018) using the partition schemes and substitution models previously selected. We applied a strict molecular clock with a Yule speciation process (Steel and McKenzie, 2001). Since we were analyzing genetic patterns at intraspecific level, we opted for a strict molecular clock, which assumes a single rate of constant evolution throughout the phylogeny (Pybus, 2006). We applied a divergence rate of 1.96% (0.0098 ± 0.0005 substitutions/site/Myr) for the mitochondrial markers (Sweet and Johnson, 2015), and a divergence rate of 0.53% (0.00265 ± 0.0005 substitutions/site/Myr) for β-fibint 7 (Johnson and Clayton, 2000). Two independent runs of 20,000,000 generations were generated starting from different random points, corroborating the −lnL values and values of the effective sample sizes (ESSs) > 200 for all of the parameters in Tracer 1.7 (Rambaut et al., 2018). We used LogCombiner for combining the output files and TreeAnnotator (Drummond and Rambaut, 2015) for obtaining a maximum credibility tree (MCCT) with a burn-in of 25% and a posterior probability limit of 0.5%.

2.4   Phylogeography

As a measure of genetic differentiation between populations, we calculated genetic fixation from the F_ST estimator, and interpreted its values based on Hartl and Clark (1997). From this index, the rate of gene flow between populations or N_m can be estimated, the higher this value the higher the gene flow (Nei, 1987). Genetic distances were measured using the Nei corrected distance parameter D_XY (Nei, 1987). Populations were identified according to the phylogeny and haplotype networks. Historical demography analyses for each clade were calculated using the R₂ test, which has been shown to have a superior statistical power in small sample sizes (Ramos-Onsins and Rozas, 2002). The significance of the R₂ test was estimated through 1000 coalescent simulations for each clade. All phylogeographic tests and genetic diversity estimates were obtained in DnaSP v.6.12.01 (Rozas et al., 2017). Given that ND2 and COI have different sample size numbers, we built different haplotype networks in Network 10.1 using a Median Joining Network algorithm (Bandelt et al., 1999). We assessed genetic variation in L. verreauxi via an Analysis of Molecular Variance (AMOVA) in Arlequin 3.5.2.2 (Excoffier et al., 2005). Significance was estimated via two independent analyses with 16,000 permutations based on pairwise differences: one included all populations, while another included only populations present in the IT, to test whether the levels of variation were consistent when only these populations are considered.

2.5   Ecological niche models

Since Pleistocene climatic changes have significantly influenced the evolution and biogeography of lowland taxa in the IT (Castillo-Chora et al., 2021; Sánchez-González et al., 2023), we estimated paleodistribution models for populations of L. verreauxi during the Middle Holocene (MH, ~6000 years ago), Last Glacial Maximum (LGM, ~21,000 years ago), and Last Interglacial (LIG, ~120,000–140,000 years ago), as well as for the Present. For model building, we used the data corresponding to occurrence records of L. verreauxi that we used in molecular analyses, supplemented with those downloaded from the Global Biodiversity Information Facility (GBIF; https://doi.org/10.15468/dl.f7wwag). Duplicated and erroneously georeferenced records were filtered and eliminated for analyses, and we performed a randomized spatial thinning with the R package “spThin” (Aiello-Lammens et al., 2015), which allowed us to build highly accurate models by reducing biases associated to data concentration in some geographic areas. Based on the genetic structure, we divided the total data into three datasets corresponding to each lineage present in the IT. Given that geographic boundaries for each lineage in this lowland taxon are not clear, we generated polygons from occurrence data corresponding to unequivocally assigned genetic samples; this procedure allowed us to reduce identification bias for records from biological collections identified below the species level (Appendix Fig. S1).

We downloaded environmental variables for the present and for past scenarios corresponding to the MH, LGM, and LIG from WorldClim 1.4 (Hijmans et al., 2005) with a resolution of 2.5 min. For the MH and the LGM models, we used layers from three different Atmosphere-Ocean General Circulation Models: CCSM4 (Community Climate System Model) (Gent et al., 2011), MIROC-ESM (Model for Interdisciplinary Research on Climate) (Hasumi and Emori, 2004), and MPI-ESM-P (Max Planck Institute) (Baehr et al., 2015); layers corresponding to the present and the LIG were downloaded from the CCSM4 model. We retained environmental variables with a greater contribution to the present distribution ranges, while variables with high collinearity were eliminated via the variance inflation factor (VIF), with a VIF >10 as implemented in the R package “corrplot” (Wei and Simko, 2021). We therefore selected nine bioclimatic variables: mean diurnal range (BIO2), isothermality (BIO3), max temperature of warmest month (BIO5), mean temperature of wettest quarter (BIO8), precipitation of wettest month (BIO13), precipitation of driest month (BIO14), precipitation seasonality (BIO15), precipitation of warmest quarter (BIO18) and precipitation of coldest quarter (BIO19).

We used the “Grinnell” v.0.0.21 R package (Machado-Stredel et al., 2021) to define an accessibility area (M) for each lineage (Appendix Fig. S2) based on estimates of the historical accessible area derived from simulations considering biological processes such as dispersal, colonization, and information on climate change over time. These M areas were further refined through comparisons with the terrestrial ecoregions (Dinerstein et al., 2017), and biogeographic provinces of the Neotropics layers (Morrone, 2014). All environmental models corresponding to the present and projections to the past were generated using the maximum entropy algorithm in Maxent v3.4.4 (Phillips et al., 2017) as implemented in “KUENM” v 1.1.9 package (Cobos et al., 2019), that is designed to strengthen the model calibration and creation process, via the generation of a set of candidate models to select one or several optimal models from certain parameterizations. Finally, to explore potential regions with low reliability for model transference to the past, we performed extrapolation risk analyses via the Mobility-oriented parity (MOP) test in “NicheToolbox” (Osorio-Olvera et al., 2020). The final models were edited in QGIS v3.18 (QGIS Development Team, 2022; http://www.qgis.org/)

2.6   Niche conservatism and range-breaking tests

Due to the absence of conspicuous geographic barriers across the lowland distributional range of L. verreauxi, we explored if genetic structure has been maintained through ecological niche differentiation in the different lineages detected in Mesoamerica. We compared the ecological niches through paired tests of niche conservatism between each lineage; additionally, we tested whether environmental differences across the IT may act as an ecological barrier separating these lineages. All these analyses were performed using “ENMTools” (Warren et al., 2021).

Niche conservatism models were generated with Maxent models (Phillips et al., 2017) to perform pairwise comparisons between the observed lineages. These analyses compare the overlapping of empirical models with a null distribution constructed from a set of pseudo-replicates using Monte Carlo methods. Empirical and pseudoreplicated values were calculated from measures of niche overlap (Schoener’s D) via 100 replicates to generate the null distribution across all tests. We test whether niches differ from what is expected by chance via a niche equivalence test (or identity test) and the symmetric background similarity test (Warren et al., 2021). Niche equivalence and similarity tests were conducted using model-based techniques test as implemented in ENMTools (Warren et al., 2021). To assess the role of environmental differences as a barrier for dispersal of the lineages across the IT, we used blob and ribbon range-breaking tests (Glor and Warren, 2011). Blob range-breaking test evaluates whether regions occupied by each lineage are more environmentally different than expected by chance; if so, the environmental differences would represent an abrupt biogeographical barrier. On the other hand, ribbon range-breaking test assesses whether unsuitable habitat divides areas of suitable habitat for each of the lineages which may function as a biogeographical boundary (Warren et al., 2021). We selected two areas that may functioning as a barrier for dispersal (the “ribbon”) in this lowland taxon: 1) the Sierra del Tolistoque, in the southern IT, and 2) the wetland region in the southern Gulf of Mexico. These two areas may serve as barriers due to apparent environmental differences on each side, maintaining lineage divergence despite confirmed overlap between the NA and NCA lineages in the southern IT. To optimize the ribbon range-breaking test, we generated a set of random points within the two selected areas due to the small number of occurrence records.

3.   Results

We obtained 48 sequences corresponding to ND2 (1041 base pairs, bp), 26 to COI (379 bp) and 16 to β-fibint 7 sequences (1136 bp). mtDNA topologies obtained from BI and ML analyses were highly similar. The tree showed two highly supported monophyletic groups for L. verreauxi (Fig. 1), separated by the IT. One group included the North American (NA) clade, distributed west of the IT, and a few from northwestern Chiapas (geographically east of the IT); the other main group, distributed east of the IT, is subdivided in two groups, each of which present further well-supported groupings. Samples from the Soconusco region in the Pacific slope in southeastern Mexico were grouped in a Northern Central American Clade (NCA), which is sister to the Yucatán Peninsula clade (YUC), while samples from Southern Central America and northern South America (SCA) formed a clade sister to west and central South America samples (SA).

Our mtDNA phylogeographic analyses showed two sympatric lineages in the southern IT region (Fig. 1). The NA lineage is distributed across the region in both the Pacific and Gulf slopes, while the NCA lineage is apparently restricted to the dry forests south of the IT in the Pacific slope. The nDNA tree (Appendix Fig. S3) did not show a clear geographic structure. Most of the samples of Leptotila species used as outgroup were included in a group sister to several L. verreauxi samples. The other group included most of samples distributed west of the IT (angelica and capitalis), two samples from NCA (bangsi), and an outgroup species, L. jamaicensis. These groups shared indels that coincide in position and length, with low nodal support (PP > 0.76 and BS > 70).

3.1   Divergence times

Our data suggest a divergence event during the Late Pliocene – Lower Pleistocene transition at approximately 2.5 Mya [1.94–3.43 Mya, 95% Height Posterior Density (HPD)] between the group west of the IT (the NA clade) and the group east of the IT. Further divergence events in the group east of the IT occurred between the clades SCA–SA and the clades NCA and YUC during the Lower Pleistocene around 1.89 Mya (1.39–2.49, HPD). The most recent divergences occurred between NCA and Yucatan at 0.9 Mya (0.62–1.34, HPD) and between SCA and SA at 0.92–1.78 Mya. A recent differentiation event between populations of the Tres Marías islands (L. v. capitalis) and rest of NA clade occurred around 100 kya (0.04–0.18, HDP).

3.2   Phylogeography

The haplotype network for the mDNA markers showed the same structure as the phylogeographic tree (Fig. 2). We found 25 haplotypes in ND2 and 8 haplotypes in COI. In both haplotype networks, we found a relatively large number of mutational steps separating two main clades roughly divided at the IT (25 in ND2 and 11 in COI), while NCA and Yucatán clades are separated by a smaller number of mutational steps (14 in ND2 and 5 in COI). We detected a star-like pattern in the NA clade; shapes are unclear for the other clades. The demographic R₂ test showed population expansion for the NA clade (P < 0.05); however, R₂ values for NCA and YUC were not statistically significant (P > 0.05 for both).

Figure  2.  Phylogenetics of Leptotila verreauxi. (A) mtDNA haplotype network. Black dots represent mutational steps, while circles indicate possible extinct or unsampled haplotypes. Circle size represent the frequency of haplotypes. (B) mtDNA tree. Numbers above nodes correspond to Bootstrap/posterior probability values; numbers below nodes indicate 95% HPD (High Posterior Density) ranges for divergence dates (in Mya). The dotted black box groups samples from the Tres Marías Islands. Red labels in the terminals correspond to sympatric samples from the Sierra del Tolistoque in the southern IT. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

The F_ST, N_m, and D_XY indices showed a strong differentiation in most of the populations (Table 1). All the F_ST values are close to 1.0, signifying high genetic fixation and nearly complete haplotypic separation, in all lineages, as expected in genetically structured populations. N_m in the other hand, indicated less than one migrant every four generations. Genetic distances, as measured by the Nei corrected distance parameter D_XY revealed high levels of genetic differentiation (>2%) for all paired comparisons, except for NCA-YUC with 0.00748 (0.7%).

Table  1.  F_ST values (below the diagonal) and D_XY values (above the diagonal) for the different L. verreauxi lineages. Values in the parentheses represent the number of migrants (N_m) per generation.

	NA	NCA	YUC
NA	–	0.03840	0.03659
NCA	0.97807 (0.01)	–	0.00748
YUC	0.97006 (0.01)	0.95563 (0.01)	–

 | Show Table

DownLoad: CSV

The AMOVA analysis (Table 2) including all populations showed that 95.9% of the genetic variation is found between populations, while only 4.06% is found within populations. The results are similar when the NA-NCA lineages are included.

Table  2.  Results of AMOVA of the L. verreauxi populations.

Populations	Source of variation	d.f.	Sum of squares	Variance components	Percentage variation	FST	P
All populations	Among populations	4	577.492	18.09475 Va	95.94	0.95936	0.0001
	Within populations	50	38.326	0.76652 Vb	4.06
NA–NCA	Among populations	1	230.745	17.76667 Va	96.60	0.96603	0.0001
	Within populations	40	24.993	0.62482 Vb	3.40

 | Show Table

DownLoad: CSV

3.3   Distributions and paleodistributions

Distribution models for the present and projections to the past showed high performance in each of the lineages (Table 3), with models that are significantly better than expected by chance (AUC ratio >1) and p-values <0.001. For the NA clade, environmental variables with highest contribution were precipitation seasonality and precipitation of wettest month (BIO 15 = 28.1%, BIO 13 = 22.3%, respectively); while for NCA were isothermality and mean diurnal range (BIO 3 = 48.3%, BIO 2 = 30.4%, respectively). In the case of YUC, precipitation of wettest month and precipitation of warmest quarter (BIO 13 = 48.4%, BIO 18 = 8.9%, respectively).

Table  3.  Summary of model performance statistics per lineage.

Lineage	Regularization multiplier	Feature classes	AUC ratio	AICc	Omission rate	Parameters
NA	0.25	qp	1.4773	21022.95	0.05	35
NCA	3	q	1.5733	5963.72	0.03	7
YUC	0.5	lq	1.2037	4721.84	0.04	15

 | Show Table

DownLoad: CSV

The areas of high environmental suitability for the three Mesoamerican lineages showed clear geographic changes in comparison to present (Fig. 3). During the LIG, NCA and YUC showed a smaller geographic distribution in small and isolated habitat patches. NA showed an increase in the geographic extension of the environmental suitability in southeastern Mexico and northern Central America. During the LGM, NA experienced a reduction in suitable areas, being restricted to three main lowland areas in both slopes and an almost complete disappearance of suitability areas south of the IT. Regions such as YUC exhibited a decrease in environmental suitability, while NCA lineage showed an increase of the suitable conditions towards the MH. It is important to observe that all the three lineages experienced an increase in the environmental suitable areas; these areas have been apparently maintained towards the Present. According to environmental suitability maps, current contact zones across IT may have not been formed until recently (at least since LGM), as most of suitable areas appear as allopatric through time intervals for which ecological modeling was conducted.

Figure  3.  Ecological niche models of each lineage during the Present, Middle Holocene, Last Glacial Maximum and Last Interglacial. Dark areas in each map in the Present correspond to areas of climatic stability.

DownLoad: Full-Size Img PowerPoint

3.4   Ecological niche comparisons

Because results of D were similar and I tend to express greater variability, for niche comparisons we focused on the Schoener’s D. The niche identity tests revealed that the niches of each of the three lineages are not identical (i.e, equivalent; empirical Schoener D values lower than the pseudoreplicated values; P < 0.01; Table 4; Appendix Fig. S4). On the other hand, the background similarity tests did not reject the null hypothesis for the NA-NCA comparison, while in the other two comparisons, the null hypothesis was rejected, indicating niche divergence (i.e., the ecological niches of these clades are less similar than expected based on their background environments).

Table  4.  Schoener’s D values of the niche comparison tests (Identity, Similarity and Blob rangebreak).

Lineages	Identity test	Similarity test	Blob rangebreak
NA–NCA	0.297**	0.271ns	0.282ns
NA–YUC	0.169**	0.148*	0.159ns
YUC–NCA	0.090**	0.103*	0.099ns
The significance of the test is shown ns, non-significant; P < 0.05; **P < 0.01.

 | Show Table

DownLoad: CSV

Range-breaking blob tests did not show empirical values with a significant deviation from the null hypothesis (P > 0.05; Table 4; Appendix Fig. S4), suggesting that no biogeographical barrier is associated with environmental differences across the IT ecological gradient. However, the ribbon range-breaking test for the NA and NCA lineages suggested the presence of a region with less suitable conditions for the NCA lineage in the contact zone (D value, P < 0.01 for the comparison of the NCA lineage vs. the ribbon; Table 5; Appendix Fig. S5). Similarly, the comparison between the NA and YUC lineages suggests the presence of a region of less suitable conditions for the YUC lineage (D value, P < 0.05 for YUC vs. the ribbon; Table 5; Appendix Fig. S6).

Table  5.  Schoener’s D values from the ribbon-rangebreak tests.

L1	L2	L1 vs. L2	L1 vs. Ribbon	L2 vs. Ribbon	Outside vs Ribbon
NA	NCA	0.297ns	0.184ns	0.049**	0.154**
NA	YUC	0.086**	0.253ns	0.263*	0.326ns
The significance of the test is shown ns, non-significant; *, P < 0.05; **, P < 0.01.

 | Show Table

DownLoad: CSV

4.   Discussion

Our study illustrates both emergent patterns found in different taxa across the study region, as well as discrepancies between phenotypically-based and molecular-based taxonomies in L. verreauxi.

4.1   Phylogeography

Our phylogenetic tree (Fig. 2) agrees with divergence events described in other studies along the Mexican Pacific slope and the IT in which previously unexpected evolutionary patterns are being revealed (Arcangeli et al., 2018; Castillo-Chora et al., 2021; Sánchez-González et al., 2021, 2023; Llanes-Quevedo et al., 2022). On the other hand, our analyses also show that traditional phenotypically-based taxonomy (Fig. 1) may either fail to correspond with evolutionary patterns as differentiation may lie further east from the IT, perhaps up to the Usumacinta River delta, as in other taxa (Howell and Webb, 1995; Bradley et al., 2008; Ramírez-Díaz et al., 2023); or partially support them, as differentiation between NA and NCA corresponds to the ecological gradient, but for which sympatry in the southern IT (southern Oaxaca) between these same lineages may be (theoretically) unexpected.

Our phylogeographic tree showed further differentiation events within each group divided by the IT (Fig. 1). In the group west of the IT, within the NA clade, all individuals from Tres Marías islands (currently in subspecies capitalis) were grouped, showing a genetic structure consistent with simple paraphyly (Omland et al., 2006), in which lineage sorting seems to be incomplete. Simple paraphyly in this case may be the result of a small number of founders that likely dispersed to the islands (180–40 kya), which is consistent with islands emergence (McCloy et al., 1988) and with dispersal times in other insular bird species (Ortiz-Ramírez et al., 2018). In the group east of IT, genetic differentiation between sister clades YUC and NCA is relatively low (0.7%), also suggesting a recent divergence event. Genetic divergence between east and west Mesoamerican clades in L. verreauxi is very deep, underlying the role of the IT as a significant barrier for differentiation in lowland taxa. Values above traditional species thresholds for species recognition (Seutin et al., 1995; Hebert et al., 2004; Fuchs et al., 2021) were found between NA–NCA (3.8%), as well as between NA and YUC (3.6%). Our analyses also found deep differentiation in the SCA and the SA despite having a small number of samples.

4.2   Biogeography

Previous research in mid-sized New World doves (Leptotila, Geotrygon and Zenaida) suggested that Leptotila dispersed to North America from South America after the closure of the Isthmus of Panama about 3.4 Mya (1.9–5.8 Mya). According to Johnson and Weckstein (2011), this dispersal event allowed for both, the evolution of a northern clade (including all L. verreauxi taxa), and its radiation following a stepping-stone fashion back into South America. However, our results (Fig. 2) suggest two main differentiation events: the first involves North America and Central America + South America, while in the second one the earliest differentiation occurred in South America, thus suggesting that the Central American taxa (NCA and YUC) are the most recent diverged lineages. Our differentiation scenario differs from the stepping-stone biogeographic scenario suggested in Johnson and Weckstein (2011).

Based on our results, lineage differentiation in Mesoamerican L. verreauxi occurred during the Late Pliocene-Pleistocene (Fig. 2) at about 2.6 Mya (1.9–3.4 Mya) with the split of lineages west (NA) and east of the IT (NCA-YUC). This divergence event is consistent with climatic changes that led to vegetation changes across the region, likely severing connections between formerly continuous ecosystems and therefore, allowing persistence of suitable conditions within small patches (Fig. 3, see also Castillo-Chora et al., 2021), which reduced gene flow between lowland populations (e.g., Sánchez-González et al., 2021, 2023). Both The star-like pattern in the haplotype network and R₂ statistic for the NA clade support the biogeographic pattern described above, suggesting recent population growth probably due to the expansion of suitable conditions across the region. Demographic patterns in the other clades are unclear from our data.

Divergence date for the west and east groups is also consistent with proposed geographic processes in the region, such as formation of a seaway crossing the IT at approximately 2–3 Mya during the late Pliocene (Barrier et al., 1998), which has been frequently invoked as the main promoter of vicariant events in the region for the differentiation in unrelated lowland taxa, such as mammals (Whitmore and Stewart, 1965), geckos (Butler et al., 2023), and montane birds (Ornelas et al., 2013). The seaway hypothesis, however, has been challenged (Wyatt Durham et al., 1955; Mulcahy et al., 2006), and alternative hypothesis have been proposed, such as the ecological niche divergence (see Sánchez-González et al., 2023). Although the hypothesis related to a seaway across the IT may explain west-east divergence pattern, further differentiation events across the north to south ecological gradient have been described in other vertebrates, such as reptiles (Parkinson et al., 2000; Zarza et al., 2008; Gray et al., 2019; Butler et al., 2023).

4.3   Secondary contact

Our analyses revealed a secondary contact zone for the NA and NCA lineages south of the IT (Fig. 1), in which genetically divergent individuals are syntopic in the seasonally dry forests foothills of the Pacific slope of the Sierra del Tolistoque, Oaxaca. Sympatric lineages with a deep genetic divergence, such as these in L. verreauxi, have been rarely reported in birds (Hogner et al., 2012; Benham and Cheviron, 2019; Pons et al., 2019), but also in other groups (Tominaga et al., 2009; Xiao et al., 2012).

Different causes may result in the establishment of a secondary contact zone (Benham and Cheviron, 2019). Sympatric overlap of distinct lineages is typically interpreted as the result of formerly isolated divergent populations that resumed secondary contact (Hogner et al., 2012; Benham and Cheviron, 2019). Pleistocene climatic fluctuations may have promoted the divergence of the west and east groups across the IT, which may have been restricted to allopatric stable climatic areas (Fig. 3). After the disappearance of the barrier (either geographical or ecological), populations are expected to expand to ecologically suitable areas. However, environmental differences in the landscape and biological interactions (e.g., competitive exclusion) between taxa with similar niches may slow down the recolonization process (Anderson et al., 2002), acting as a reinforcement that may help in maintaining differentiation over time.

4.4   Ecological niche differentiation

Estimating ecological niches using intraspecific lineages or subspecies as units can represent their niche more accurately, without ignoring local adaptive responses as would occur when estimating models at the species level (Smith et al., 2019). Our comparisons showed that the ecological niches of the three L. verreauxi lineages are non-equivalent. However, our analyses detected niche similarity in the NA-NCA comparison, but niches were not similar for the clades found to the east of the IT. In the case of NA–NCA, both lineages are occupying different environments and climatic conditions, but ecological niches are similar, and therefore may evolve slowly (Peterson et al., 1999).

The similarity in the ecological niches of the NA and NCA may explain, in part, their sympatric distribution in the southern IT. For the clades eastern of the IT, comparisons in NA-YUC and YUC-NCA indicated that niches are not equivalent nor similar, suggesting niche divergence.

According to range-breaking tests, environmental differences across the north to south ecological gradient in the IT do not represent an abrupt environmental barrier for dispersal of the populations of L. verreauxi between habitats. Thus, our results suggest that the secondary contact zone in southern IT may have less suitable conditions for dispersal of NCA, than for NA. Thus, NCA populations are seemingly less tolerant to dry forest climatic conditions on Pacific slope of the IT, which likely prevents its dispersal further west. NA populations, however, seem to be tolerant to both humid and dry conditions, suggesting that other factors, such as the role of biological interactions at local scales (e.g., competition) between both lineages should be considered in maintaining this limited secondary contact zone. A further possibility is related to the seasonal differences in the precipitation between the dry forests of the southern IT and the mesic forests of the Soconusco region in coastal Chiapas (Binford, 1989): the seasonal increase of humidity would allow the dispersal of the mesic-adapted NCA further north to be in contact with NA in the southern IT. We suggest that this secondary contact zone, has been maintained through limited ecological suitability for NCA to the west of the IT, while to the east, divergence may have been maintained via ecological interactions or seasonal movements. We lack data on lineages north of the IT, where NA and YUC may overlap due to the absence of conspicuous barriers; a shallow genetic divergence however between these two lineages, suggests a corridor across the Usumacinta River drainage.

The Isthmus of Tehuantepec is an area long-recognized as a geographic limit for several North American bird families (Navarro-Sigüenza et al., 2014), as well as a promoter of speciation and ecological differentiation, particularly for montane bird taxa, which have different species divided by the Tehuantepec Valley (Barber and Klicka, 2010). However, its role as a barrier for lowland bird taxa is still poorly understood (see Binford, 1989). Patterns of divergence in lowland populations in the Isthmus region should be viewed with caution due to a dynamic biotic history, which due to changes in ecological conditions has likely promoted ecological niche differentiation, as suggested in several bird taxa (Castillo-Chora et al. 2021; Sánchez-González et al, 2021, 2023; Llanes-Quevedo et al., 2022).

Declaration of competing interest

The authors declare no conflict of interest.

Ethics statement

Specimens and tissues used in this project were collected under permits granted by the Secretaría del Medio Ambiente y Recursos Naturales (SEMARNAT) of the Mexican Government.

CRediT authorship contribution statement

Orlando J. Espinosa-Chávez: Writing – original draft, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. Adolfo G. Navarro-Sigüenza: Writing – review & editing, Supervision, Data curation. Hernando Rodríguez-Correa: Writing – review & editing, Supervision. Luis A. Sánchez-González: Writing – review & editing, Writing – original draft, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

Acknowledgments

This paper is part of the fulfillments of OJE-C master’s degree project in the Posgrado en Ciencias Biológicas, UNAM, who also had a CGEP Scholarship for UNAM postgraduate studies. Many thanks to the financial support for this study by PAPIIT- UNAM grant (IN222817) to LAS-G, and the François Vuilleumier Fund for Neotropical Bird Research from the Neotropical Ornithological Society (NOS) awarded to OJE-C. Special tanks to Uri O. García-Vázquez and David A. Prieto-Torres for comments to early drafts. We appreciate comments by two anonymous reviewers, which greatly improved our paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.avrs.2024.100160.