Global climate change is associated with changes in precipitation patterns and an increase in extreme weather events, which might shift the geographic distribution of species. Despite the importance of this topic, information is lacking for many species, particularly tropical birds. Here, we developed species distribution models (SDMs) to evaluate future projections of the distribution of the widespread Buff-bellied Hummingbird (Amazilia yucatanensis) and for each of the recognized subspecies (A. y. yucatanensis, A. y. cerviniventris, A. y. chalconota), under climate change scenarios. Using SDMs we evaluate current and future projections of their potential distribution for four Representative Concentration Pathway (RCPs) for the years 2050 and 2070. We also calculated the subspecies climatic niche breadth to test the relationship between their area of distribution and climatic niche breadth and their niche overlap. Future climate-change models suggested a small increase in the potential distribution of the species and the subspecies A. y. yucatanensis, but the predicted potential geographic range decreased in A. y. chalconota and remained unaffected in A. y. cerviniventris. The climatic niche of A. y. cerviniventris contained part niche space of A. y. yucatanensis and part of A. y. chalconota, but the climatic niches of A. y. yucatanensis and A. y. chalconota did not overlap. Our study highlights the importance of correctly choosing the taxonomic unit to be analyzed because subspecies will respond in a different manner to future climate change; therefore, conservation actions must consider intrinsic requirements of subspecies and the environmental drivers that shape their distributions.
Introduction
Anthropogenic climate change is a complex phenomenon that affects natural ecosystems and a biodiversity threat in several regions in the world (Freeman et al., 2019; Townsend Peterson et al., 2001; Thomas et al., 2004). Increase in temperature, shifts in precipitation patterns, and the increasing occurrence of extreme weather events jeopardize the persistence of species and their ecosystems due to its potential to affect areas far from human settlements (IPCC, 2014; Malcolm et al., 2006). Current climate change is already altering the geographic distribution of plant and animal species and several others are expected to increase/decrease and/or to shift their geographic distribution in response to climate change (Borzée et al., 2019; Freeman et al., 2019; Mendoza-González et al., 2013; Ornelas et al., 2018). Studies carried out with distribution data (1967 to 2002) of 56 bird species showed that birds with northern distributions across the temperate northern hemisphere are not expanding their ranges southward due to global warming (Hitch & Leberg, 2007; Thomas & Lennon, 1999). However, birds with a southern distribution showed a 2.35-km/year northward expansion in their northern limit (Hitch & Leberg, 2007; Thomas & Lennon, 1999). Despite the high vulnerability of tropical bird species to global climate change, studies show that these have received less attention than birds from temperate regions (Freeman et al., 2019; Harris et al., 2011). In the tropics, the general trend shows that species may be displaced to higher elevations, decreasing its distribution area as the climate warms in the lowlands (Chen et al., 2011; Harris et al., 2011; Peh, 2007; Prieto-Torres et al., 2020).
Species Distribution Models (SDMs) are spatial representations of species presence probability or abundance and are excellent tools to understand the factors that affect the potential distribution of organisms at different scales (Borzée et al., 2019; Pearson & Dawson, 2003). SDMs are generated using responses to environmental predictor variables, spatial relationships such as convex hulls or interpolation, or a combination of both. These are projected across landscapes, and therefore, might be helpful when acquiring occurrence data on the distribution of species is not achievable, or when the ecological variables related to the distribution of the species had changed (Borzée et al., 2019; Freeman et al., 2019; Guisan & Thuiller, 2005). SDMs have been widely used to predict the potential distribution as a function of climate change in several organisms (Borzée et al., 2019; Kafash et al., 2018; Ornelas et al., 2018; Ramírez-Preciado et al., 2019), including birds, in which SDMs were effective to predict their current and future distribution (Atauchi et al., 2020; Freeman et al., 2019; Townsend Peterson et al., 2001; Prieto-Torres et al., 2020; Şekercioğlu et al., 2012).
Several studies have evaluated the effects to climate change on current distribution of bird species and showed that the potential distributions in some species is reduced (Atauchi et al., 2020; Knudsen et al., 2011), expand northwards (Hitch & Leberg, 2007; Lara et al., 2012), shifted and contracted when exacerbated by habitat loss (Colyn et al., 2020), and in some cases the impacts of climate change on current distributions of birds were minimal, particularly in range-restricted species (Freeman et al., 2019; Lara et al., 2012). Hummingbirds, as pollinators, are often immersed in highly dependent interactions with plants (Correa-Lima et al., 2019; Ornelas et al., 2004), and constitute the most diverse ensemble of specialized nectarivorous birds in the Americas (Sonne et al., 2016; Zanata et al., 2017). Due to this characteristic, any environmental factor that affects the distribution of hummingbirds would also affect the distribution of their plants, and vice versa (Buermann et al., 2011). For instance, synchronization of the flowering period with the movements and foraging activity of hummingbirds is important for the maintenance of ecological networks, given that plants provide nectar resources for hummingbirds (Junker et al., 2013; López-Segoviano et al., 2018; McKinney et al., 2012). In comparison to range-restricted species, widely distributed hummingbird species facing broad environmental conditions that extend their geographical range might migrate tracking new floral nectar resources and colonizing new adequate habitat under future scenarios more readily than range-restricted species (Licona-Vera & Ornelas, 2017; Ornelas et al., 2015).
Here, we assessed future potential distribution of the widespread Buff-bellied Hummingbird (Amazilia yucatanensis) under climate change scenarios. Because predictions below the species level can lead to very different scenarios than those expected when evaluated at the species level, particularly for widespread taxa (Gonzalez et al., 2011; Mota-Vargas & Rojas-Soto, 2016; Smith et al., 2019), the effects of climate change were also evaluated at the subspecies level. Specifically, the objectives of this study were to: (1) estimate the current potential distribution of Amazilia yucatanensis as a full species and to predict the potential impact of climate change on its geographic distribution using different climate change scenarios for the years 2050 and 2070; (2) estimate the current potential distribution of each of the subspecies individually (A. y. yucatanensis, A. y. cerviniventris and A. y. chalconota), and predict the potential impact of climate change on their geographic distribution using different climate change scenarios for the years 2050 and 2070, as to know the prediction differences between the different models; (3) test whether ecological niches have diverged significantly among the three distributional areas of these subspecies; and (4) to assess whether the breadth of climatic niche related to the predicted future changes in the geographical distribution of subspecies. We hypothesized that due to morphological differences between the subspecies and the environmental variation across their geographical distribution, the distribution ranges in future scenarios will differ when evaluated at the subspecies level than when evaluated at the species level and therefore their climatic niche will be different.
Methods
Species Target
The Buff-bellied Hummingbird (Amazilia yucatanensis) is distributed on the slope of the Gulf of Mexico from Tamaulipas and S Texas to the Yucatán Peninsula, NW Guatemala and N Belize (Howell & Webb, 1995; Schuchmann, 1999). Despite the complicated nomenclatural and taxonomic history of Amazilia and allies (Stiles, Piacentini, et al., 2017; Stiles, Remsen, et al., 2017), the subclade composed of Amazilia rutila, A. yucatanensis and A. tzacatl is retrieved as monophyletic in molecular phylogenies (McGuire et al., 2014; Ornelas et al., 2014). Amazilia yucatanensis is cataloged as a Least Concern species due to its wide geographic distribution (International Union for Conservation of Nature [IUCN], 2018). It inhabits different types of environments, including seasonally deciduous tropical dry forest, sub-humid forest, bushes, and scrubs from almost sea level to 1200 m (Howell & Webb, 1995).
Three subspecies are currently recognized based on distribution and phenotypic differences, with belly from pale greyish buff to cinnamon (Howell & Webb, 1995; Figure 1A–C). Amazilia y. yucatanensis (Figure 1A) is distributed on the Yucatán Peninsula (Tabasco, Chiapas, Campeche, Quintana Roo, Yucatán) in Mexico, Petén of N Guatemala and N Belize. It inhabits lowland areas associated with seasonally deciduous tropical dry forest of the Yucatán Peninsula, and also occurs in secondary growth, clearings and forest edges. Amazilia y. cerviniventris (Figure 1B) is distributed in Chiapas, Veracruz, Puebla, Oaxaca, SE San Luis Potosí, and S Tamaulipas (Tampico), in Mexico. It inhabits arid as well as shady areas of humid forests, and thickets along the rivers (Schuchmann, 1999). Finally, A. y. chalconota (Figure 1C) is historically distributed from Tamaulipas, Nuevo León and San Luis Potosí in NE Mexico to S Texas in the US (Schuchmann, 1999), inhabiting open areas of pine-oak forest, semiarid coastal shrub, secondary growth, and patches of Texas palmetto (Sabal texana). Recent records suggest that A. y. chalconota has extended its distribution area northwards, from Texas to Louisiana or Florida (Brush et al., 2020; Chávez-Ramírez & Moreno-Valdez, 2015). These records are attributed to irregular migration in the autumn or winter northeast along the Gulf of Mexico coast, reaching Louisiana, Mississippi and Florida; rarely even farther north. In this area, the species is considered nonbreeding (Schuchmann, 1999), so the occurrence in this area where reproduction rate is negative was not considered for SDM (sink populations).
Figure 1.
Amazilia yucatanensis Subspecies and Current Geographic Distribution. The distribution areas were developed according to biogeographic regions proposed by Olson et al. (2001), and occurrence unique records used for modeling current species distribution. A: A. y. yucatanensis. B: A. y. cerviniventris. C: A. y. chalconota. D: Distribution map of A. yucatanensis subspecies. Different colors show the mask for each subspecies differentiated (BIOs value extraction); the sum of the masks represents the mask of the species (Figure S1) based on the ecoregions proposed by Olson et al. (2001). Photo by A: Pedro Yuit, B and C: Naturalista ( http://www.naturalista.mx/taxa/43155-RAmazilia-yucatanensis).
Data Collection and Study Area
Hummingbird occurrence data were assembled mainly from online records in the public databases Global Biodiversity Information Facility (GBIF.org (17 June 2020) GBIF Occurrence downloads, https://doi.org/10.15468/dl.erz6xy, https://doi.org/10.15468/dl.xqs645, https://doi.org/10.15468/dl.p6jwyk, https://doi.org/10.15468/dl.snztga), eBird (Sullivan et al., 2009), VertNet ( http://portal.vertnet.org/search?q=Amazilia±yucatanensis), and supplemented with our georeferenced records from field collection.
We selected records from museum collections located within a known distribution. Each database was downloaded and refined separately and in case of subspecies datasets some georeferenced records were eliminated when it was impossible to determine subspecies identity. We assembled data from the four databases (A. yucatanensis and its three subspecies). Georeferenced data were checked for errors and data consistency for geographic coordinates. The datasets were verified spatially to remove duplicate points using SDM Toolbox 2.2b in ArcMap 10.5 (Environmental System Research Institute [ESRI], 2016), excluding duplicate occurrence records or in close proximity to each other (ca. 1 km2) to reduce the effects of spatial autocorrelation. After careful verification of every data location, excluding duplicate occurrence records, we restricted the dataset to 310 unique presence records for the analysis for A. yucatanensis. Since the subspecies databases were downloaded individually to avoid identification errors at the limits of the distribution areas, the number of occurrence records decreased (n = 293), 81 for A. y. yucatanensis, 152 for A. y. cerviniventris and 60 for A. y. chalconota. These presence records were used to generate current and future SDMs.
We used species distribution modeling (SDM; Elith et al., 2011) to predict the distribution of A. yucatanensis and three subspecies (A. y. yucatanensis, A. y. cerviniventris, A. y. chalconota). An important step in ecological niche modeling is to delineate a realistic calibration region (‘M’, BAM diagram; Soberón & Townsend Peterson, 2005); that is, the set of sites accessible to a species over which models are calibrated (Atauchi et al., 2020; Barve et al., 2011; Freeman et al., 2019; Soberón & Townsend Peterson, 2005). In this study, we calibrated the SDMs for A. yucatanensis and the three subspecies (in practice a mask or GIS polygon) using a geographical clipping based on the ecoregions of eastern Mexico (Figure 1D and S1) proposed by Olson et al. (2001), representing potential boundaries on the landscape to dispersal (Barve et al., 2011). For A. yucatanensis and its subspecies, we buffered the ecoregions with highest concentrations of occurrence points (Townsend Peterson et al., 2011). Following calibration, we masked all environmental variables to the extent of the ‘M’ in ArcMap 10.5 (ESRI, 2016). For present projections, we use A. yucatanensis area M, and for future climate projections, we transferred our models to an extent of the Gulf of Mexico slope (Figure S1).
Environmental Data
To carry out SDMs environmental landscapes, we used 19 climatic variables summarizing data of precipitation and temperature (BIO1–BIO19 variables; Table S1) as climate layers from WorldClim 1.4 (Booth et al., 2014; Hijmans et al., 2005) at 1 km [ca. 30 sec of arc (arcmin)] spatial resolution with observed average values for mean annual temperature and total annual precipitation for the period 1961–1990. We converted all climatic variables in ASCII format with a 30-arc sec resolution (0.93 × 0.93 km = 0.86 km2 at the equator). All analyses, conversions and calculations were made in ArcMap.
To exclude highly correlated variables and multicollinearity, we ran a Pearson correlation test among the 19 climatic variables for each dataset using the “ggplot2” (Wickham, 2016) and “corrplot” (Wei & Simko, 2017) libraries in R 4.0.0 (R Development Core Team, 2017). When Pearson’s correlation coefficients were higher than 0.8 the variables were considered highly correlated, whereas variables with correlation coefficients less than 0.8 were selected to represent climatic limitations (Pearson et al., 2007). However, we are aware that hazards of multicollinearity to the recommended >0.7 threshold are possible (Dormann et al., 2013). After removing one of the highly correlated variables, six variables were used in the final analyses for both current and future scenarios of A. yucatanensis (BIO3, BIO4, BIO7, BIO13, BIO15, BIO17), six variables for A. y. yucatanensis (BIO1, BIO3, BIO4, BIO9, BIO13, BIO17), six variables for A. y. cerviniventris (BIO3, BIO4, BIO7, BIO11, BIO13, BIO15), and six variables for A. y. chalconota (BIO3, BIO4, BIO7, BIO9, BIO14, BIO17).
Selection of Climate Models to Estimate Future Distribution
We selected three different global climate models (GCMs) because predicted future species distributions may exhibit differences among them (Collevatti et al., 2013; Goberville et al., 2015; Heikkinen et al., 2006). To do this, we used the web application GCM compareR (Fajardo et al., 2020) evaluating the differences into the 32 GCMs projections from the Coupled Model Intercomparison Project Phase 5 (CMIP5) developed by the Intergovernmental Panel on Climate Change (IPCC, 2013). Then, we selected three GCMs that showed: 1) temperature and precipitation close to the average ensemble projection (the Community Climate System Model 4.0 [CCSM4]), 2) high temperature and low precipitation compared to the ensemble projection (the Met Office Hadley Centre Global Environment Model 2 - Earth System [HadGEM2-ES]) and, 3) high temperature and high precipitation compared to the ensemble projection (the Model for Interdisciplinary Research On Climate [MIROC5]) (see Borzée et al., 2019; Gent et al., 2011; Meinshausen et al., 2011). All models were developed for the years 2050 (average for 2041–2060) and 2070 (average for 2061–2080) using four contrasting climate change scenarios (Representative Concentration Pathways; RCPs) based on atmospheric CO2 trajectories: 2.6, 4.5, 6.0, and 8.5 (IPCC, 2013). The first scenario (RCP 2.6) represents an optimistic projection characterized by a very low concentration and emissions levels of greenhouse gases, whereas the last scenario (RCP 8.5) is a pessimistic projection with high levels of concentrations and emissions of greenhouse gases (Meinshausen et al., 2011; Moss et al., 2008; Riahi et al., 2011).
Spatial Model Procedures and Validations
We constructed SDMs with MaxEnt 3.4.1 (Phillips et al., 2017), a maximum entropy approach, using occurrence records and the 19 climatic variables. MaxEnt is a machine learning algorithm that allows SDMs to be generated using presence-only data, making it an effective tool to predict species distribution when obtaining presence-absence data is logistically impractical (Merow et al., 2013; Phillips et al., 2006, 2009, 2017). Although recent studies have shown variation in forecasted species distributions depending on the algorithm employed (Heikkinen et al., 2006; Qiao et al., 2015), we used MaxEnt over other available methods as it has been proved its high performance and suitability for presence-only data (Elith et al., 2011; Morales et al., 2017). Moreover, it has been reported that results from Bioensembles tend to overfit suitable areas, increasing the omission errors in the data (Prieto-Torres et al., 2020), Therefore, we discarded this last approach from our study.
Final models for A. yucatanensis and its subspecies were constructed using MaxEnt and were run with no extrapolation and no clamping to avoid artificial projections of extreme values of ecological variables (Elith et al., 2011; Owens et al., 2013; Prieto-Torres et al., 2020), while other MaxEnt parameters were set to default for convergence threshold (10−5) and 500 iterations, ensuring only one locality per grid cell. Here we assumed that the species was unable (by abiotic and/or biotic conditions) to disperse into those potential new areas and would inhabit only portions of the present distribution that remain habitable in the future. We also assumed no evolution in niche characteristics and do not consider specific interactions among species (Atauchi et al., 2020). The models were run with 10 cross-validation replicates with a 30% occurrence records for model evaluation, which is considered appropriate for estimating presence probability (Phillips et al., 2017), using the 10th percentile training presence logistic threshold (T10LT). The averages of all runs were used as final models, and jackknife analysis was used to determine the factors contributing to the greatest amount to habitat suitability (Borzée et al., 2019). The obtained maps were subsequently converted into binary presence-absence data and overlapping was performed in ArcMap 10.5 (ESRI, 2016). Occurrence data often exhibit strong spatial biases in survey effort meaning simply that some sites are more likely to be surveyed than others; such potential sampling biases are typically spatially autocorrelated and that might introduce bias in our results (Phillips et al., 2009). However, the occurrence records of hummingbird species from museums and collections are supported by thousands of occurrences data from databases such as eBird or Naturalista and supplemented with our georeferenced records from field collection.
The MaxEnt final models were evaluated by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC) curve (Elith et al., 2011), a threshold-independent statistic varying from 0 to 1, in which values around 0.5 represent distribution models no better than random and those around 1 represent a perfect fit between the observed and the predicted species distribution; acceptable models are those with > 0.7 AUC values (Phillips et al., 2006). However, several criticisms have been associated with this approach (e.g., Cobos et al., 2019; Lobo et al., 2008; Merow et al., 2014; Townsend Peterson et al., 2008), including, that the two error components (omission and commission) are inappropriately weighted equally. Therefore, the statistical performance of the models was also evaluated using the Partial-ROC test (Townsend Peterson et al., 2008). Within a value range from 0 to 2, values over 1 suggest a better performance than chance, by analyzing the presences versus the absence against the total area predicted by MaxEnt (Osorio-Olvera et al., 2020). We also used the true skill statistic (TSS), which is a threshold-dependent measure of model performance, to evaluate the accuracy of predictive maps generated by presence-only data (Allouche et al., 2006; Liu et al., 2013), where TSS values ranging between 0.4 and 0.8 are considered useful (Fielding & Bell, 1997). For each replicate, TSS was calculated using the T10LT and then TSS values were averaged among replicates.
To assess the extrapolation risk, we performed a Multivariate Environmental Similarity Surface (MESS) analysis to determine novel climatic conditions under future climate scenarios (Elith et al., 2006, 2010). MESS analysis compares the environmental similarity of variables and identifies areas where one or more environmental variables are outside of the training range using “ntbox” R package (Osorio-Olvera et al., 2020). Negative values indicate a novel climate, and the magnitude indicates the degree to which a point is out of range from its predictors. Positive values indicate climate similarity and are scored out of 100, with a score of 100 indicating that a value is entirely non-novel (Elith et al., 2010).
Lastly, the distributional areas of the species and subspecies for current and future projections were calculated (in pixels) for all different climate change scenarios (RCPs) for the years of 2050 and 2070 for global climate models. The final binary maps were projected at UTM 14 coordinates and the distribution area was calculated considering only the pixels with presence values. We then quantified the changes in the predicted surfaces, comparing all the current and future presence/absence maps of the species and the subspecies. All the presence/absence maps indexed the environmental suitability with values of 0 (unsuitable) and 1 (suitable). We converted the presence/absence maps to raster layers with float data-type using ArcMap, then we calculated the number of cells (pixels) among projected climatic extent using zonal statistics in spatial analyst tools in ArcMap. Finally, we converted the differences in the mean number of cells of potential habitats to surface area (km2).
Climatic Niches
To test the relationship between ecological niches and multivariate environmental and geographic spaces, we calculated the climatic niche breadth using NicheA 3.0 (Qiao et al., 2016), and assuming that the fundamental ecological niches of species are convex in shape (Soberón & Nakamura, 2009). The operationalized niches as Minimum-Volume Ellipsoids (MVE; Van Aelst & Rousseeuw, 2009) were then quantified, and similarities between ecological niche models for virtual species generated in NicheA were contrasted in terms of overlap in n-dimensional environmental spaces (Qiao et al., 2016). We extracted climatic data from the six climatic variables previously used in SDMs for each subspecies and performed a principal components analysis (PCA) with 10 variables total (BIO1, BIO3, BIO4, BIO7, BIO9, BIO11, BIO13–15, BIO16; Table S2). Using the first three axes, we created the environmental background, yielding more than 2 million points. The climatic niche breadth (MVE) of each subspecies was generated using the current time’s binary model.
Additionally, to reinforce the NicheA analysis, a comparative analysis of niche overlap was performed using ENMTools 1.3 (Warren et al., 2010). The assessment of niche overlap allows quantifying the niche shared by subspecies. Here, niche overlap between subspecies was computed by means of the Schoener’s D statistic directly from ecological niche space (Schoener, 1968; Warren et al., 2008). The value of D ranges between 0, which means that niches are completely dissimilar, to 1, which means that niches completely overlap (Broennimann et al., 2012).
Results
Current SDMs
The current distribution of A. yucatanensis was supported by high predictive power (mean ± SD; AUC = 0.906 ± 0.005, TSS = 0.654 ± 0.036; Figure 2A). The partial ROC test (1.889 ± 0.037) showed that models were statistically significant (p < 0.01). Thus, performance values for the model assessment approach indicated that the distribution model was statistically accurate. The SDM revealed that the distribution of suitable habitat for A. yucatanensis was continuous, with a range of 384,838.82 km2 (T10LT = 0.3685). The contributions of temperature and precipitation variables influencing performance of A. yucatanensis model were 48.7% and 43.4%, respectively. Based on the model’s predictions, A. yucatanensis occurred in an area with high variation in precipitation and warm temperatures (Figure S2).
In the subspecies analysis, the current distribution model of A. y. yucatanensis, A. y. cerviniventris and A. y. chalconota were also highly supported (mean ± SD; A. y. yucatanensis: AUC = 0.956 ± 0.004, TSS = 0.720 ± 0.083, partial ROC = 1.673 ± 0.277, p < 0.01); A. y. cerviniventris: AUC = 0.965 ± 0.003, TSS = 0.795 ± 0.055, partial ROC = 1.939 ± 0.022, p < 0.01); A. y. chalconota: AUC = 0.976 ± 0.004, TSS = 0.807 ± 0.098; partial ROC = 1.968 ± 0.008, p < 0.01). The distributions of all three subspecies were continuous, with a range of 162,793.7 km2 of suitable habitat for A. y. yucatanensis (T10LT = 0.4347; Figure 2B), 132,238.76 km2 for A. y. cerviniventris (T10LT = 0.2789; Figure 2C) and 117,488.04 km2 of suitable habitat for A. y. chalconota (T10LT = 0.2644; Figure 2D).
Temperature variables were the main factors influencing the model performance of A. y. yucatanensis (81.7%) and A. y. chalconota (81.8%), whereas the model performance of A. y. cerviniventris was influenced by both precipitation and temperature variables (62.7% and 23.1%, respectively). Based on model’s predictions, A. y. yucatanensis occurred in areas with high temperature and low precipitation (Figure S3), A. y. cerviniventris in areas with high precipitation (Figure S4), and A. y. chalconota in dry areas with low precipitation and stable temperature (Figure S5).
MESS analysis identified multiple areas where no analogs or novel climates were present, which varied among models, years, model types, specie and subspecies (Figure S14–S17). Dissimilarity values were relatively low within suitable areas for A. yucatanensis and A. y. chalconota. Nonetheless, novel conditions are not predicted in the most suitable areas for A. y. yucatanensis and A. y cerviniventris (Figure S14-S17).
Figure 2.
Predicted Current Distribution of Amazilia yucatanensis and Subspecies. A: Amazilia yucatanensis. B: Amazilia yucatanensis. yucatanensis. C: Amazilia yucatanensis cerveniventris. D: Amazilia yucatnensis chalconota. E: Environmental space in NicheA: Three-dimensional visualization of ecological niches in terms of principal components (PC); ecological niche models were displayed as Minimum Volume Ellipsoids used to illustrate limits of environmental distributions of subspecies (Green ellipsoid: A. y. yucatanensis; Mint ellipsoid: A. y. cerviniventris; Blue ellipsoid: A. y. chalconota). Gray dots represent the environmental background produced by the Principal Components Analysis (PCA) of the 10 selected variables (Table S2). F: Geographic polygons of tree subspecies. (Green A. y. yucatanensis; Mint A. y. cerviniventris; Blue A. y. chalconota; in red show overlap areas.
Future SDMs
The projections of the overlap SDMs for the 2050 and 2070 climate-change scenarios in A. yucatanensis showed that the current distribution of suitable habitat for the species increases (.10 to 9.11%) or decreases slightly (–2.19 to –12.16%), depending on the climate change scenario (Table 1, Figure 3). However, the current distribution suitable habitat for the species increases in most GCMs projections (Table S3, S4, S5; Figure S6 and S10). A little increase in the area is observed to the north (Nuevo Leon, Tamaulipas, Texas) and towards the coast (Veracruz), and a higher increase to SE Mexico (Oaxaca, Veracruz, Tabasco, Chiapas) and the Yucatán Peninsula (Campeche, Quintana Roo), Belize, Guatemala, and Honduras. Moreover, suitable area was predicted to the Florida Peninsula (Figure 3 and S6, S10, S14).
Table 1.
Predicted Increase or Decrease (%) in the Extent Of Suitable Areas (km2) for Amazilia yucatanensis and Its Subspecies, Summation of GCCs (CCSM4, HadGEM2-ES and MIROC5) and Climate-Change Scenarios (RCPs) for the Years 2050 and 2070.
The same response was observed for A. y. yucatanensis. The SDM projections of the overlaps for the 2050 and 2070 climate-change scenarios showed that the current distribution increases slightly (Table 1, Figure 4). However, the current distribution of suitable habitat increases or decreases depending on the RCP or GCM model (Table S3, S4, S5; Figure S7 and S11). The increase is showed to SE Mexico (Veracruz, Chiapas) and the Yucatán Peninsula (Campeche, Quintana Roo), Belize, Guatemala, and Honduras. The overlaps projected distribution of suitable habitat for A. y. cerviniventris showed decreases (–1.39 to –31.79%; Table 1, Figure 4). However, the current distribution of suitable habitat increases or decreases depending on the RCP or GCM model. A little increase in the area is observed to the north (Tamaulipas and Nuevo León), SE Mexico (Veracruz, Tabasco, Chiapas, Campeche, and Yucatán), Belize, Guatemala, Honduras, and the Florida Peninsula (Table S3, S4, S5; Figure S8 and S12). Lastly the current distribution of suitable habitat for A. y. chalconota is fragmented and decreases in Mexico (San Luis Potosí, Nuevo León, Tamaulipas) and Texas in all projections for the 2050 and 2070 climate-change scenarios (Table 1, S3, S4, S5; Figure 4, S9, S13).
Climate Niches
The ecological niche based on a minimum-volume ellipsoid varied among subspecies (Figure 2E), in which A. y. cerviniventris had the widest climatic niche and A. y. yucatanensis had the narrowest one. Niche overlap was observed among subspecies, in which the A. y. cerviniventris ecological niche contained all A. y. yucatanensis niche space and a portion of the A. y. chalconota niche space, but the ecological niches of A. y. yucatanensis and A. y. chalconota do not overlap (Figure 2E and F, Supplementary Material, Video 1). The PCA of raw GIS data indicated three main niche axes that together explained 91% of the variation (Table S2). The first niche axis (61%) was associated with isothermality (BIO3) and precipitation of the wettest month (BIO13). The second niche axis (23%) was associated with temperature annual range (BIO7) and mean temperature of coldest quarter (BIO11) while the third niche axis (7%) was associated with precipitation of driest quarter (BIO17) and annual mean temperature (BIO1).
Niche overlap results suggest a great variability in the environmental space inhabited by the three subspecies (Table S6) and match with the results obtained with NicheA analysis. Allopatric A. y. yucatanensis and A. y. chalconota occupy considerably different environmental niches (D = 0.021). While A. y. cerviniventris at the center of the species distribution occupies an environmental niche relatively more related to the A. y. yucatanensis and A. y. chalconota niches (D = 0.18 and D = 0.22, respectively; Table S6). By quantitatively assessing the niche overlap between subspecies (Schoener’s D statistic), our results suggest that the distribution of one subspecies cannot be implied by the distribution of another one.
Discussion
In this study, we assessed the possible impact of the global climate change on hummingbirds to taxonomic categories below the species level (i.e., subspecies). When modeling at the subspecies level, we found that climate change will modify the location and spatial extent of their suitable habitat, whereas the current distribution area remains more stable under future climate-change scenarios when modeling at the species level.
Our results showed high predictive values (AUC, TSS and partial ROC) for the current potential distributions of A. yucatanensis and its subspecies according to SDMs. While our results accurately predicted the known distribution for A. yucatanensis, they also suggest that the area of distribution of suitable habitat for the species will increase a little in the three GCMs and under four climate projections as compared to the more drastic increases or decreases observed when modeling was performed below the species level. These results suggest that partitioning the species occurrence data into subspecies and creating separate niche models for each subspecies best capture the relationship between environmental and geographical spaces and can give us a better understanding of the effects of future climate change (Mota-Vargas & Rojas-Soto, 2016; Qiao et al., 2016; Smith et al., 2019), particularly when modeling the distribution of widespread species whose distributions span large environmental or habitat variation (Gonzalez et al., 2011).
Figure 3.
Predicted Distribution of Amazilia yucatanensis Hummingbirds Under Climate-Change Scenarios (Overlap Global Climate Models), From Optimistic RCP 2.6 to Extreme RCP 8.5 Concentration Pathways Conditions, for 2050 and 2070.
Figure 4.
Predicted Distribution of Amazilia yucatanensis, Amazilia yucatanensis yucatanensis, Amazilia yucatanensis cerveniventris and Amazilia yucatnensis chalconota Under Climate-Change Scenarios (Overlap Global Climate Models), From Optimistic RCP 2.6 to Extreme RCP 8.5 Concentration Pathways Conditions, for 2050 and 2070.
Unlike expectations of species range shifts toward higher latitudes and elevations as temperatures rise (Buermann et al., 2011; Chen et al., 2011; Parmesan, 2006), when evaluating the response of A. yucatanensis to future climate-change scenarios the SDMs predicted an increased species range including a southward expansion but the expanded area was small as compared to the magnitude of increases or decreases in predicted areas observed for the subspecies. Increased suitable area was predicted for A. y. yucatanensis in most of the RCPs scenarios for 2050 and 2070 depending on GCMs, whereas the greatest decrease in distribution area was observed in A. y. chalconota. Although decreases were observed in A. y. cerviniventris in the RCPs for 2070, increases of suitable area were predicted in the RCPs for 2050. Altogether these results suggest regional variation in habitat suitability and subspecies differences in their vulnerability to anticipated climate change (D’Amen et al., 2013; Murphy & Lovett-Doust, 2007), which it would be scrutinized with evidence of genetic differentiation, morphological variation and phylogeographic structure among subspecies (Vásquez-Aguilar et al., unpublished data).
The increase in the predicted distribution of suitable habitat for A. y. yucatanensis might be a response related to historical changes in the distribution of the seasonally dry tropical deciduous forest on the Yucatán Peninsula due to climate changes during the Pleistocene, which increases in temperature favored the expansion of populations of different species and the area currently covered by seasonally dry tropical deciduous forest since the Last Glacial Maximum (Licona-Vera et al., 2018; Ortiz-Rodriguez et al., 2020). The seasonally dry tropical deciduous forest is globally considered the tropical forest most threatened due to deforestation and land conversion for agriculture, and by climate change caused by an estimated reduction of precipitation and a higher and more homogeneous temperature in its areas of distribution (Sánchez-Azofeifa et al., 2013). The potential impact of climate change on the seasonally dry tropical deciduous forest is likely to be regionally different, and changes due to warmer temperature and increased variation in precipitation will potentially alter its distribution and phenological patterns and, consequently, the distribution and seasonality of the animal species that inhabit it (Miles et al., 2006; Sánchez-Azofeifa et al., 2013). Future distribution changes of the Mexican seasonally dry tropical deciduous forest were forecasted by Prieto-Torres et al. (2016) considering two RPC scenarios for 2050 and 2070. The projection of climate-change scenarios showed increases in area of 3.0–10.0% and 3.0–9.0% for 2050 and 2070, respectively, with particularly increases in area and shifts into new and higher areas in the Tamaulipas-Veracruz region and the Yucatán Peninsula (Prieto-Torres et al., 2016). However, suitable areas for the seasonally dry tropical deciduous forest were not predicted in the more arid regions in northern Tamaulipas and Texas, as observed for A. y. chalconota. Even though our results show different scenarios when evaluated at the species level than at the subspecies level, it is difficult to predict how species or subspecies will respond to climates that do not exist at present because the adaptive potential is determined by evolutionary rate and ability to respond to environmental change, which depends directly on factors such as reproductive rates, geographic range size and dispersal ability of the species (Prieto-Torres et al., 2020). Changes in the physical variables along the distribution could constrain the spread of individuals, the different suitability values and fragmented populations could be interpreted as potential geographic and ecological barriers for the dispersal capacity of A. yucatanensis subspecies as shown in other bird species (Atauchi et al., 2020). Therefore, more studies analyzing alternative hypotheses to the niche conservatism for species are needed to measure this variation and, if possible, increase our understanding of dynamic responses of species and their adaptation for future climate change (Prieto-Torres et al., 2020).
Regarding the current climate niche, the three subspecies appear to inhabit distinct ecological niches. Occurrence data suggested niche overlap between A. y. yucatanensis and A. y. cerviniventris and between A. y. cerviniventris and A. y. chalconota, but no overlap was observed between the allopatrically distributed A. y. yucatanensis and A. y. chalconota (Figure 2). However, in the geographical space the overlap areas show that not all areas were accessible to all subspecies, indicating that there are geographical or ecological barriers preventing dispersal and limiting contact between subspecies populations (e.g., Tocchio et al., 2015). Given the variation in environmental conditions where the three subspecies occur, it is not surprising the low niche overlap particularly between allopatrically distributed subspecies (Table S6), indicating their different environmental constrains. These results are important because show the differences between niches of the three subspecies supporting the hypothesis that each subspecies will be affected differently by climate change.
Contrary to studies showing that habitat loss under future climate change scenarios correlates with climatic niche size (e.g., Ramírez-Preciado et al., 2019), our data showed that A. y. yucatanensis has the largest distribution and the distribution area of suitable habitat increased in all RCPs scenarios for 2050 and 2070 but the smallest niche breath (Figure 2). Although the niche of A. y. chalconota was broader than that of A. y. yucatanensis, the predicted area of suitable habitat was predicted to decrease in future RCPs scenarios. This suggests that A. y. yucatanensis with broader geographic distribution, as compared to the other subspecies, is more sensitive to precipitation and temperature variation (Figure S2–S4) than its broad range might otherwise suggest and, consequently, might be threatened by climate change.
Future geographical distribution might be related to the tolerance to a wide range of climatic conditions that facilitates the occupation of large geographic areas (Atauchi et al., 2020; Freeman et al., 2019; Slatyer et al., 2013). Although range-restricted species are predicted to be more vulnerable to global warming (Colyn et al., 2020; Urban, 2015), our results showed that the effects of climate change are not restricted to range-restricted or habitat specialist species (Dirnböck et al., 2011; Malcolm et al., 2006), and depends on multiple factors, as shown for A. y. yucatanensis and A. y. chalconota. In this study, the subspecies with strong cinnamon coloration (A. y. yucatanensis) showed the highest increase in future distribution area, while the predicted area for the subspecies with less cinnamon coloration (A. y. chalconota) decreased in all future climate-change scenarios. In the future, however, a large part of the ranges they currently occupy will not retain the set of climatic conditions suitable for the presence of at least one subspecies (A. y. chalconota). Ecological interactions might also influence changes in their geographic range, especially in regions in which these hummingbirds have a specialized interaction with plant species from which they obtain floral nectar (Correa-Lima et al., 2019; Greig et al., 2017; Junker et al., 2013; Ornelas et al., 2004). Lastly, populations of A. y. yucatanensis might expand under future climate change scenarios because new areas with suitable climatic conditions would be available for dispersal, as the predicted increase for the seasonally tropical dry forest on the Yucatán Peninsula (Prieto-Torres et al., 2016), whereas populations of the other subspecies would decrease or remain stable under future conditions (Buermann et al., 2011; Courter et al., 2013; Greig et al., 2017; Lara et al., 2012).
Even though the climatic conditions are not encouraging for some hummingbird species, the distribution ranges of some species are expanding due to urbanization and anthropogenic changes (Greig et al., 2017; Jones, 2011). Although this does not mean that the climatic and environmental conditions are at the best in urban areas, supplementary feeding has become popular particularly in the United States (Orros & Fellowes, 2015), thus causing northward expansion (Courter et al., 2013; Greig et al., 2017; Plummer et al., 2015) and increases in the relative abundance of A. yucatanensis in urban settings of southern Texas (Brush et al., 2020). However, a northward expansion in the distribution of A. y. chalconota into northern Texas, or in the distribution of the seasonally dry tropical deciduous forest, was not detected under future climate-change scenarios (Prieto-Torres et al., 2016; this study). Therefore, it is likely that the current northward expansion in the distribution of A. y. chalconota in the United States is due to urbanization and supplementary feeding, which would favor the presence of medium-sized hummingbirds with relatively long bills as shown in urban environments of central Mexico (Puga-Caballero et al., 2020).
Implications for Conservation
The determination of species distribution area is a complex task involving various theoretical aspects and perspectives at different scales (Gonzalez et al., 2011; Mota-Vargas & Rojas-Soto, 2016; Smith et al., 2019). Our results showed that depending on the taxonomical level (species or subspecies), future distribution projections can vary drastically, placing at risk a widespread species that is not considered vulnerable to the effects of climate change, when in fact lineages below the species level might be threatened by climate change. Our study might guide further modeling studies to analyze the impact of climate changes on hummingbird species and correctly choosing the taxonomic unit to be analyzed, especially when the taxonomic information or the phylogeographic patterns are not yet clear. The differences in the environmental niche shown in our study offer insights on the ecological niches of Amazilia yucatanensis subspecies and the interaction between environmental and geographical spaces. Conservation actions particularly of currently endangered hummingbird species must consider intrinsic requirements of subspecies and the main environmental drivers that shape their distributions.
Acknowledgements
We thank Mariana Hernández-Soto, Diego F. Angulo, Yuyini Licona-Vera, Andres E. Ortiz-Rodríguez and Jorge Antonio Gomez-Diaz for their advice on data analysis and Dolores Hernández-Rodríguez for useful comments. This work constitutes partial fulfilment of A.A.V-A. ’s doctorate in the Tropical Ecology program at the Centro de Investigaciones Tropicales (CITRO), Universidad Veracruzana.
Declaration of Conflicting Interests The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was financially supported by research funds from the Departamento de Biología Evolutiva, Instituto de Ecología, A.C. (INECOL) awarded to J.F.O. (20030/10563) and the Centro de Investigaciones Tropicales (CITRO), Universidad Veracruzana. Publication of this article was funded by funds from the Centro de Investigaciones Tropicales (CITRO), Universidad Veracruzana.
ORCID iDs Antonio Acini Vásquez-Aguilar https://orcid.org/0000-0002-4920-7278
Juan Francisco Ornelas https://orcid.org/0000-0002-1124-1163Flor Rodríguez Gómez https://orcid.org/0000-0001-8613-3418
M. Cristina MacSwiney https://orcid.org/0000-0002-9007-4622
Supplemental MaterialSupplemental material for this article is available online.
