Understanding the role of ecological divergence in the evolution of isolated populations in the Arremonops rufivirgatus species complex across Mesoamerica

Abstract

Background

The topographic complexity and wide range of environmental conditions of the Neotropical region have allowed the evolution of the most diverse avifauna in the world. Distributional patterns within this avian diversity mirror this complexity, and many species show allopatric distributions in environmentally continuous regions. Here, we used environmental variables and historical presence records to understand the evolution of the distribution of three isolated groups (Gulf, Pacific, and Yucatan Peninsula) of the Olive Sparrow (Arremonops rufivirgatus) species complex. We assessed the role of environmental factors underlying geographic distribution patterns in the complex based on ecological niche modeling and performed paleoclimatic reconstructions to assess distributional changes based on suitable areas during the Late Pleistocene.

Results

Niche similarity was not rejected in the Pacific/Yucatan comparison, but the Gulf/Pacific and Gulf/Yucatan comparisons showed niche differentiation. We found regions with low climatic suitability representing a biogeographic barrier for the Pacific and the Yucatan groups, but not for the Yucatan and the Gulf groups, suggesting that biotic factors, such as competition with ecologically similar species, may be involved in geographic isolation.

Conclusions

Our results suggest that allopatric distributions in the three groups within the A. rufivirgatus complex probably evolved due to biotic interactions with ecologically similar species in the relatively environmentally continuous areas across the Gulf Slope, but to range contractions leading to isolation in the Yucatan and the Pacific groups.

Background

The Neotropical region holds the most diverse avifauna in the world [1, 2]. Some of the factors driving this species’ richness are related to topographic complexity and the geological history of the region (e.g [3]).,. These features have configured a complex range of environmental conditions in which species may differentiate due to macrohabitat preferences and climatic tolerances [4], but at the same time providing them with opportunities for dispersal and colonization [5].

Despite being the highest in the world, avian Neotropical diversity is likely underestimated (see [6]), especially due to the difficulties in establishing species limits and phylogenetic relationships, and, until recently, an almost exclusively qualitative-based assessment of environmental niche evolution and local adaptative responses [7]. Recent integrative studies, including new theoretical and quantitative tools for delimiting taxonomic units and niche evolution, have provided valuable information for a better understanding of the evolutionary patterns underlying biodiversity [8, 9].

Although geographic speciation is likely dominant in Neotropical birds [10, 11], the quantitative study of climatic niche evolution in geographical isolation and its relationship with species diversification under a niche modelling framework, has been relatively recently addressed [5, 12]. In closely-related species, niche conservatism may remain after geographic isolation [13, 14], since it is expected that isolated populations retain ancestral environmental conditions due to shared ancestry, therefore showing a limited fitness outside the ancestral niche [15, 16]. However, populations may shift to different ecological niches over time, resulting in adaptation to local ecological conditions [17,18,19]. This scenario implies both differences in niche characteristics due to environmental divergence among the accessible areas for each species, as well as the maintenance of isolated distributions due to unsuitable areas among populations (e.g [14, 20, 21]).

Species Distribution Models (SDMs) and Ecological Niche Models (ENMs) are essential tools in understanding how species are distributed across the geographic and environmental space, as well as in the role of the evolution of environmental niches in the speciation process [12, 20]. Modern approaches involving ENM include the quantifying of the environmental niche overlap between taxa, allowing us to understand the significance that ecological and geographical factors have in biotic diversification, as well as the assessment of biogeographic boundaries separating closely related taxa [22]. Recent research in Mesoamerican birds has shown the relationship between lineage divergence and ecological heterogeneity in different geographic settings. In montane birds of the Chestnut-capped Brushfinch (Arremon brunneinucha) complex, most lineages showed low climatic niche overlap, suggesting niche divergence in allopatric populations (e.g [20]), while in lowland birds, historical changes in the climatic conditions along the seasonally tropical dry forest in the Pacific slope led to the geographic restriction that probably promoted the adaptation and differentiation in different bird species (e.g [23, 24]).

The Olive Sparrow (Arremonops rufivirgatus; Passerifomes, Passerellidae), is a sedentary Neotropical bird with geographically isolated populations across the lowlands of Mesoamerica (Fig. 1). Based on subtle plumage variations and an allopatric distribution [25, 26], some taxonomic proposals group their populations in putative evolutionary lineages, recognizing two (rufivirgatus, superciliosus; [27, 28]) or three species with a disjunct distribution (rufivirgatus, from extreme SE Texas to the northern Isthmus of Tehuantepec in the Gulf slope; superciliosus, from Sinaloa in western Mexico to northwestern Costa Rica in the Pacific slope; and verticalis from the Yucatan Peninsula; [29,30,31,32,33]). More recently, the study of the acoustic divergence in the songs and behavioral responses also supported three groups, suggesting behavioral isolation [34]. Differentiation may be expected due to the widespread but allopatric distribution of this species complex, which covers several climates and habitats, implying that environmental variation may have played a role in the divergence and evolution of these groups.

Fig. 1

Occurrence records for three studied groups within Arremonops rufivirgatus complex used to perform the Ecological Niche Models (ENMs). Colored dots represent each group, purple = rufivirgatus (Gulf), red = superciliosus (Pacific), light orange = verticalis (Yucatan). Two biomes are shown on the map, representing the Tropical Dry Forests (orange) and the Tropical Moist Forests (green)

Here, we assessed the significance of niche environmental variables in shaping the historical and current geographical distribution of three isolated groups within the Arremonops rufivirgatus complex through the development and analyses of SDMs. We also hypothesized ecological niche differentiation due to the geographic disruption of suitable environments among these groups. Finally, we discuss the implications of environmental drivers and geographic patterns as potential reinforcers of the divergence of isolated populations in this species complex.

Materials and methods

Defining the units of analysis

Given that the geographic distribution of the Olive Sparrow spans three major lowland Mesoamerican regions with contemporary ecological and geographical continuity [26, 33], we used these areas for our ecological analyses. This criterion is further supported and congruent with behavioral isolation [34] and taxonomical proposals recognizing three groups [29]. Thus, in this study, we considered three isolated geographic groups within the Olive Sparrow complex: the Gulf group (including rufivirgatus, ridgwayi and crassirostris), the Pacific group (including superciliosus, sumichrasti, sinaloae and chiapensis), and the Yucatan group (including verticalis and rhyptothorax).

Occurrence data gathering

We compiled a database of all available records of A. rufivirgatus in (1) the ‘Atlas of Distribution of Mexican Birds’ [35, 36]; (2) The Global Biodiversity Information Facility (GBIF, 2023; available in: http://www.gbif.org); (3) observational occurrences in eBIRD (The Cornell Lab of Ornithology-The Audubon Society, ebird.org); and (4) sound records from private and public acoustic libraries, such as Macaulay Library, Cornell Laboratory (http://macaulaylibrary.org), Biblioteca de Sonidos de Aves del Museo de Zoología ‘Alfonso L. Herrera’ (Facultad de Ciencias, Universidad Nacional Autónoma de México), Banco de Sonidos de las Aves de México BSAM (Instituto de Ecología, A.C., México), Xeno-Canto (Xeno-canto Foundation, http://xeno-canto.org), and Sounds of Nature SONAT private collection (A. Celis-Murillo). We complemented all this information with records from our field surveys carried out in México from 2016 to 2020 [34, 37] and 2022 to 2023. In total, we obtained 129,212 occurrence records.

We performed a data cleaning procedure on all compiled records through manual inspections and software tools. We used OpenRefine (http://openrefine.org) for cleaning and transforming databases from different sources, which allowed us to fix inconsistencies and merge and match similar values. We also validated the information corresponding to date records and georeferencing accuracy. For this latter, we visualized the congruence of coordinates with the geography using QGIS v.3.28 [38]. We used the Species Matching tool (https://www.gbif.org) for taxonomic validation, normalization, and correction of the binomial nomenclature and synonymies based on the GBIF Backbone Taxonomy.

We verified the correspondence between the historical range and the species distribution models using the Mesoamerican birds’ digital maps set [39], as well as in the species distribution maps of the birds of the world [40]. For records in citizen science platforms (e.g., eBird), we retained only those supported with photographic, video, or audio evidence. We removed ambiguous information related to the locality and records with problematic or inaccurate species’ occurrence. We also followed the BirdLife International species datasheet [41] to remove those records beyond the known upper elevation limits of each group, as well as those records without information on climatic data.

We further performed spatial and environmental data cleaning using climatic and elevation layers (see Environmental data). We retained occurrence records within the timeframe covered by the climatic variables and evaluated the climatic suitability to include records from 1979 to 2013 (see Environmental data). Additionally, for locality records from 2014 to 2022, we performed an outlier exclusion procedure in the environmental space outside of the interquartile range of three environmental variables, representing the annual average values of temperature and humidity (annual mean temperature [bio 01], annual precipitation [bio 12]), and precipitation seasonality [bio 15] for occurrences from 2014 to 2022 [42]. This procedure allowed us to both reduce the sampling bias and improve data quality [43]. These steps are important to identify problematic or imprecise species’ occurrences with incorrect climatic values because the choice of a climatic baseline and the reduction of sampling biases influence model performance [44, 45].

Finally, we used “spThin” [46] for R [47] to correct for spatial biases due to the high density of occurrences in some areas (see [13]). For this step, we estimated the most appropriate distance among the presence records and the nearest neighbor to retain only those records separated by at least 1 km, which also corresponds to the spatial resolution of environmental layers used in the analyses. After these steps, we obtained a total of 5,491 occurrence records for the three groups, which are divided as follows: Gulf group = 3,543; Pacific group = 1,151; and Yucatan group = 797. All geographic coordinates were transformed to decimal degrees based on the WGS84 datum.

Defining the accessible area (M)

The concept of the area ‘M’ (defined as the historical region accessible to dispersal or colonization by the species over some relevant time interval [13, 48]) is considered as an important step in ecological niche modeling, since it supports the outcome of calibration, evaluation, and comparison of models [49]. Therefore, we defined M according to the BAM (biotic and abiotic environmental conditions, and movement) diagram [13, 48]. To do this, we delimited M area for each modeled group by the overlapping of occurrence records with biogeographic regions and terrestrial ecoregions [50]. We assumed that the M area may be reached by dispersal from existing populations, thus probably representing the species’ tolerance limits and, therefore, the historical and ecological barriers to dispersal. Based on previous studies with marked individuals reporting limited movements and high natal philopatry in A. rufivirgatus [33, 34, 51], we draw a buffer of 13 km around the area M of each group to account for potential dispersal events (Fig. 2).

Fig. 2

Representation of the geographical and environmental space in the Olive Sparrow species complex. Maps show pairwise comparison and overlapping area (diagonal line patterns) between each group’s estimated distributions. A transparent colored polygon represents the accessible area (M) of each group. Plots along the diagonal (A, E, I) represent occurrence density in the environmental space for each group, with the color gradient indicating the suitability of environmental conditions. The shared area (M) is represented by a continuous line in each map (the dashed line represents 50% of the available environment)

Environmental data

We used bioclimatic and vegetation variables to characterize the environmental niches for each group. We downloaded 19 bioclimatic variables from CHELSA v.1.2 [52] with a spatial resolution of 30 arc sec (≈ 1 km²) describing temperature and precipitation, seasonality, and other limiting environmental factors in a period ranging from 1979 to 2013. Additionally, for the ecological niche comparison in the present, we used a total of four vegetation variables related to habitat heterogeneity based on the textural features of the Enhanced Vegetation Index (EVI), acquired by the Moderate Resolution Imaging Spectroradiometer (NASA-MODIS/Terra; available at http://earthenv.org). The EVI metrics obtained were: (1) coefficient of variation, which measures the normalized dispersion; (2) standard deviation, which measures the dispersion; and (3) dissimilarity, which measures the difference between adjacent pixels. We also used a layer of percent tree cover (PTC) representing the density of trees on the ground (available at https://globalmaps.github.io/ptc.html). Finally, to avoid collinearity biases due to the inclusion of redundant environmental variables, we selected and retained 16 variables (Table 1) based on the Variance Inflation Factor (VIF < 10), as implemented in ‘usdm’ [53] for R.

Table 1 Environmental variables retained and used to perform ENMs in each group of the Arremonops rufivirgatus complex. For each group, we show the variance inflation factor (VIF value). Shared variables among the three groups in bold

Geographic spatial projection from ecological niche modeling

We obtained potential distribution models for each group using a maximum entropy algorithm, as implemented in MaxEnt v.3.4.4 [54], which calculates the most probable distributional range of a taxon based on the occurrence (presence-only data) and the environmental variables of the accessibility area. We followed a calibration protocol to configure a set of parameters to access model complexity, and evaluated and selected the best model based on the omission rate range using ‘kuenm’ v.1.1.9 [55] for R. We defined the following arguments for calibration: 15 regularization multipliers (0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.5, 3, 3.5, 4, 5) and feature classes (l = linear, q = quadratic, p = product, t = threshold, and h = hinge) in the basic combination (“l”, “lq”, “lqp”, “lqpt”, “lqpth”). This configuration generated 75 candidate models per group, which were evaluated via the partial receiver operating characteristic (partial-ROC; [56]) by analyzing the proportion of 10 bootstrap replicates with the area under the curve (AUC-ratio). We selected all statistically significant models meeting the omission rate criteria (less than 5%) and with the least complexity based on the Akaike information criterion corrected for small sample sizes (AICc).

We then created presence/absence maps (continuous probability from 0 to 1; [54]) from the logistic values of the suitability maps, using the minimum training presence (MTP) as the threshold value [57]. To test the environmental equivalence among groups, we used binary maps for inter-prediction (degree of geographical matching) of the different group models. We expected that if the environmental niche is similar among species, models should forecast similar potential distribution ranges despite barriers to dispersal (see [7]). Later, for each group, we described the environmental profiles that defined the potential areas and calculated the ranges for each of the variables considered in sites where models predicted areas currently suitable [58]. Then, we extracted the values of the environmental variables for each group and performed a Discriminant Function Analysis (DFA) using the ‘MASS’ package in R [59] to assess the conditions defining how suitability areas differ among groups.

Niche overlap in the ecological space

To test niche similarity/differentiation among the three groups in the Olive Sparrow complex, we performed a Principal Component Analysis (PCA) using the retained environmental variables shared among the three groups and analyzed their ecological overlapping using PCA-env, an ordination technique quantifying niche overlap calibrated on the entire environmental space occupied by the studied taxa [60]. We then calculated pairwise niche overlap values using Schoener’s D index [61], which measures the similarity of suitable habitats for a pair of models by assessing the overlap along the principal component scores discriminating differences between the environmental space occupied by the studied taxa. The index ranges from 0 to 1 and is interpreted as follows: 0–0.2, no or very limited overlap; 0.2–0.4, low overlap; 0.4–0.6, moderate overlap; 0.6–0.8, high overlap; and 0.8–1.0, very high overlap [62]. As in the A. rufivirgatus groups, species showing fragmented distributions may have different accessible environments that may reflect a different environmental distribution unrelated to their ecological preferences [12]. Therefore, we performed a similarity test to correct for the effects of differential suites of environments available for each group on the niche differentiation [17]. We then tested the null hypothesis of niche similarity compared to background environments. From this perspective, niche similarity is accepted if the observed D values are significantly greater than expected from 1,000 pseudoreplicated data sets, allowing random shifts [17, 60].

Ecological barriers and paleodistributions

To test whether species’ range limits are bounded by abrupt environmental conditions that may act as ecological barriers, we performed the ribbon range-breaking test, as implemented in ENMTools [63], which allowed us to explore whether the ranges of two taxa are divided by unsuitable habitats representing an environmental transition within the drawn limits [22]. We generated 100 pseudo-replicates of the ribbon for this analysis through the combined range of each pairwise group comparison. We kept constant the width of the ribbon (about 70 km), estimated using QGIS v.3.28 ( [38]) based on the width of the continuity and overlapping area predicted in the distribution models among each paired comparison in the three Olive Sparrow groups. We located the first ribbon at the Isthmus of Tehuantepec (Pacific/Gulf comparison) and a second ribbon at the base of the Yucatan peninsula (Gulf/Yucatan and Gulf/Pacific comparisons) (Fig. S1). We calculated Schoener’s D and Warren’s I (a measure to quantify the dissimilarity between two probability distributions based on the Hellinger distance) values from the pseudo-replicates to generate null distributions for comparison with empirical values. In the ribbon range-breaking test, ENMs were generated for each set of occurrences on each side of the line and the ribbon using a maximum entropy algorithm.

Finally, to identify the historical range extent and stability areas of each group, we performed a paleoclimatic reconstruction of the distribution of suitable areas during the Late Pleistocene climatic fluctuations using the same methodological approach as for the ENM (see above). We used the environmental variables from CHELSA v.1.2 [52] with a spatial resolution of 2.5 min (~ 5 km²) to perform the niche models at the present. We projected models to suitable climatic conditions during the Mid-Holocene (Mid-Hol, 6000 Ya) and to the Late Pleistocene climate fluctuations, including the Last Glacial Maximum (LGM; 21,000 Ya) and the Last Interglacial (LIG, ~ 120,000–140,000 Ya). Since vegetation variables are not available in paleo-temporal scales, we used the retained climatic variables shared among the three groups (bio 03, bio 08, bio 14, bio 18), which we downloaded from PALEOCLIM (Paleclim.org [64]), with a spatial resolution of 2.5 min (~ 5 km). We created presence/absence maps (from the logistic values of the suitability maps) using the minimum 10-percentile training presence as the threshold value. To assess the stability and reliability of the transference of the models to the past, we performed an extrapolation risk analysis via the Mobility-Oriented Party (MOP) test as implemented in “NicheToolbox” [65]. This procedure allowed us to identify regions where environmental conditions in past climates deviated significantly from those in the calibration period, thus flagging areas where model predictions may be less reliable [66]. We performed the analyses considering random sampling of 25% for each area M.

Results

Present geographical distributions

We selected a set of significant environmental variables for model building for each group (Table 1). Overall, selected final models were used to perform 10 replicates, which showed significant partial ROC test values (> 1.5, Table 2A).

Table 2 Summarized model selection results for (A) climate-vegetation current distribution and (B) climate paleodistribution

Based on the MTP threshold and our binary maps, the largest geographic range was estimated for the Gulf group, which covers the lowlands of the Gulf of Mexico and is limited by the elevations of the Sierra Madre Oriental; a clear bottleneck is shown in the northeastern slope of the Trans-Mexican Volcanic Belt in central Veracruz. According to our geographic model, currently non-occupied but suitable areas were highlighted at the northern and central Isthmus of Tehuantepec, as well as in the southern Yucatan peninsula (Fig. 2A, B). The potential distribution estimated for the Pacific group runs almost continuously from northwestern Mexico (southern Sinaloa), where the predicted distribution is limited by the highlands of the Sierra Madre Occidental and Sierra Madre del Sur, across the Pacific slope to the lowlands of Costa Rica. There is a distributional gap likely due to the humid lowlands of southeastern Mexico (Soconusco) and southern Guatemala (Fig. 2C, D). The smallest potential distribution range was estimated for the Yucatan group, which covers the Yucatan peninsula and small and isolated areas in the central Chiapas depression (Fig. 2E, F).

We calculated the geographical distribution overlap as an indirect estimation of the ecological niche overlap. In general, we found that predicted distributions showed a limited (< 17%) geographical matching among groups (Table S1). The three groups showed a potentially suitable geographical distribution overlap in the southern Yucatan Peninsula; suitable areas for two of the groups were in the central Isthmus of Tehuantepec and central Veracruz and northern Oaxaca (Fig. 2). We found the highest geographic overlap between the Gulf and the Yucatan groups (17%, Fig. 2E), in contrast, the smallest overlapping area was between Pacific and Yucatan (3%, Fig. 2D).

Additionally, only seven of the environmental variables retained based on the VIF criteria were shared by all three species groups. These are related to variation in the ranges of temperature, precipitation, and vegetation (Table 1). The DFA showed high accuracy in the classification of the different Arremonops groups, suggesting that conditions defining the suitability areas among groups differ (0.821, 95% CI = 0.820,0.821, Kappa = 0.7, P < 0.001; Table 3, S2, Fig. S2).

Table 3 Confusion matrix for prediction of Arremonops groups based on the environmental ranges of shared environmental variables

Niche overlap in the ecological space

The PCA analysis indicated that the two first principal components explain 54.72% of the environmental variation (Fig. S3). The first principal component (PC1) explained a 31.61% of the total variation and was strongly associated with the habitat structure (percent of tree cover [tree_ptc]) and the relation of annual and daily temperature oscillations (isothermality [bio_03]), while the second principal component (PC2) explained a 23.11% of the total variation and was strongly associated with the canopy structure and its variation (coefficient of variation of Enhanced Vegetation Index [evi_cv] and dissimilarity of the Enhanced Vegetation Index [evi_dis]).

Our results indicated no or very limited ecological niche overlap (Schoener’s D < 0.2, Table S3) for each paired comparison among geographic groups. According to the background tests in the Gulf/Pacific and Gulf/Yucatan comparisons, results indicated a low niche similarity, suggesting that the ecological niche in each group is not more similar than random expectations in both directions. Such results allowed us to reject niche similarity among these groups. However, we could not reject niche similarity in the Pacific/Yucatan comparison (Fig. 3).

Fig. 3

Results of the background niche similarity test for each paired comparison among groups, measured by Schoener’s D. Histograms correspond to the expected distribution based on the randomization of the data. A vertical red line with a diamond represents niche overlap values compared to a null distribution of background divergence

Biogeographic boundaries and paleodistributions

Schoener’s D values from the ribbon range-breaking test suggest a region with unsuitable environmental conditions separating the Pacific group from the Yucatan and the Gulf groups (P < 0.05; Table 4; Fig. 4), but not among the Gulf and Yucatan groups (P = 0.139). Additionally, the environmental conditions of the ribbon in the flanking regions are no more different from one another than expected by chance for all the paired comparisons (P < 0.05, Table 4; Fig. 4).

Table 4 Schoener’s D values and significance of ribbon range-breaking test

Fig. 4

Results from the ribbon range-breaking test between pairwise groups. Vertical bar plots show values of Schoener’s D in flanking regions. Significance is shown as ns = non-significant and *=P < 0.05. Vertical density plots indicate the distribution of Warren’s I values from 100 unique pseudo-replicates calculated between the flanking regions and the ribbon

All paleodistribution model projections showed high performance (ROC partial > 1.4, Table 2B), indicating that models are significantly better than expected by chance. High environmental suitability areas for the three groups in the Olive Sparrow complex showed changes in size and location in comparison to the present (Fig. 5). Thus, for the Gulf group, LIG indicated widespread suitability across the central and southern Gulf slope from extreme southern USA and northwestern Mexico to Panama. Climatic changes produced range contractions during the LGM, with suitable areas mostly restricted to the Gulf Slope and the northern Isthmus of Tehuantepec, across the Pacific slope and the northeastern Yucatan peninsula. Later, in the Mid-Holocene, the geographic range expanded towards the northern Gulf slope, most of which has been maintained up to the present.

Fig. 5

Presence-absence maps for the three Olive Sparrow groups based on ecological niche models projected onto the geographic areas for the Gulf group (violet), Pacific group (red), and Yucatan group (orange). Maps depict suitable historical distribution areas for climate niches across the Present, the Holocene, and Last Glacial Maximum (LGM) and Last Interglacial (LIG) periods of the late Pleistocene. For the Present, maps show the distribution within the accessible area. Darker colors depict regions of historical climatic stability in a pairwise comparison with the Present. Smaller areas for Yucatan group in LGM are pointed with arrows

For the Pacific group, LIG showed suitable environments isolated in three main regions throughout the Pacific slope: (1) from northern Sinaloa to eastern Guerrero in Mexico, (2) the southern Isthmus of Tehuantepec in southeastern Mexico, and (3) from southern Guatemala to western Nicaragua. This widespread distribution severely contracted to highly isolated habitat patches across the Pacific slope and the southern Yucatan peninsula. Most of these isolated patches disappeared during the Mid-Holocene, and suitable environmental conditions were restricted to northern Mexico and contracted even to smaller regions across the Pacific slope. Suitable areas in the present resemble the distribution during the LIG, suggesting recolonization from these small patches.

For the Yucatan group, LIG showed isolated suitable areas in Mesoamerica, which severely contracted in the LGM, maintaining small and isolated regions in the base of the Yucatan peninsula, southern Central America (probably in the Volcanic Arc in El Salvador), and in the Pacific slope (southern Sinaloa). However, suitable environmental conditions in the Mid-Holocene recovered in the Yucatan Peninsula, from which this taxon expanded to most of the areas in the present. In general, according to environmental suitability maps, the historically suitable areas for all three groups appear as mainly allopatric, except in some parts of the Yucatan Peninsula in different time intervals (Fig. 5).

Finally, our MOP analysis showed high congruence in the climatic conditions among past and present ENM models but showed some regions with non-analogous climatic conditions for all groups during the LIG (e.g. northern Central America and the southwestern Mexico; Fig. S4).

Discussion

Species’ distributional patterns and their environmental tolerances allow the understanding of the role of ecological factors in speciation and niche evolution [67]. In polytypic taxa, using lower taxonomic rank categories to estimate potential distributions and ecological niche similarities (or differences) has been recommended as a more accurate approach [7]. This approach allowed us to find niche divergence in the three groups analyzed within the Olive Sparrow complex, suggesting that Grinnellian niches are not conserved, awarding ecological differentiation among groups.

The niche conservatism hypothesis predicts that the distribution of closely-related species is limited by contrasting environments due to physiological constraints [68, 69]. This pattern may be expected as closely related taxa share ancestral environmental suitability (see [13, 70]). Nevertheless, our similarity tests indicated that there are niches with low environmental similarity for most of our pairwise comparisons. Although we found low similarity in the ecological niches of the Pacific and the Yucatan groups, we cannot reject the niche conservatism hypothesis [12, 60]. Niche conservatism has been found in other birds with similar allopatric distributional patterns in seasonally dry tropical forests in the Yucatan peninsula and the Mesoamerican Pacific slope, suggesting that currently disjunct geographic distribution may be due to range contraction [14]. However, the Yucatan group inhabits a wider range of environmental conditions, covering both dry and moist forests in the Yucatan peninsula, suggesting that the occupation of some habitats in the region may be due to ecological release, probably promoted by niche differentiation and further expansion [71]. Thus, whether the lack of significant niche similarity between the Yucatan and the Pacific groups may be due to a widespread ecological distribution or to historical range contraction remains to be tested.

Overall, the geographic distribution for the three groups in the Olive Sparrow species complex mirrors emergent evolutionary patterns of the Mesoamerican lowland forests related to the historical climate changes that may have promoted the evolution of different lineages in these biomes (e.g [14, 23, 72, 73]).,. Despite the fragmented distribution and differences in environmental conditions, the projection of the niche models on each other’s geographic distribution showed some overlapping areas, as in the northern Isthmus of Tehuantepec (Gulf/Pacific; Fig. 2) and the southern Yucatan peninsula (Gulf/Pacific, Gulf/Yucatan, Pacific/Yucatan; Fig. 2). The current absence of occurrence records in these areas may be associated either with historical isolation among populations due to environmental barriers (see below), or to the role that biotic interactions may have in maintaining and reinforcing geographic isolation and divergence among groups.

The niche divergence found between the Gulf and the Yucatan groups suggests local adaptation and tolerance to different environments despite the apparent geographic continuity of the lowlands. Similar scenarios have been proposed in codistributed bird taxa, such as the Buff-bellied hummingbird Amazilia yucatanensis [74], and the White-tipped Dove Leptotila verreauxi [24], which showed concordant phylogeographic breaks in the southern Gulf of Mexico slope, separating populations between northern (Gulf slope) and southern (Yucatan peninsula) groups with non-equivalent niches. Our tests on biogeographic boundaries revealed environmental differences between most groups due to the potential presence of unsuitable areas in the southern Yucatan peninsula and the Isthmus of Tehuantepec. This latter region seems to act as a barrier for the Pacific population due to environmental differences found in the comparisons between the Gulf/Pacific and Pacific/Yucatan groups, which may limit dispersal in these populations. However, the area in the southern Yucatan Peninsula apparently does not represent a barrier of unsuitable climatic conditions for the Gulf and Yucatan groups. Additionally, the barrier did not show differences with the involved taxa on both sides in all comparisons, suggesting that probably biotic factors have maintained the isolation in these groups. A similar scenario has been reported in three Leptotila lineages in Mesoamerica, in which the environmental differences among lineages and the presence of a potential unsuitable region in the Isthmus of Tehuantepec have likely promoted secondary contact [24]. Our results may also indicate that other factors such as biological interactions with closely related taxa with similar niche requirements, may have a significant role in maintaining differentiation and isolation in the geographically isolated Olive Sparrow populations.

According to our results, ecological differentiation has maintained the allopatric distribution in the Olive Sparrow species complex. Based on phylogenetic analyses in New World Sparrows [75, 76], the Arremonops northern clade (A. rufivirgatus and A. chloronotus) diverged relatively recently (~ 3.03 Ma) in comparison to the clade in southern Central and South America (A. conirostris and A. tocuyensis). These results suggest that the northern Arremonops clade colonized new regions during the Late Pliocene and Pleistocene climate cycles. Late Pleistocene climatic suitability indicates that the distribution of the Gulf group has expanded and maintained since the Mid-Holocene to northeastern Mexico and Texas. The vanishing suitable areas detected in southern México (mainly in the Yucatan Peninsula) and northern Central America may be due to biological interactions involving a potential competitive exclusion with A. chloronotus and A. conirostris a local scale. The dynamic biotic history between these taxa and their related environmental preferences may sustain differentiation over time, thereby delaying and limiting opportunities for colonization, similar to findings in other species [77]. The high fragmentation detected in suitable areas across the Pacific slope since the LIG suggests that the Pacific group should be evaluated in more detail, as some populations in this group, which includes four apparently isolated subspecies, may be restricted by ecological barriers, the competitive interactions, or geographic isolation as has been found in similarly codistributed bird taxa [21, 23, 73, 78, 79].

Subsequent climatic changes and stable climatic areas during the LGM may have promoted opportunities for allopatric diverge in isolation, as in other bird lineages (e.g [78].; Columbidae: Leptotila verreauxi [24], Picidae: Melanerpes santacruzi [80]; Tyrannidae: Attila spadiceus [81]; Cuculidae: Piaya cayana [23]; Corvidae: Calocitta [21]. A detailed reconstruction of the historic paleoclimatic distribution in the Olive Sparrow and the Arremonops genus with support of multidisciplinary approaches using molecular, genetic, and morphological analysis is needed to understand the historical and ecological process associated with historical climatic stability that may support these divergence patterns.