Stairway to extinction? Influence of anthropogenic climate change on distribution patterns of montane Strigiformes in Mesoamerica

. Although anthropogenic climate change (ACC) is a global phenomenon affecting all ecosystems, its effects are especially relevant in certain ecosystems, such as tropical montane forests. Responses of montane species to ACC in Mesoamerica remain unclear, limiting our ability to assess their vulnerability and the impacts on these ecosystems overall. To understand mechanisms underlying the distribution and vulnerability of montane faunas, we analyzed the influence of ACC on the geographic distribution of owls (order Strigiformes), which are a group of top avian predators distributed in montane forests. Using ecological niche models, we estimated the potential distributions of 35 species at present and under projected future climates (2050 and 2070) and analyzed changes in distributional patterns in terms of range size and altitudinal distribution for each species, as well as spatio-temporal patterns of species richness. Most of our projections (~86%) were consistent with the widely accepted hypothesis of species range shift to higher altitudes combined with reduction in distribution, as species try to track their climatic preferences. Interestingly, the mid-elevation species emerge as the most strongly affected by ACC, showing the highest rates of change. All climate scenarios produced a similar pattern of change in owl species richness, but they differed in the total number of species, a loss of 11 species and a maximum gain of seven species. Species richness remained relatively constant at mid elevations, whereas the greatest losses were in the highlands and the contiguous lowlands. Overall, our results suggest a severe impact of ACC in the coming decades for most owls of Mesoamerican montane forests.


INTRODUCTION
Anthropogenic climate change (ACC) is one of the most challenging threats to biodiversity in the 21st century, leading to an increase in temperature, changes in precipitation regimes, and more frequent extreme weather events because of human activities (Beaumont et al. 2011, IPCC 2014. Although ACC is a global phenomenon, its effects are not homogeneous; it has especially strong effects on ecosystems located at the extremes of climatic gradients, such as tropical montane forests (Pacifici et al. 2015, Brotons et al. 2019. Tropical montane forests are ecosystems that extend from the base of mountains to the tree line, and they are highly vulnerable to ACC because of their  (Dinerstein et al. 2017). ecological, topographical, and historical characteristics (Price et al. 2011, Freeman et al. 2018, Bender et al. 2019, Rahbek et al. 2019a. In these regions, the combination of temperature and topography generates a marked climatic zonation with differing temperatures across short distances, favoring heterogeneous habitat configuration as altitude increases (Rzedowski 2006).
In response to climatic alterations, species in montane regions are expected to search for suitable sites by moving to higher altitudes, depending on topography, ecological specialization, and their climatic preferences (Rojas-Soto et al. 2012, Şekercioĝlu et al. 2012, Freeman et al. 2018, Bender et al. 2019. Thus, the first expected change is a modification of the distribution area, which may also produce alterations in demography, synchrony of events (e.g., flowering, migration, or breeding), biotic interactions (e.g., prey availability), and community dynamics (Pacifici et al. 2015, Brotons et al. 2019. Species that are distributed in environments that limit their movements, such as those at higher elevations, could be more strongly impacted and may even face local extinction (Lenoir andSvenning 2015, Foden et al. 2019). However, beyond these general theoretical expectations, details about the response of montane species to the accelerated climatic alterations are still poorly known. Therefore, we need to study species' distribution dynamics to improve our understanding of potential shifts in species and communities (e.g., Prieto-Torres et al. 2021).
The Mesoamerican montane forests (MMF), comprising the highlands of central Mexico to Panama (Fig. 1), are considered one of the most diverse and relevant regions in Mesoamerica. They lie in a transition zone between the Neotropical and the Nearctic realms, representing an area of biotic interaction between faunas with different affinities. Nearctic lineages tend to occupy the upper mountains, Neotropical lineages are distributed in the lower parts, and the intermediate altitude zone shows mixed biotas . Geographically, species richness reaches its maximum in the contact zones between lowlands and mountains along the eastern Mexico slope and decreases as altitude increases (Koleff et al. 2008, Navarro-Sigüenza et al. 2014. Given this high diversity, the MMF are recognized as a center of diversification, endemism, and biogeographic transition for many taxa (Navarro-Sigüenza et al. 2007, Sánchez-Ramos et al. 2018, Moreno-Contreras et al. 2020, Morrone 2020, Ramírez-Albores et al. 2020. Temperature and precipitation are the main factors controlling the range limits of montane species. Recent studies show spatially heterogeneous temperature and precipitation patterns in Mesoamerica over the last three decades, with a slight increase in mean annual temperature and more pronounced dry and wet seasons (Cuervo-Robayo et al. 2020, Wootton et al. 2022. Combined with the complex topography, these are key factors that contribute to a narrow range of suitable conditions for montane species (Price et al. 2011, Payne et al. 2017, Freeman et al. 2018). In addition, the orientation of the main mountain ranges is crucial in determining precipitation and humidity patterns (Rzedowski 2006, Challenger and (Enríquez 2017).
Although owls continue to be one of the least studied bird taxa in Mesoamerica , recent studies show a growing interest in exploring ecological traits, population trends, and geographic distributions (e.g., Valencia-Herverth et al. 2012, Vázquez-Pérez and Enríquez 2016, Enríquez 2017, Fernández Martínez 2017, Ramírez-Santos et al. 2018, Ayma et al. 2019. One approach to exploring the potential change in the distribution of owls is the use of ecological niche models (ENMs). These models are mathematical abstractions that estimate a species' distribution based on correlative relationships between species' occurrences and environmental factors (Peterson et al. 2011, Zurell andEngler 2019). When models are projected onto future climates, they can suggest geographic distributional shifts driven by alternative ACC scenarios (Foden et al. 2019). This information can then be used to test hypotheses about the dynamics of the species distribution under climate variations.
Here, our main goal was to analyze the potential effect of anthropogenic climate change on the geographic distribution of Strigiformes inhabiting the MMF. To achieve this, we proposed to (1) characterize current distributional patterns in terms of range size and altitudinal distribution of the species, (2) determine expected distributional changes of the species under alternative ACC scenarios for the years 2050 and 2070, and (3) examine how patterns of owl species richness might change with ACC. We hypothesized a shift of most species toward higher altitudes coupled with a reduction in their range size, expecting the most severe changes in species located at altitudes above 1500 m. For species below 1500 m, the direction of their movements is uncertain, dependent on local conditions and individual species' climatic preferences, but if they move to lower altitudes the range size will increase. Species with a wide altitudinal distribution are not expected to show significant change because they are able to tolerate a broad range of climatic conditions. Finally, we expected an altitudinal shift with species moving from lower elevations to the highlands, resulting in lower owl species richness in the lowlands and a similar number of species at mid-elevations. Lower owl species richness was also expected in the highlands because of a reduction in the range size of the species present there.

Study area
Montane forests are herein defined as the areas extending from the base of the mountain to its tree line and comprise oak forests, pine forests, mixed coniferous forests, cloud forests, and other less extensive vegetation types, such as spruce-fir forests (Rzedowski 2006). These areas exhibit a humid subtropical to temperate climate with annual temperatures of 5 °C to 25 °C and annual precipitation of 600 mm to 1200 mm, although some areas record over 3000 mm. Overall, the altitudinal range of the considered vegetation types vary from 800 m to 3600 m (Rzedowski 2006, Challenger and. The Mesoamerican montane forests were delimited according to the biogeographical regionalization of the Neotropics (Morrone 2014) and cropped with the terrestrial ecoregions classification (Dinerstein et al. 2017) including the aforementioned vegetation types. The study area is found from north-central Mexico to the northern tip of Panama and grouped into six biogeographic provinces with discontinuities in the Isthmus of Tehuantepec and the Nicaraguan Depression ( Fig. 1). We kept all historical records to calibrate the models, because the records for most species were from between 1900 and 1965, and the climate in northern Mesoamerica did not vary drastically throughout the 20th century and the first decade of 21st century (Cuervo-Robayo et al. 2020). Third, to reduce possible sampling biases and model overfitting, we retained only the records corresponding to localities separated by at least the mean distance between occurrence records (Appendix 1) with the "spThin" package (Aiello-Lammens et al. 2015) in R (R Core Team 2021). We selected 35 of the 38 Mesoamerican owl species (scientific and common names provided in Appendix 1), considering (1) distribution, functional role, and altitudinal patterns in the MMF (Valencia-Herverth et al. 2012, Vázquez-Pérez and Enríquez 2016, Enríquez 2017, Fernández Martínez 2017, Billerman et al. 2022and (2) presence records in at least 20 localities in the study area (the number of localities proposed to be sufficient to create and validate the model; Pearson et al. 2007). Although some species can be more strongly associated with other vegetation types (e.g., Burrowing Owl, Athene cunicularia), all were included because there were reliable records of their presence within the delimited area of analysis. Three species (Glaucidium gnoma, G. hoskinsii, and Megascops lambi) were excluded because they did not meet the criterion of a minimum of 20 records.

Species sampling, locality records, and climatic variables
A subset of potential predictor variables was selected (described below) from the bioclimatic variables provided by WorldClim version 2.1 (Fick and Hijmans 2017) at a resolution of 2.5′ (∼5 km²), which are derived from monthly temperature and rainfall values. From the same repository, an elevation layer with the same resolution was downloaded for spatial analysis of the distributions. The model calibration area was estimated for each species using the intersection of occurrence records with the terrestrial ecoregions (Dinerstein et al. 2017) and the biogeographical provinces of the Neotropics (Morrone 2014), considering that both approaches can address their historical and ecological limits (Soberón and Peterson 2005, Prieto-Torres et al. The climate projections used were the RCP 4.5 and RCP 8.5 scenarios for the years 2050 and 2070, based on CMIP5 data available at the WorldClim portal. RCP 4.5 is an intermediate stabilization scenario in which emissions peak around 2040 and then decline, whereas RCP 8.5 represents very high emissions of greenhouse gasses and few climate change mitigation policies (IPCC 2014). Thus, for each species we considered the present plus four future scenarios. Given that the choice of different general circulation models (GCMs) has been identified as a source of variability between distribution models (Zappa andShepherd 2017, Fajardo et al. 2020), we employed the GCM-CompareR platform (Fajardo et al. 2020), which adopts a "storyline" approach to classify GCMs into narratives representing future climate conditions (Zappa and Shepherd 2017). Projections can also be refined by identifying how the future estimates might depend on initial experimental conditions, algorithms, and model biases (Guevara et al. 2018). MIROC5 was the selected model because it incorporates a realistic simulation of El Niño-Southern Oscillation, on the basis of an improved estimation of precipitation, equatorial ocean surface temperature, and zonal mean atmospheric fields (Watanabe et al. 2010). Our analyses were performed under the expectation of higher temperatures and slightly less precipitation than present for both 2050 and 2070.

Ecological niche modeling
Ecological niche models were constructed separately for each species by using Maxent version 3.4.3 (Phillips et al. 2017) with the "kuenm" package (Cobos et al. 2019). Maxent has shown good performance with presence-only data, whereas kuenm allows to generate candidate models with different parameterization (regularization factors and features). To minimize collinearity of predictor variables and model overfitting, a subset of potential predictor variables was constructed on the basis of a Pearson correlation coefficient (r < 0.8) and the Variance inflation factor (VIF < 10). The parameters that optimize the balance between goodness-of-fit and complexity (Appendix 2) were identified by testing a variety of regularization multipliers (0.1 to 1 every 0.1, and 1 to 4 at intervals of 1) and several combinations of five feature types ("basic" combination of linear, quadratic, product, threshold, and hinge; sensu Phillips et al. 2017). The candidate models' performance was evaluated via the partial area under the ROC curve (partial ROC), an omission rate (E) lower than 5%, and least complexity according to the corrected Akaike's information criterion (AICc), selecting those that met all three criteria (Appendix 2). The models were constructed using a random sample of 70% of the locality records as training data and the remaining 30% as validation data, 50,000 background points and 5 replicates.
Future climate projections might include non-analogous conditions (i.e., conditions beyond those available in the model calibration), which can lead to uncertainty in geographical predictions (Owens et al. 2013, Guevara et al. 2018. To identify areas where extrapolation risks could be expected, constructed models were transferred under two assumptions (unconstrained extrapolation and clamping), and for each model the response curves and the geographical prediction were evaluated. Also, a mobility-oriented parity test (Owens et al. 2013) was performed to identify sites with a high degree of environmental dissimilarity. Finally, to generate distribution maps, the median values of the replicates were considered. Presence-absence maps were obtained by converting the cloglog output format (a continuous scale ranging from 0 to 1 of environmental suitability values) to binary, using the 10-percentile training presence threshold values (Appendix 2) with the "raster" package (Hijmans 2021). This threshold omits the 10% of records with the lowest suitability values (under the assumption that these sites are not representative of the species' requirements), to exclude outliers (Escalante et al. 2013) and thus minimize commission errors (i.e., predicting a species as present when it is not).

Spatial analysis of distributions
Potential impacts of ACC on the species were based on changes in three aspects of their distributions: extent of geographic distribution, altitudinal shift, and species richness. These changes were calculated by subtracting future from current potential distributions, and all post-modeling calculations were performed with the "raster" package (Hijmans 2021). The extent of geographic distribution was characterized by quantifying the range size (km²) and presence within the study area (percentage). The range size was classified into three categories according to the total area (km²) in which the species was predicted to be present in the current scenario (Appendix 3): large (upper quartile > 180,000 km²), intermediate (< 180,000 km², > 20,000 km²), and small (lower quartile < 20,000 km²). The altitudinal distribution was classified into four categories on the basis of the quartile ordination of the current species distribution (Appendix 3; Prieto-Torres et al. 2021) and the median value of each species: highlands (upper quartile 2100 m), mid-elevations (< 2100 m, > 1100 m), lowlands (lower quartile 1100 m), and generalists (difference between Q3 and Q1 > 1000 m). Also, a Kruskal-Wallis test was implemented to estimate differences among altitudinal and range size categories on the basis of species' range changes. Last, current and future patterns of owl species richness at the regional scale were estimated by adding all the binary distribution maps obtained and converting them to a standardized raster (0-1 values). The patterns obtained for 2050 and 2070 under the two RCP scenarios were compared by subtracting the owl species richness per site with respect to the current distribution. All analyses were performed in R.

Species models and current distributional patterns
The selected species belonged to 12 genera. Seven genera are represented by a single species, whereas Glaucidium and Megascops are represented by seven and 10 species, respectively (Appendix 1). The models exhibited good performance (i.e., far more accurate than expected by chance), with significant values for the partial ROC test (1.08-1.95, p < 0.05), and low omission rates (< 10%); model parameters and performance metrics are detailed in Appendix 1, and individual responses for each modeled species are in Appendix 4.
The area of the predicted distribution varied considerably among species (Appendix 3), from 3000 km² for the Costa Rican Pygmy-Owl (Glaucidium costaricanum) to 375,000 km² for the American Barn Owl (Tyto alba). According to our range size categories, nine  (Fig. 2a,  Appendix 3). Categorization by altitudinal range resulted in 13 lowland, 13 mid-elevation, two highland, and seven generalist species (Appendix 3). The altitudinal distribution showed species throughout the altitudinal gradient (Fig. 2b). The Stygian Owl (Asio stygius) had the broadest altitudinal range, from sea level to 3200 m, and the narrowest altitudinal range was shown by the Striped Owl (Asio clamator) from sea level to 1250 m. In addition, even though some species had broad altitudinal ranges, their distributions were generally skewed toward one end of the altitudinal gradient. For example, although the Northern Saw- whet Owl (Aegolius acadicus) had its lower altitudinal limit around 1000 m, it was mostly distributed in the upper mountains, and although the Crested Owl (Lophostrix cristata) was mostly predicted at elevations below 250 m, some areas of its potential distribution reached up to 2500 m. Sixty percent (n = 9) of the species found in highlands and those found in mid-elevations were highly restricted to their respective altitudinal ranges, with more than 75% of their potential distribution occurring in that altitudinal range (Fig. 3)

Impacts of ACC
Areas with non-analogous climates were identified by the mobility-oriented parity (MOP) test. However, they represented a low proportion of our predictions (< 10% on average) and thus were included in subsequent analysis, although we treated them with greater caution (especially in RCP 8.5). These areas were located mainly in the lowlands near montane areas, such as the northern lowlands of the Pacific and Gulf of Mexico slopes, the Nicaraguan Depression, and the lowlands of southern Central America. Considering the MOP results, the response curves, and the fit of geographic predictions, a single extrapolation method, extrapolation by clamping, was selected for all scenarios (Appendix 2).
When considering the proportion of the species' total distribution range that falls within the MMF as defined at present (Fig. 3), 74% (26) of the owl species showed a higher proportion of their potential distribution within the MMF, likely because of their altitudinal shifts. Especially, the proportion of the distribution of the Tamaulipas

Species richness patterns
The estimated species richness varied from five to 21 species under the current climate scenario. Sites of high species richness (> 10 species, dark-colored sites in Fig. 3) were mostly located in the contact zones between highlands and lowlands as well as between biogeographic provinces: for example, between the Sierra Madre Oriental and adjacent lowlands, between the Transmexican Volcanic Belt with the Sierra Madre Occidental, and in the southern part of the Sierra Madre del Sur. The lowest species richness was observed in northern Mexico, especially in the northwest of the Sierra Madre Occidental and in the western part of the Nuclear Central American Highlands (Fig. 3). Among biogeographic provinces, the Transmexican Volcanic Belt had the highest species richness with 23, followed by the Sierra Madre del Sur with 20; the Sierra Madre Occidental had the lowest richness, with 17.
All climate scenarios produced a similar pattern of change in owl species richness, but they differed in the total number of species (Fig. 5). The main species losses occurred in the northern part of the Sierra Madre Occidental, the central part of the Sierra Madre Oriental, the southern portion of the Transmexican Volcanic Belt and Sierra Madre del Sur, and the northwestern part of the Nuclear Central American Highlands. The RCP 8.5 for 2050 and 2070 showed particularly widespread losses in these areas. The areas with species gain were in the southern Sierra Madre Occidental, northern Sierra Madre Oriental, and western Transmexican Volcanic Belt. Overall, the projected patterns showed a loss of 11 species and a gain of six, except for RCP 4.5 in 2050, where losses of 10 species were estimated, and for RCP 8.5 in 2070, which showed a maximum gain of seven species.

DISCUSSION
According to the models generated, ACC is predicted to have severe impacts in the coming decades for most of the owls of the Mesoamerican montane forests. Our results are consistent with the widely accepted hypothesis of species range shift to higher altitudes combined with reduction in distribution as species try to track their climatic preferences (Table 1; Parmesan 2006, Pacifici et al. 2015, Bender et al. 2019. We found this pattern in most of the owl species (~86%), and although the scenarios in 2070 show the most severe results, the projections for 2050 require more attention. The mid-elevation species emerge as one of the groups that are most vulnerable to ACC, showing the highest rates of change (Table 1, Fig. 4; Lenoir and Svenning 2015). This is the case with the Tamaulipas Pygmy-Owl and Eastern Screech-Owl, which show reductions of up to 73% in range size and altitudinal shifts of around 650 m. This pattern may be explained by the complex climatic patterns and notable changes in conditions over relatively short distances in mid-altitude areas (Rahbek et al. 2019a(Rahbek et al. , 2019b. Moreover, because highland species cannot expand their distributions beyond the mountain habitable zone (Şekercioglu et al. 2012(Şekercioglu et al. , Freeman et al. 2018, these species are also projected to be very strongly affected, potentially threatening their persistence in the region (Şekercioglu et al. 2007(Şekercioglu et al. , Urban 2015. MMF have undergone significant transformation because of land use changes, and in many cases have undergone accelerated biodiversity loss over the last 50 years (Challenger andSoberón 2008, Enríquez 2017). The synergistic effect of land use changes and ACC could limit species' ability to follow their climatic preferences (Rojas-Soto et al. 2012) and result in even greater negative impacts (Jetz et al. 2007, Beaumont et al. 2011).
As hypothesized, species richness remained relatively constant at middle elevations and the largest losses were in the highlands and the contiguous lowlands ( Fig. 4; Bender et al. 2019). When a species shifts upward, it can be replaced by species from lower elevations (Bender et al. 2019), so altitudinal shifts do not necessarily lead to a net decrease in species richness at midelevations (Colwell et al. 2008, Rahbek et al. 2019a). However, species richness is expected to decrease at both ends of the altitudinal gradient, though for different reasons at each end. The lowlands surrounding the MMF of northern and central Mexico might show a reduction because there are no longer species to replace the emigrants in those areas (Dunn and Møller 2019). Meanwhile, highland species (e.g., Aegolius spp., Saw-Whet Owls) cannot expand their distributions beyond the mountain habitable zone, so their distribution areas are expected to decrease (Şekercioglu et al. 2012(Şekercioglu et al. , Freeman et al. 2018. Moreover, the gains and losses in owl species richness suggest that the highlands are climatically different from the lowlands (Fig 4; Rahbek et al. 2019a), favoring geographic isolation of remaining sites with high species richness (Şekercioĝlu et al. 2012(Şekercioĝlu et al. , Payne et al. 2017. Special attention should be paid to the Transmexican Volcanic Belt and Sierra Madre del Sur, where the greatest reduction in owl species richness is expected to occur. Although our results are interesting and informative, they should be interpreted with caution. Species distribution is shaped by factors other than climate that our models are unable to process but should be considered (Peterson et al. 2011, Guisan et al. 2017). These include ecological factors like density-dependent interactions (e.g., prey availability) and population dynamics (Parmesan 2006, Dunn andMøller 2019), or historical processes such as speciation or colonization (Espinosa and Ocegueda 2008, Rahbek et al. 2019b). Moreover, multiple compensatory mechanisms such as modifying activity patterns or physiological adjustment have been recorded when climatic changes exceed the natural variation of the region (Prieto-Torres et al. 2021). Mechanisms such as diet modification, reduction of body mass, demographic responses, or use of secondary vegetation may buffer the effects of ACC (Newton 2003, Dunn andMøller 2019), although ecological information is scarce for most Mesoamerican owls. Characterization of species' responses to environmental changes is clearly a complex task, but information on their distributions is a critical starting point to understand underlying mechanisms (Foden et al. 2019).
Species with small niches are typically specialists, both ecologically and in terms of distribution, so in theory they are less likely to be able to adapt to new climatic conditions (Wiens et al. 2010). Our findings are in line with this notion because the species with the largest changes were specialists, restricted to a particular altitudinal range or area. Moreover, the projection of mostly loser (specialists) and few winner species (generalists) under new climatic conditions supports the hypothesis that the Avian Conservation and Ecology 17(2): 37 http://www.ace-eco.org/vol17/iss2/art37/ ecological specialization of a species is a key attribute influencing species vulnerability (Thuiller et al. 2005, Pacifici et al. 2015, Guisan et al. 2017). The pattern found here (range shift to higher altitudes and reduction in range size) supports previous studies conducted on several study systems. For example, similar trends were found within a single province, the Sierra Madre Oriental (Rojas-Soto et al. 2012); indeed, just like our study, the central part of the Sierra Madre Oriental was found to be the most strongly affected. These results have also been found when assessing the effects of projected climate changes on endemic bird species in the main montane regions of Mexico (Sierra-Morales et al. 2021) and tropical mountain regions worldwide (Freeman et al. 2018).
Consistent with previous studies (Sánchez-Ramos et al. 2018, Bender et al. 2019), we found that species that inhabit higher altitudes had smaller range sizes. This was evident in owls with median elevations above 2000 m, whose distributions were smaller than 50,000 km² and for which ≥ 90% of the range was contained within the MMF. The distribution limits of a species are mostly Avian Conservation and Ecology 17(2): 37 http://www.ace-eco.org/vol17/iss2/art37/ influenced by large-scale environmental variables, such as temperature, precipitation, humidity, or vegetation type (Sexton et al. 2009, Freeman et al. 2018, Rahbek et al. 2019b Urban 2015). Nine of the owl species we analyzed (26%) have a range size that could justify including them in some category of risk, although they are not currently recognized (IUCN 2021, https://www.iucnredlist.org). For example, two species (Glaucidium costaricanum and Megascops clarkii) had an estimated distribution of less than 5000 km², which could qualify them as Endangered, whereas seven could be Vulnerable given their distributions smaller than 20,000 km² (BirdLife International 2022). From this conservation perspective, further attention is also required to explore changes in species turnover (e.g., Ochoa-Ochoa et al. 2014) and changes at the level of individual species, communities, or ecosystem functioning (Pacifici et al. 2015, Brotons et al. 2019). These potential effects may provide a picture of the possible impact of human-induced changes not only on charismatic owls, but on the already threatened whole biota of our planet. .753 † Feature classes are here abbreviated as follows: l= linear, q = quadratic, p = product, t = threshold, and h = hinge. † Binarization threshold refers to the minimum environmental suitability value used to convert continuous model output into a presence-absence estimation (Peterson et al. 2011). Given each species´ climatic preferences, the suitability values change and thus, the choice of the threshold.