Niche conservatism and convergence in birds of three cenocrons in the Mexican Transition Zone

Study area

The MTZ includes the mountain areas of Mexico, Guatemala, Honduras, and El Salvador (Halffter, 1987; Halffter & Morrone, 2017; Morrone, 2020; De Mendonça & Ebach, 2020). We used a polygon of the MTZ (Fig. 1, left) that includes the Sierra Madre Occidental, Sierra Madre Oriental, Trans-Mexican Volcanic Belt, Sierra Madre del Sur, and Chiapas Highlands biogeographic provinces (Morrone, 2020), and adjacent lowlands. In the mountain bird database, Quintero & Jetz (2018) named this area as the Mexico and Sierra Madre de Chiapas mountain ranges.

To understand the environmental conditions within this area (also called environmental space), we performed a principal component analysis to the bioclim layers from WorldClim 2.0 (Fick & Hijmans, 2017) with a resolution of 2.5 arc min. The details of the analysis and results of the PCA are available in Appendix 1. The first two components explain 80% of the environmental variation (56.8%, and 23.2%, respectively, Table 2) and were used to visualise the existing climate conditions in the MTZ (Fig. 1, right).

Phylogenetic data

To include the relationships among the bird species for the species selection and phylogenetic diversity measures we used Jetz et al.’s (2012) dated phylogeny with Hackett et al.’s (2008) backbone. The original phylogeny includes all the 9,993 species described by the time of its publication (available at https://birdtree.org/downloads). We used the last set of stage two trees with Hackett’s backbone, which includes 1,000 trees to incorporate phylogenetic uncertainty. We pruned the phylogeny to keep only the selected species, listed in Appendix 2 and plotted in Fig. S1. All the analyses were run through the 1,000 trees set to incorporate phylogenetic uncertainty.

Species selection

The process of species selection involved several steps that resulted in a total of 10 orders, 47 families, 326 genera and 808 species (see Supplementary Material Appendix 2). Most of these species are distributed in the MTZ and Mesoamerica (286 species) or have more than 80% of their distribution in the MTZ (246 species). We removed the following species: (1) those of the aquatic bird orders (Suliformes, Phaethontiformes, Procellariiformes, Charadriiformes, Anseriformes, Podicipediformes, Gaviiformes, Ciconiiformes, Pelecaniformes, Gruiformes, and Eurypygiformes), for not being representative of the Mexican avifauna; (2) migratory families (Parulidae, Hirundinidae, and Motacillidae), since they have an additional complexity due to the niche-switching and niche tracking phenomena (Nakazawa et al., 2004); (3) widespread families (Accipitridae, Cathartidae, Columbidae, Falconidae, Tytonidae and Strigidae), because wide-ranging species drive the variation in richness or average species attribute values (Quintero & Jetz, 2018); and (4) families not distributed in the Americas (Alaudidae, Certhiidae, Otididae, and Prunellidae), for not being relevant for this work.

Then, we selected species listed in the MTZ mountain ranges (Mexico and the Sierra Madre de Chiapas mountain ranges, (Quintero & Jetz, 2018) whose distribution is certain and have native distribution (breeding season and resident) according to the Handbook of the Birds of the World (BirdLife International, 2020). Then, we identified the sister species or clades of each of the species distributed in the study area by finding the parental node and its descendants in all the trees to ensure that only monophyletic groups are analysed using Jetz et al.’s (2012) dated phylogeny with Hackett et al. (2008). All the analyses were repeated through the 1,000 trees included in the set to incorporate phylogenetic uncertainty. Finally, we pruned the phylogeny to keep only the selected species, listed in Appendix 2 and plotted in Fig. S1.

Cenocron assignment

The original concept of cenocron was based on the differential richness of clades as a proxy of the ancestral area (Juárez-Barrera et al., 2020). So, we used Morrone’s (2014) regionalization for the Americas and the ecoregions’ realms for the rest of the world, recoding these areas as the northern hemisphere (Nearctic, Palearctic, or Saharo-Arabian), the southern hemisphere (Afrotropics or Neotropics), or Mesoamerican (Mexican Transition Zone and Mesoamerican dominion). To indicate shared evolutionary history, we assigned each family to a cenocron using these rules: (1) Nearctic if the family has >80% of its species distributed in the Northern hemisphere or if it has more species in the Northern hemisphere than in the Southern hemisphere or Mesoamerica, and phylogenetically related species are Nearctic; (2) Mountain Mesoamerican if a family has more species in Mesoamerica than in the Southern and Northern Hemispheres, and if phylogenetically related species are Neotropical; and (3) Typical Neotropical, if a family has >80% of its species distributed in South America and phylogenetically related species are Neotropical (Fig. 2). A detailed description and code of this section are included in Appendix 1.

Collection, cleaning, and thinning of occurrence data

We obtained the occurrence data for the years 1950 to 2020 from CONABIO (2021) and GBIF.org (2021a); GBIF.org (2021b); GBIF.org (2021c); GBIF.org (2021d); GBIF.org (2021e); GBIF.org (2021f); GBIF.org (2021g); GBIF.org (2021h). The occurrence query in GBIF included native occurrence data from specimens and observations, marked as without spatial issues, with less than 10 km of uncertainty and a coordinate precision >0.01. Since the data set was too large, it was split into several data sets that can be found on GBIF.org (2021a); GBIF.org (2021b); GBIF.org (2021c); GBIF.org (2021d); GBIF.org (2021e); GBIF.org (2021f); GBIF.org (2021g); GBIF.org (2021h). The database was cleaned using the CoordinateCleaner R package (Zizka et al., 2019) to automatically remove duplicated records, records at sea, and near geographic centroids. We also omitted records outside known breeding distribution according to the Handbook of the Birds of the World (BirdLife International, 2020). Species with more than 5,000 records were down-sampled to 5,000 records and then thinned the data to ensure >5 km of distance between records using the spThin R package (Aiello-Lammens et al., 2015). This distance ensures a minimum separation of two cells between records in the 2.5 arc-minute raster layers utilized in our modelling, minimizing the potential for record clustering and avoiding overfitting. To construct the model, we split the data into training (75%) and testing (25%) for the model construction and validation.

Environmental data and calibration area

We included both ecological and historical conditions to build species distribution models. As ecological factors, we selected variables of high biological importance from WorldClim 2.0 (Fick & Hijmans, 2017) and ENVIREM (Title & Bemmels, 2018). The selected variables fall into three categories: (1) Topographic: elevation, slope, and topographic roughness index, the elevation data from WorldClim were used to create the slope and Terrain Ruggedness Index; (2) Moisture: Thornthwaite’s aridity index, climatic moisture index, and annual precipitation; and (3) Temperature: maximum temperature of the coldest month, minimum temperature of the coldest month, the maximum temperature of the warmest month, minimum temperature of the warmest month. We used a resolution of 2.5 arc min, which was the optimal suitable resolution for conducting a comprehensive macroecological analysis within an environmentally heterogeneous region, given the available computational resources. Each species was modelled with a subset of these variables assuring that the correlation between them was less than 0.5.

To include historical factors, climate layers were cropped to include only those areas that have been accessible to each species over a relevant time (Soberón & Nakamura, 2009). We obtained this area for each species, known as the accessible area or M (Barve et al., 2011), by masking the raster environmental data with the breeding occurrence records and the biogeographic provinces of the Nearctic region (Escalante, Rodríguez-Tapia & Morrone, 2021) and the Neotropical region (Morrone, 2014; Morrone, 2017; Morrone & Ebach, 2022), then added a 1 km buffer, this extra area is added to avoid edge effects and reduce the risk of biased predictions near the distribution boundaries.

Species distribution modelling and hyperparameter optimization

Species distribution models (SDMs) for each selected bird species were constructed to improve the distributional data predictive value and its integration with spatial analysis. We used the MaxEnt algorithm version 3.4.4 (Phillips, Anderson & Schapire, 2006) with the SDMTune package in R (Vignali et al., 2020). This package allows trainng and tuning the hyperparameters of Maxent in an automated way. For each species, we started with a “default model” (regularization multiplier = 1, feature classes = “lqph”, and 500 iterations). Then, we undertook a data-driven variable selection process to get a unique set for each species using the varSel function, selecting a correlation threshold of 0.5, and using AICc as metric to discriminate between models. This function iterates through all the variables starting from the one with the highest contribution and performs a Jackknife test among the correlated variables to remove those that improve the model when omitted. Five species could not be modelled using these hyperparameter tunning and the default model was used.

Then, we performed a grid search to select the best combination of hyperparameters for the model constructed with the selected variables. The parameters tested were regulator multiplier (from 0.8 to 1.2 by a 0.1 increase) and feature classes (“lq”, “lp”, “qp”, and “lqp”). To select the final model, we chose the top five models according to the delta AICc. Then we evaluated the models using the partial ROC (using ntbox Osorio-Olvera et al., 2020), and the omission rate at 5% and selected the model with the highest pRoc value. Evaluation results are included in Appendix 3.

The predicted presence map was obtained with a “clog log” output. To obtain a binary map, we used a Maximum test sensitivity plus specificity threshold (Liu, Newell & White, 2016), since it produces consistent results when using either presence/absence or presence-only datasets, and between rare and common species. This binary map was used for the subsequent analyses. Furthermore, we confirmed visually and manually that our model outputs matched the known distribution according to the Handbook of the Birds of the World (BirdLife International, 2020). A detailed description and code of this section is included in Appendix 1 and the description of the hyperparameters chosen for each species is included in Appendix 3.

Assemblage measures

We calculated various assemblage measures for each cenocron to understand their differences and correlations with environmental variables using the distribution models to calculate raster data of alpha diversity (species richness), age (Mya), number of nodes, phylogenetic diversity (PD), and weighted endemism. Alpha diversity was obtained as the sum of the binary distribution maps of each species. Age is the mean branch length value of all the species found per cell, and number of nodes is the number of nodes from tip to root of all the species found per cell. We made a community matrix from 50 thousand randomly selected points and cells of 2.5 arc min to create a community matrix for the calculation of phylogenetic diversity and weighted endemism. Phylogenetic diversity was calculated with correction for changes in species richness using the tip shuffle model (Laffan & Crisp, 2003), which accounts for changes in species richness when comparing phylogenetic diversity across different communities that could result in skewed results. It randomizes the tips of the tree while preserving the branching patterns and the PD of each community is recalculated. This results in a distribution of phylogenetic diversity values for each community, considering the effects of species richness standardization. Weighted endemism is species richness inversely weighted by species ranges. Both indices were obtained using the package phyloregion R package (Daru, Karunarathne & Schliep, 2020).

Cenocron assemblage’s response to environmental variables

To compare whether these assemblage measures are different between cenocrons, we used a one-way ANOVA with Tukey’s HSD post hoc test, the effect size was labelled according to Makowski et al. (2023). Then, we evaluated whether the assemblage measures are correlated with the environmental variables used to model the species distribution (elevation, slope, topographic roughness index, Thornthwaite’s aridity index, climatic moisture index, annual precipitation, maximum temperature of the coldest month, minimum temperature of the coldest month, maximum temperature of the warmest month, minimum temperature of the warmest month). We applied a Pearson correlation test and compared the coefficient values between cenocrons to detect whether their assemblage responds to environmental variables and the direction of this response (positive vs. negative), and tests with p < 0.05 were accepted as significant. All the analyses were conducted using the R Statistical language (version 4.1.0; R Core Team, 2021).

Diversity in environmental and geographic space

To explore the overlap of species richness in the geographic space, we obtained the sum of the predicted distribution models by cenocron and the total of species. The overlap was also evaluated in environmental space, described as the two principal components that resulted from a principal component analysis to the bioclim layers from WorldClim 2.0 (Fick & Hijmans, 2017) with a resolution of 2.5 arc min, which explains 80% of the environmental variation (56.8%, and 23.2%, respectively, Fig. 1, right). We extracted the values of PC1 and PC2 of the predicted distribution for each species classified by cenocron. Next, we grouped the data by unique coordinates and obtained the richness of each point in environmental space by cenocron and total species. We also counted the number of cenocrons overlapping in each unique coordinate, and selected those coordinates where the three cenocrons were present. The points within the environmental space are interpreted as the environmental transition zone.

Finally, we projected the environmental transition zone into geographic space. This was done by joining the coordinates in the PCA with their original geographical coordinates from the PCA raster layers to visualize the environmental transition zone in a map and compare it with the MTZ (sensu Morrone, 2020, Fig. 1 left).

Cenocron assignment

We assigned species to the cenocrons using simple rules. While previous studies assigned groups to cenocrons based on expert knowledge of the MTZ (Halffter’s original works), used multivariate analyses (Lobo, 2007), or manually linked distribution and phylogenetic data (Roig-Juñent et al., 2018), our approach introduces a more systematic and data-driven method. One of the advantages of this approach is its ability to reduce subjectivity associated with expert knowledge-based classification. Instead, it relies on a set of clear, rule-based criteria that makes it reproducible, as it only uses the species distributions and a biogeographical regionalization.

Our study stands out as the first to test cenocron assemblages for distinct ecological (richness, endemism, and PD) and evolutionary (age and number of nodes) characteristics on a macroecological scale, supporting the idea that cenocrons are not arbitrary groupings but genuine ecological and evolutionary biotic units.

Our method may serve as a preliminary approach for studying cenocrons in the MTZ. While it offers several advantages, the classification method still needs to be improved. Researchers can build upon this methodology by incorporating additional data sources, such as strict phylogenetic information, geological events, or diversification rates.

Conserved niches

Niche conservatism predicts higher species richness within environments occupied by species’ ancestors (Wiens & Donoghue, 2004), and lower age and PD within environments that differ from it (Qian & Ricklefs, 2016). This pattern is observed in the Nearctic and Typical Neotropical cenocrons. The higher species richness within these regions is usually explained by the extended period of colonization, enabling the accumulation of richness over time (Fine, 2015; Pyron et al., 2015; Machac, 2020). Nevertheless, these cenocrons were the last to disperse into the MTZ (Miocene-Pliocene, and Pliocene-Pleistocene, respectively, (O’Dea et al., 2016; Morrone, 2020). Therefore, the species had fewer opportunities to shift niches in comparison with the Mountain Mesoamerican. During this time, as summarized by Mastretta-Yanes et al. (2015), there was a global cooling trend followed by climate fluctuations. In consequence, the distributions of temperate to cold-affinity taxa underwent elevational changes. This is a ‘sky-island’ dynamic, that occurs when species track their preferred habitat in response to environmental changes, which is the main response of a species to environmental changes according to the niche conservatism hypothesis (Lobo, 1999; Stigall, 2014) and a driver of climate-induced vicariance (Claramunt & Cracraft, 2015).

It is important to mention that before the Miocene, South America was an island continent, which resulted in a large diversification and high levels of endemism at continental and bioregional scales (Ricklefs, 2002) and a high specialization to tropical environments. After the Panamanian Bridge lifted and allowed the dispersal of the Typical Neotropical cenocron into the MTZ, birds dispersed as described for the well-understood mammalian interchange (Smith & Klicka, 2010) and their dispersal was also asymmetric. Unlike the northern taxa, the avian families in the Neotropics remain almost entirely restricted to that region (Smith & Klicka, 2010; Smith et al., 2012). For instance, Furnariidae represents the largest continental radiation of Neotropical endemics (Villalobos, Pinto-Ledezma & Diniz-Filho, 2020).

Dispersal from north to south has been more frequent in the MTZ. This is not only because the MTZ had a continuous connection with North America throughout its history (Smith & Klicka, 2010), but also because its mountains are higher when coming from the Neotropics (Janzen, 1967). Therefore, mountain passes are physiological, not topographic, barriers to dispersal (Ghalambor, 2006) and act as such only when there is little or no climatic overlap between the lowlands and the mountains. In temperate zones, the warm season in high elevation can be equivalent to the cold season in the lowlands, which allows species to track the preferred niche seasonally and promotes niche conservatism. This does not happen in the tropics where there is not much seasonal variation in temperature. This hypothesis is linked to the suggestion that seasonal migration promoted colonization of the tropics from the north, which is the mechanism that emberizoid passerines (superfamily Emberizoidea) seem to have followed to colonize the tropics from the north temperate region and is consistent with niche conservatism (Winger, Barker & Ree, 2014). As a result, tropical mountains can be generally structured into distinct thermal zones, with relatively little thermal overlap between low and high-elevation sites across seasons (Muñoz & Bodensteiner, 2019).

Convergent niches

Having discussed the evidence for niche conservatism in the MTZ, we now turn our attention to the Mountain Mesoamerican cenocron, which exhibits patterns consistent with niche convergence. The Mountain Mesoamerican cenocron does not follow the pattern expected when niches are conserved, since it does not correlate with the environmental variables. Its species are distributed at the lower and more variable temperature extremes of the Typical Neotropical cenocron’s distribution (quadrant IV in the environmental space), suggesting that the Mountain Mesoamerican cenocron may represent a subset of the Typical Neotropical cenocron in multivariate space.

The Mountain Mesoamerican cenocron, representing the oldest dispersal in our analysis, has been subject to many climatic fluctuations and landscape shifts. Its species dispersed from South America to Central America during the Oligocene and diversified in the mountains. They further extended their range into Mexico during the Pliocene (Rocha-Méndez et al., 2019; Morrone, 2020). This dispersal took place during Earth’s most recent period of sustained global warmth (Dekens, Ravelo & McCarthy, 2007), characterized by warmer temperatures at high latitudes and reduced temperature differences between the Equator and the poles (Fedorov et al., 2013). At that time, mountains where the Mountain Mesoamerican cenocron dispersed were covered by warm-tempered forests (Mastretta-Yanes et al., 2015). Additionally, there was a highland corridor connecting the mountains of Central America and Mexico through the Isthmus of Tehuantepec (Mastretta-Yanes et al., 2015), which played a role in the geographical setting of the Mountain Mesoamerican cenocron’s dispersal (Halffter, Favila & Arellano, 1995; Halffter & Morrone, 2017; Rocha-Méndez et al., 2019; Morrone, 2020). The Mountain Mesoamerican cenocron spread northwards from Panama across the Central American lowlands, and its lineages started to diverge during the Pliocene at 5.8 Ma (Rocha-Méndez et al., 2019). Subsequently, a period when intense ice age cycles began (Dekens, Ravelo & McCarthy, 2007). In response to these cycles, according to the niche conservatism hypothesis, the main response is dispersal (Lobo, 1999; Donoghue, 2008). However, when it is not possible to disperse back, adaptation to novel climatic niches promotes speciation and/or impedes extinction (Cooney, Seddon & Tobias, 2016).

The fragmentation of the once continuous cloud forest due to Pleistocene climatic fluctuations (Luna Vega et al., 1999), the partial loss of the connection along the Isthmus of Tehuantepec (Mastretta-Yanes et al., 2015), and the lift of high stratovolcanoes in the Trans-Mexican Volcanic Belt (Ferrari et al., 2012) promoted semi-permeable barriers to dispersal (Rocha-Méndez et al., 2019). These events fragmented the cloud forest into areas of high endemism (Luna Vega et al., 1999) and great genetic diversity (Rocha-Méndez et al., 2019). A vicariance model of the Mountain Mesoamerican forests and sky island dynamics explain biogeographic patterns (Luna Vega et al., 1999; Sánchez-González, Morrone & Navarro-Sigüenza, 2008; Fjeldså, Bowie & Rahbek, 2012; Sosa & Loera, 2017) and intraspecific genetic structure (Rocha-Méndez et al., 2019). Furthermore, the current diversity gradients are related to climate (Sosa & Loera, 2017). This agrees with Peterson, Soberón & Sánchez-Cordero (1999), who noticed the importance of separating geographical and environmental space and adding a temporal scale when evaluating niche conservatism. Overall, in situ endemic radiations associated with mountain uplifts and climatic fluctuations (Hagen et al., 2021) may explain the vicariance pattern in the Mountain Mesoamerican cenocron, along with our results that are consistent with a niche convergence diversification model.

Environmental Transition Zone

The Mexican Transition Zone’s provinces overlap with temperate and tropical conditions (Fig. 1), which, according to the environmental overlap hypothesis (Janzen, 1967), make these mountains a semi-permeable barrier. This agrees with the definition of a transition zone as partial barriers or filters that restrict differentially the distribution of each biotic component (Ferro & Morrone, 2014). The area found in our results corresponds to the MTZ (Morrone, 2020), but goes beyond the Nicaraguan Depression. This region represents the distributional boundary for many bird taxa (Sánchez-Ramos et al., 2018), and is considered the southern limit of the MTZ (Halffter & Morrone, 2017). The simpler explanation for this could be the greater dispersal abilities of birds. Nevertheless, the environmental overlap model suggests that the MTZ has a larger overlap with the Neotropics, making it easier for populations to disperse upwards in the mountains or northwards.

Implications for conservation

Conservation efforts should aim to create and preserve corridors that allow species to move and adapt to changing conditions. The environmental transition zone we have described is the place where species from various evolutionary backgrounds coexist, face changing climates by tracking their preferred niches, and encounter isolated habitats that may foster speciation. This underscores the vital importance of preserving this region.

The MTZ is not just a biodiversity hotspot; its diverse species and ecological niches might offer resilience against climate change. Furthermore, the MTZ serves as a valuable natural laboratory for studying ecological and evolutionary processes (Morrone, 2020). Research conducted in the MTZ has the potential to enhance global biodiversity understanding and guide conservation strategies especially when considering the historical context. For instance, Mexican conservation strategies should consider the protection of Neotropical, Nearctic, and species originated in situ in the MTZ.

Study limitations

While our study provides new insights on how to study MTZ’s diversity patterns, we must recognize its limitations. The cenocron concept, which we consider crucial for studying the MTZ, might not be well-known among the broader biogeography and macroecology community. Our study introduces categorization criteria for cenocrons, to make them more accessible and understandable. Yet, there is still potential for improving the methods keeping in mind the importance of incorporating both evolutionary relationships (phylogenetic factors) and environmental influences (ecological factors). Alternative ways to do so are spatial phylogenetics (Mishler, 2023) or the comparative phylogenetic method applied to niche evolution (Evans et al., 2009).

A special challenge in our study was choosing Species Distribution Models (SDM) as data input. since the modeling of species distributions is a complex task with no simple solutions. There are many variations in the way of constructing a model (e.g., data cleaning, calibration area choosing, variable selection or hyperparameter tunning), common problems (like overparameterization), and multiple ways to evaluate the result. Fortunately, there are automatized methods, like the one we used, that allow to perform data-driven variable selection and enable the automation of model hyperparameter tuning. In essence, SDMs are valuable tools for overcoming the Wallacean shortfall, which refers to the incomplete knowledge of species distributions (Bini et al., 2006). SDMs provide the means to compare species with varying data quality and quantity, thus preventing the overrepresentation of species boasting thousands of records while accommodating those with more limited data.

Our primary objective was to understand the broader impact of the environment on diversity within the MTZ. Although not trivial, the methods to model each species will not affect noticeably the results since gradients like the ones described herein have been identified using many data types (points, polygons, lines, grids, transects). Yet, the use of SDM aims to enhance the integration of analyses with spatial methods, and the increased resolution facilitates the identification of gradients in environmentally heterogeneous areas, like the MTZ.