Large-scale phylogeographical patterns and the underlying factors driving species divergence in Mesoamerica are poorly understood, but it is widely documented that tectonic events and Pleistocene climate changes play an important role in determining species diversification. As glaciations in Mesoamerica developed only around high mountains, one hypothesis is that the known effects of the Last Glacial Maximum on the geographical distribution and genetic diversity of bird populations, producing the contraction/ expansion latitudinal pattern observed in temperate bird species, should be largely undetected in resident bird populations inhabiting environmentally more stable habitats. To gain insight into the effects of Quaternary habitat and climate stability on the genetic diversity, we use ecological niche modelling and generalised linear modelling to determine the role of changes in habitat stability on the genetic diversity in eight widespread or range restricted hummingbird species. We found lesser changes in suitable habitat from past to present in most of the species than those predicted by palaeodistribution models at northern temperate regions. Contemporary seasonal precipitation, Quaternary habitat and climate stability had superior explanatory power, but the magnitude and directionality of their effects on genetic diversity varied between range-restricted and widely distributed species. We observed that the species studied have not responded equally to changes in climate stability in this complex region, suggesting that habitat differences and/or the altitudinal range of the hummingbird species influenced genetic diversity, and that the species-specific responses are not only linked to habitat stability in the region but also to contemporary seasonality associated with the availability of floral resources.
Introduction
The distributions of Nearctic and Neotropical biota come into contact in Mexico, a broad area of overlap known as the Mexican Transition Zone (Halffter 1987) or Megamexico (Rzedowski 1991). This complex area including the southern United States, Mexico and adjacent areas of Central America up to the Nicaraguan Depression is considered a critical area for understanding diversification in the region (Rodríguez-Gómez & Ornelas 2015). The mixed biota of this area has also received attention in conservation and biodiversity plans because of its high diversity and high degree of endemism of plant and animal taxa. The fragmented mountainous region spanning southern Arizona to western Texas in U.S.A. and northern Coahuila in Mexico south to the Isthmus of Panama encompasses two biodiversity hotspots (Stattersfield et al. 1998 , Myers et al. 2000, Mittermeier et al. 2005). The Madrean Pine-Oak Woodlands Biodiversity Hotspot includes Mexico's main mountain chains (Sierra Madre Occidental, Sierra Madre Oriental, Trans-Mexican Volcanic Belt, Sierra Madre del Sur, Sierra Norte de Oaxaca), as well as isolated mountaintop islands in Baja California and the Madrean Sky Islands in the U.S.A. (Fig. 1). The Trans-Mexican Volcanic Belt, which runs from west to east across central Mexico, and serves as bridge connecting the Sierra Madre Occidental and the Sierra Madre Oriental, is the highest mountain chain in the hotspot (Myers et al. 2000). This hotspot is characterized by pine-oak woodlands, ranging from monospecific stands of pines (Pinus), Douglasfirs (Pseudosuga), firs (Abies) or oaks (Quercus) to mixed stands of species in different regions. The pine-oak woodlands have an insular-type distribution by virtue of being surrounded by more extensive floristic provinces, which is generally tropical or arid (González-Elizondo et al. 2013). This feature is particularly noticeable in the northern Mexican Highlands and the Madrean Sky Island Archipelago.
The Mesoamerica Biodiversity Hotspot encompasses all subtropical and tropical ecosystems from central Mexico and south to the Isthmus of Panama, spanning most of Central America. In Mexico, the hotspot extends as far north as northern Sinaloa (Río Fuerte) on the Pacific Coast and as far as the middle of the Sierra Madre Oriental (west of Tampico) on the coast of the Gulf of Mexico (Myers et al. 2000, Fig. 1). The hotspot does encompass a complex mosaic of dry forests, lowland moist forests, and broad-leaved and coniferous forests at higher altitudes. In the midto the southern part of the hotspot, broad-leaved premontane and montane hardwood forests occupy steep and cloud-shrouded slopes.
The high degree of endemism in these two areas may reflect their imposing geological history and the glacial oscillations that occurred during the Pliocene and Pleistocene. A growing number of phylogeographical studies continue characterizing biotic transition patterns across Mexico, particularly as they relate to several major historical (late Neogene) events in this topographically complex landscape that appear to have had broad impacts on lineage diversification among co-distributed species (Castoe et al. 2009, Barber & Klicka 2010, Daza et al. 2010, Bryson et al. 2011, McCormack et al. 2011, Ornelas et al. 2013). The formation of the mountain systems, mountain ridges, and valleys probably caused the isolation of populations and, subsequently leading to their diversification and speciation (e.g. deep divergently phylogenetic lineages; Bryson et al. 2011, McCormarck et al. 2011). The glaciations that took place during the Pleistocene covering most of North America were likely another contributor to diversification in Mexico. Although the ice sheets only covered the highest mountain peaks of Mexico during these glaciations (Caballero et al. 2010), the area of suitable habitats in the lowlands potentially retreated by the cooler and drier conditions, and populations contracted to few environmentally stable refugia, where their long-term isolation may have driven genetic divergence. In contrast, the area of suitable habitats in the highlands experienced downslope expansion and displacement of montane colder-adapted species (pine forests and alpine vegetation) by several hundred metres (Caballero et al. 2010, Ramírez-Barahona & Eguiarte 2013). The displacement would have increased the area occupied by the species, creating opportunities for population expansion and increased gene flow (Ramírez-Barahona & Eguiarte 2013). Although Quaternary packrat middens from SW United States and paleontological studies from NW Mexico have documented movement of preferred habitats for species in the region (e.g. Holmgren et al. 2003, Coats et al. 2008), the predicted responses to Pleistocene climate changes, however, seem species-specific (e.g. Jezkova et al. 2015). Using phylogeographical data and species distribution modelling, studies of species investigated so far in the Mexican Transition Zone have described various responses to Pleistocene climate changes, including population fragmentation and displacement to southern refugia during glacial periods and postglacial range expansion to previously occupied northern areas after the glaciers receded (e.g. Milá et al. 2000, 2006, Smith et al. 2011, Pulgarin-R & Burg 2012, Manthey et al. 2014, Jezkova et al. 2015) or postglacial population expansion and the evolution of long-distance migration (Milá et al. 2006, Malpica & Ornelas 2014). Other species responses include southward displacement of populations with postglacial isolation and divergence (McCormack et al. 2008, Cavender-Bares et al. 2011, Gugger et al. 2011), downslope range expansion and population connectivity during the cold glacial periods, followed by the fragmentation into high altitude populations during warm interglacials (Ramírez-Barahona & Eguiarte 2014, Ornelas & Rodríguez-Gómez 2015), or no demographic changes with persistence across the geographical range throughout glacial cycles (Arbeláez-Cortés et al. 2014, Manthey et al. 2014, Ornelas & González 2014, Ruiz-Sanchez & Ornelas 2014) or range shifts through sites with suitable habitats (Rodríguez-Gómez & Ornelas 2015). Given the ambiguity of phylogeographical data for the last glacial period in the region, information on the genetic diversity of species, combined with palaeoclimate modelling, is needed to infer the effect of past climatic changes on poorly known species distributed in nonglaciated areas.
Here we used mitochondrial DNA (mtDNA) sequences and ecological niche modelling to test whether the genetic diversity of eight hummingbird species is affected by contemporary climate and Quaternary habitat and climate stability in the region. Hummingbirds (Trochilidae) are a particularly good model to investigate the effects of climate change on their distributional ranges because their populations are subject to different selection regimes imposed by differences in the abundance of floral nectar resources at different latitudes (Graham et al. 2009, Malpica & Ornelas 2014), phenological flowering shifts linked to climate changes (McKinney et at. 2012), and because they seem to follow a particular climate niche throughout the year (Malpica & Ornelas 2014). Lastly, Dalsgaard et al. (2011) used network analysis to determine the degree of specialization in plant-hummingbird mutualistic networks and showed that contemporary plant-hummingbird specialisation is determined by contemporary climate (precipitation) and climate-change velocity since the Quaternary. Previous phylogeographical analyses of mtDNA variation in the studied hummingbird species identified high intraspecific variation and genetic structure, with estimates of divergence times between genetic lineages spanning the Pleistocene (2.39-0.031 MYA; González et al. 2011, Arbeláez-Cortés & Navarro-Sigüenza 2013, Rodríguez-Gómez et al. 2013, Licona-Vera & Ornelas 2014, Malpica & Ornelas 2014, Rodríguez-Gómez & Ornelas 2014). These divergence time estimates suggest that more than one single glacial/interglacial cycle can be invoked to explain all of the recently described genetic divergence. Despite the fact that all species are hummingbirds and thus share somewhat similar dispersal characteristics and life history traits, the lineages included in this study do contain a diverse sampling of ecological groupings that vary in the size of their altitudinal and distributional ranges (restricted vs. wide distribution) and should be capable of providing a broad perspective of the effects of climate change on their genetic divergence in Mesoamerica.
Fig. 1.
Map of North America and Central America showing the Madrean pine-oak woodlands and Mesoamerica biodiversity hotspots (a) and hummingbird species distribution: (b) Phaethornis longirostris, (c) Campylopterus curvipennis, (d) Campylopterus rufus, (e) Amazilia cyanocephala, (f) Amazilia viridifrons, (g) Lampornis amethystinus, (h) Doricha eliza, and (i) Selasphorus platycercus. Dots correspond to occurrence unique records used for modelling current species distribution and the brown traces representing the distributional ranges according to Howell & Webb (1995).
Our aim is to use genetic diversity measures of mitochondrial DNA and ecological niche modelling to determine the link between population genetic diversity and Quaternary habitat and climate stability in eight hummingbird species. Specifically, we analysed geographical patterns of genetic diversity and tested whether hummingbird species distributed in areas with high habitat and climate stability sustain increased levels of local genetic variation (Dalsgaard et al. 2011, Ramírez-Barahona & Eguiarte 2014, Ortego et al. 2015). In comparison to species with a geographically restricted area of distribution, widely distributed species are not as constrained because they have broad environmental conditions that extend their geographical range and can migrate in order to track new floral nectar resources colonizing new adequate habitat more rapidly as floral nectar resources become available, and have larger population sizes and a longer time-lag to respond to climate changes. Therefore, we expect genetic diversity to be positively correlated with habitat stability and with contemporary and Quaternary climate stability in widely distributed species, and these relationships must be negative in range-restricted species.
Material and Methods
Study species
Our dataset includes mtDNA of eight hummingbird species, Phaethornis longirostris, Campylopterus curvipennis, Campylopterus rufus, Amazilia cyanocephala, Amazilia viridifrons, Lampornis amethystinus, Doricha eliza, and Selasphorus platycercus. The hermit P. longirostris occur at lowto-mid altitudes, from c. 100 to 1900 meters above sea level (m a.s.l.), and are mostly restricted to tropical dry forests and subtropical montane forests (Arbeláez-Cortés & Navarro-Sigüenza 2013). The emeralds have broad altitudinal ranges: C. curvipennis commonly found in cloud forests, humid evergreen forests and edges from 270 to 1590 m a.s.l. (González et al. 2011); C. rufus restricted to humid evergreen forests, forest edges and plantations from 900 to 2000 m a.s.l. (Howell & Webb 1995); A. cyanocephala commonly found in cloud forests, humid evergreen forests and second growth from 900 to 2200 m a.s.l. (Rodríguez-Gómez et al. 2013); and A. viridifrons restricted to arid to semihumid scrub, thorn forests, cloud forests, and lowland evergreen forest edges from 50 to 2100 m a.s.l. (Rodríguez-Gómez & Ornelas 2015). The mountain gems L. amethystinus occur at higher altitudes, from 900 to 3000 m a.s.l. (Howell & Webb 1995), and are mostly restricted to cloud forests and pine-oak forests. The bees D. eliza and S. platycercus occur at low coastal and dry deciduous forests (c. from sea level to 900 m a.s.l.) and high premontane areas and pine-oak forests (c. from 1600 to 3300 m a.s.l.), respectively (Licona-Vera & Ornelas 2014, Malpica & Ornelas 2014).
Sequencing
DNA samples were primarily obtained through our own collecting efforts between 2005–2014 mainly along the Sierra Madre Oriental, Sierra de Los Tuxtlas, Sierra de Miahuatlán, Sierra de Juárez, Sierra de Manantlán, and in the highlands of Chiapas and Guatemala. DNA data sets were obtained mainly as part of on-going projects investigating the evolutionary history of cloud forests in eastern Mexico (Ornelas et al. 2013). Sequence data from the mitochondrial NADH dehydrogenase subunit 2 gene (ND2) and subunit 5 gene (ND5), cytochrome b (cyt b), control region (CR), and ATPase 6 and 8 (ATP6, ATP8) were obtained from DNA extracted from tissue samples or tail feathers from birds captured in mist nets as fully described elsewhere (e.g. Malpica & Ornelas 2014). In most cases at least two mitochondrial fragments from the above were sequenced per species (Table 1). Samples were collected using required permits and approved animal welfare protocols. Methods for DNA preparation, PCR amplification, and sequence generation are given elsewhere (Cortés-Rodríguez et al. 2008, González et al. 2011, Arbeláez-Cortés & Navarro-Sigüenza 2013, Rodríguez-Gómez et al. 2013, Licona-Vera & Ornelas 2014, Malpica & Ornelas 2014, Rodríguez-Gómez & Ornelas 2014, 2015). All newly acquired sequences reported in this paper have been deposited at GenBank with the accession numbers and markers by species specified in Table 1.
Population genetic diversity
We calculated haplotype diversity (h) and nucleotide diversity (p) for each population as described elsewhere (González et al. 2011, Arbeláez-Cortés & Navarro-Sigüenza 2013, Rodríguez-Gómez et al. 2013, Licona-Vera & Ornelas 2014, Malpica & Ornelas 2014) using 16000 permutations in ARLEQUIN 3.01 (Excoffier et al. 2005). Samples (number of populations and individuals per species) are given in Table 1.
Species distribution modelling
We constructed a species distribution model (SDM, Elith et al. 2011) to evaluate the potential distribution of each hummingbird species under current climatic conditions and to predict where the suitable conditions resided during the Last Glacial Maximum (LGM, 21000-18000 years ago) and Last Interglacial (LIG, 140000-120000 years ago). We assembled a data set of occurrences from museum specimens according to Atlas de las Aves de México (Navarro et al. 2003) and obtained through http://vertnet.org and the Global Biodiversity Information Facility (GBIF, http://data.gbif.org/species/browse/taxon). The obtained data set was supplemented with our georeferenced records from field collection. After careful verification of every data location and removing duplicate occurrence records, we restricted the data sets to unique records for the analyses, leaving 50 unique presence records for P. longirostris, 282 for C. curvipennis, 185 for C. rufus, 130 for A. cyanocephala, 76 for A. viridifrons, 100 for L. amethystinus, 121 for D. eliza, and 201 for S. platycercus. These localities sample the entire distribution range of each species. Distribution records were input into and analysed with the maximum entropy algorithm in MaxEnt v. 3.2.2 (Phillips et al. 2006) using the dismo v. 1.0–5 package (Hijmans et al. 2014) in R v. 3.0.3 (R Development Core Team, http://www.r-project.org/) to infer the SDMs. Presentday temperature and precipitation data (BIO1-BIO19 variables) were drawn as climate layers from the WorldClim database at a spatial resolution of c. at 1 km2 in each raster (Hijmans et al. 2005) and masked the data to include only 125° to 75° W and 0° to 50° N. This spatial coverage encompassed the actual ranges occupied by the eight hummingbird species. To reduce the number of climatic variables and to minimize collinearity, a principal components analysis (PCA) using SPSS v. 17 for Mac (SPSS, Armonk, NY, U.S.A.) was carried out for each of the species data sets (Table 1). We then ran a correlation analysis to eliminate correlated environmental variables using the program PAST v. 2.12 (Hammer 2011). When the correlation coefficient was higher than 0.5 the variables were considered highly correlated, and for each pair of correlated variables we selected the ones with the highest loadings on the first PC components (Table 1). After removing the highly correlated variables, the remaining were used to generate the SDM models under current climate conditions using MaxEnt with the default parameters for convergence threshold (10-5) and 500 iterations, ensuring only one locality per grid cell. We evaluated model performance of each of the species using a random test percentage set to 30 % of the occurrence points to test the model and the other 70 % (training data) to generate the predictive models, and then the area under the receiver operating curve (AUC) of the threshold-independent receiving operating characteristic curve (ROC) calculated. Resulting species distributions under current climate conditions were projected onto past climate scenarios, at the LGM (at c. 2.5 km2) and LIG (at 1 km2), using the dismo package (Hijmans et al. 2014) in R. Past climate layers were also drawn from the WorldClim webpage for two LGM scenarios developed by the Paleoclimate Modelling Intercomparison Project Phase II (Braconnot et al. 2007): the Community Climate System Model (CCSM, Collins et al. 2004) and the Model for Interdisciplinary Research on Climate (MIROC, Hasumi & Emori 2004), and the LIG (Otto-Bliesner et al. 2006). The CCSM and MIROC climate models simulate different climate conditions, with cooler conditions assumed in CCSM than in MIROC (Otto-Bliesner et al. 2007).
Table 1.
Samples, DNA markers, and GenBank accession numbers of samples used in this study.
Quaternary habitat and climate stability and genetic Diversity
Indices of habitat stability were calculated from the present-day and past distribution models following Dalsgaard et al. (2011), Ramírez-Barahona & Eguiarte (2014), and Ortego et al. (2015). Three indices were calculated for LIG to present and LGM (CCSM) or LGM (MIROC) to present separately, by subtracting the corresponding values of each pixel for subsequent time periods: T1 = LIG to present, T2CCSM = LGM (CCSM) to present, and T2MIROC = LGM (MIROC) to present. Because the calculation was done from past to present, positive values indicated decreasing suitability from past to present.
To test the hypothesis that genetic diversity was negatively correlated with habitat stability in rangerestricted hummingbird species (P. longirostris, C. rufus, A. viridifrons, D. eliza) and positively associated in widespread species (C. curvipennis, A. cyanocephala, L. amethystinus, S. platycercus), we used a generalised linear model (GLM) approach in R v. 3.0.3 (R Development Core Team, http://www.rproject.org/). Haplotype diversity (h) and nucleotide diversity (p) were analysed as continuous response variables. Two separate full GLM models included the factor species or the factor distributional range as defined above, and the indices of habitat stability (T1, T2CCSM, T2MIROC) treated as continuous covariates and their interactions with species or distributional range. To evaluate the impact of Quaternary climate stability, the climate-change velocities (VT = mean annual temperature, VP = mean annual precipitation) were calculated for LIG to present and LGM (CCSM) or LGM (MIROC) to present. We used a method similar to those used by Loarie et al. (2009), Sandel et al. (2011), and Dalsgaard et al. (2011). The climatechange velocities (hereafter climate stability) for each species were simply calculated by subtracting the corresponding values of subsequent time periods for each pixel and divided by the elapsed time: LIG to present = 140000 years, and LGM (CCSM) and LGM (MIROC) to present = 20000 years. To test the association between genetic diversity and climate stability in range-restricted and widespread hummingbird species as defined above, we used the same GLM approach in R, in which a full GLM model included the factors species or their distributional range, the indices of climate stability (VPLIG to present, VTLIG to present, VPLGM (CCSM) to present, VTLGM (CCSM) to present, VPLGM (MIROC) to present, VTLGM (MIROC) to present) treated as continuous covariates and their interactions with distributional range or species factors, and haplotype diversity (h) and nucleotide diversity (p) as continuous response variables.
Lastly, to evaluate the impact of contemporary climate, mean annual temperature (BIO1), temperature seasonality (BIO4), mean annual precipitation (BIO12), precipitation seasonality (BIO15) and elevation (meters above sea level) were tested against genetic diversity measures (h, p) using linear regression models in R.
Results
Species distribution modeling
The SDMs yielded a good fit for the current geographical distributions of all eight-hummingbird species (Figs. S1 and S2). The area under the curve (AUC) values were 0.981, 0.983, 0.996, 0.976, 0.989, 0.979, 0.996 and 0.912 in P. longirostris, C. curvipennis, C. rufus, A. cyanocephala, A. viridifrons, L. amethystinus, D. eliza and S. platycercus, respectively. Although some over prediction was observed (Figs. S1 and S2) along certain geographic areas (Sierra Madre del Sur in C. rufus and A. viridifrons, north of the Rocky Mountains in S. platycercus and Caribbean islands in D. eliza that were predicted by the models with a moderately high probability), the predicted models closely match the current known distributions of all eight hummingbird species. The estimated potential distribution during the LGM indicates that some species expanded their distribution, with a lower probability applying CCSM than MIROC simulations (Figs. S1 and S2). A pattern of expansion/contraction of suitable habitat conditions is observed for populations of P. longirostris and A. cyanocephala, with expansion during LGM glacial period and contraction during interglacials (Fig. S1). However, our models also suggest some important changes in distribution since LIG. Briefly, the estimated potential distribution of L. amethystinus indicates that this species had a relatively stable distribution range during the last 140000 years (Fig. S2). For sedentary populations of S. platycercus, the models predicted that suitable habitat conditions for populations in eastern and southern Mexico expanded under LGM conditions along the Sierra Madre Occidental mountain range and south of the Rocky Mountains (Fig. S2). The comparison of the distribution models of C. curvipennis, C. rufus, A. viridifrons and D. eliza indicates that suitable conditions expanded during the LGM; range shifts of C. curvipennis towards the Yucatan Peninsula (Fig. S1), C. rufus north and south of current distribution (Fig. S1) and A. viridifrons towards the Trans-Mexican Volcanic Belt (Fig. S2) from LIG to LGM, and from central Veracruz for populations of D. eliza from LGM to present (Fig. S2).
Fig. 2.
Relationships between haplotype diversity (h) and nucleotide diversity (π) and indices of habitat stability, T1 = LIG to present, T2CCSM = LGM (CCSM) to present, and T2MIROC = LGM (MIROC) to present. See Table 2 for significant interaction terms. Relationships are shown for widely (thin lines, open circles) and range-restricted species (thick lines, filled circles).
Fig. 3.
Significant relationships between haplotype diversity (h) and nucleotide diversity (π) and climate stability measures, VT = mean annual temperature and VP = mean annual precipitation were calculated for LIG to present and LGM (CCSM) or LGM (MIROC) to present. See Table 3 for significant interaction terms. Relationships are shown for widely (thin lines, open circles) and range-restricted species (thick lines, filled circles).
Table 2.
The effects of distributional range or species and habitat stability from LIG to present and LGM (CCSM) or LGM (MIROC) to present on haplotype diversity (h) and nucleotide diversity (π) for widely distributed species (Campylopterus curvipennis, Amazilia cyanocephala, Lampornis amethystinus, Selasphorus platycercus) and range-restricted species (Phaethornis longirostris, Campylopterus rufus, Amazilia viridifrons, Doricha eliza) using generalised linear models. (a) by range, (b) by species. * P < 0.05, ** P < 0.01, *** P < 0.001.
Quaternary habitat and climate stability, and genetic diversity
For localities with three or more sampled individuals, haplotype diversity (h) ranged from zero to 1 and nucleotide diversity (π) from zero to 0.17 (Table 2). Haplotype diversity (h) was not associated with indices of habitat stability when analysed within a GLM framework (Table 2). However, haplotype diversity was significantly affected by range. The significant interaction between range and T1 and T2CCSM suggests that haplotype diversity (h) varied between widely distributed species and those with restricted areas of distribution (Table 2). For rangerestricted species, haplotype diversity was negatively associated with T1 (estimate ± SE: -4.47 ± 0.83, t = -5.37, P < 0.0001), and with T2CCSM (estimate ± SE: -2.38 ± 0.80, t = -2.96, P = 0.0037; Fig. 2). The range of the species had also a significant effect on nucleotide diversity (Table 2) but only significantly and negatively associated with T1 in range-restricted species (estimate ± SE: -0.03 ± 0.007, t = -3.62, P = 0.0004; Fig. 2). When full models were conducted by species, haplotype and nucleotide diversity measures were significantly affected by species identity (Table 2) but none of the interactions between indices of habitat stability and species were significant, except for the marginally significant interaction between species and T2MIROC (Table 2) in which haplotype diversity was positively associated with T2MIROC for the widely distributed S. platycercus (estimate ± SE: 2.25 ± 0.98, t = 2.29, P = 0.0239).
Table 3.
The effects of distributional range or species and climate stability (VT = total annual temperature, VP = total annual temperature) from LIG to present and LGM (CCSM) or LGM (MIROC) to present on haplotype diversity (h) and nucleotide diversity (π) for widely distributed species (Campylopterus curvipennis, Amazilia cyanocephala, Lampornis amethystinus, Selasphorus platycercus) and rangerestricted species (Phaethornis longirostris, Campylopterus rufus, Amazilia viridifrons, Doricha eliza) using generalised linear models. (a) by range, (b) by species. * P < 0.05, ** P < 0.01, *** P < 0.001.
Genetic diversity values were not significantly affected by climate stability measures, except VPLIG to present and VTLGM (MIROC) to present in haplotype diversity and nucleotide diversity, respectively (Table 3). However, range had a significant effect in both haplotype and nucleotide diversity values. In the full GLM models, for haplotype diversity (h) the interactions between distributional range and VTLIG to present and VPCCSM to present were statistically significant (Table 3). Haplotype diversity was positively related to climate stability from LIG to present (VT), and negatively from CCSM to present conditions (VP) in range-restricted species; the opposite pattern was observed for widespread species (Fig. 3). For nucleotide diversity (p), significant interactions with range were observed for VTLIG to present, VPLIG to present, VTMIROC to present and VPMIROC to present (Table 3). Nucleotide diversity was negatively related to VT and VP from MIROC to present conditions in range-restricted species; no significant relationships were observed for widespread species (Fig. 3). When models were conducted by species, genetic diversity values were again not affected by climate stability measures, except VTLGM (CCSM) to present in haplotype diversity (Table 3). However, haplotype diversity and nucleotide diversity were significantly affected by species (Table 3). The interactions between species and VTLIG to present and VTLGM (CCSM) to present for haplotype diversity and VTLIG to present, VPLIG to present, VTLGM (CCSM) to present, VPLGM (CCSM) to present and VTLGM (MIROC) to present for nucleotide diversity, were statistically significant (Table 3). Haplotype diversity was negatively associated with VTLIG to present and VPLIG to present in A. viridifrons (VT, estimate ± SE: -9.6e + 0.3 ± 2.8e + 03, t = -3.35, P = 0.0.0013; VP, estimate ± SE: -2.9e + 0.3 ± 9.3e + 01, t = -3.0, P = 0.0.0032), and positively associated with VTCCSM to present in A. viridifrons (estimate ± SE: 3.79e + 02 ± 1.3e + 02, t = 3.0, P = 0.0037) and C. rufus (estimate ± SE: 2.28e + 00 ± 8.9e - 01, t = 2.5, P = 0.0127). The same trend was observed in L. amethystinus but the effects were marginally significant (estimate ± SE: 6.7e - 01 ± 3.6e - 01, t = 1.85, P = 0.068). Nucleotide diversity was negatively associated with VTLIG to present and VPLIG to present in A. viridifrons (VT, estimate ± SE: -2.1e + 01 ± 2.4e + 00, t = - 8.88, P < 0.0001; VP, estimate ± SE: -1.7e - 01 ± 1.6e - 02, t = -10.81, P < 0.0001), positively associated with VTCCSM to present and VPCCSM to present in A. viridifrons (VT, estimate ± SE: 3.7e - 02 ± 1.2e - 02, t = 2.85, P = 0.005; VP, estimate ± SE: 1.1e - 03 ± 1.4e - 04, t = 7.96, P < 0.0001), and marginally significant for D. eliza (VT, estimate ± SE: 1.1e - 01 ± 6.3e - 02, t = 1.79, P = 0.077; VP, estimate ± SE: 1.5e - 03 ± 8.4e - 04, t = 1.81, P = 0.075) and L. amethystinus (estimate ± SE: 4.7e - 03 ± 2.6e - 03, t = 1.77, P = 0.081). From MIROC to present, nucleotide diversity and VT were positively associated in A. viridifrons (estimate ± SE: 1.4e - 01 ± 2.4e - 02, t = 5.66, P < 0.0001) and VT and VP were negatively and marginally associated in D. eliza (estimate ± SE: -1.3e - 01 ± 7.3e - 02, t = -1.77, P = 0.081) and L. amethystinus (estimate ± SE: -5.2e - 05 ± 3.0e - 05, t = -1.71, P = 0.092), respectively.
Haplotype diversity (h) and nucleotide diversity (π) were negatively associated with contemporary precipitation seasonality (BIO15) in widely distributed species (h, estimate ± SE: -9.1e - 03 ± 2.5e - 03, t = -3.57, P = 0.0005; p, estimate ± SE: -5.9e - 05 ± 2.7e - 05, t = -2.16, P = 0.033), and temperature seasonality (BIO4) and mean annual precipitation (BIO12) were positively associated with nucleotide diversity (π) in range-restricted species (BIO4, estimate ± SE: 7.6e -06 ± 2.5e - 06, t = 3.03, P = 0.0031; BIO12, estimate ± SE: 3.3e - 06 ± 9.9e - 07, t = 3.34, P = 0.0011). None of the other correlations were statistically significant (all P-values > 0.05; haplotype diversity, multiple R-squared: 0.27, F11, 106 = 3.61, P < 0.0002; nucleotide diversity, multiple R-squared: 0.32, F11, 106 = 4.62, P < 0.0001).
Discussion
Pleistocene climatic fluctuations had profound effects on the geographical distribution and genetic diversity of extant species, most notably in the temperate regions of northern Europe and North America (Hewitt 2000). During the LGM, a time when ice sheets covered much of present-day temperate Europe and North America, entire biotas were fragmented and displaced to more southern latitudes (i.e. southern refugia) (e.g. Pulgarin-R & Burg 2012). After the glaciers receded certain species expanded north and colonized these newly available areas by longdistance dispersal (i.e. postglacial expansion; Hewitt 1996, 2000, Milá et al. 2000, 2006, Pulgarin-R & Burg 2012, Malpica & Ornelas 2014). In contrast, glaciations in Mesoamerica developed only around high mountains (e.g. highest volcanoes along the Trans-Mexican Volcanic Belt, Caballero et al. 2010), not on the plateau surface, or occurred asynchronously, where the effects on the geographical distribution and genetic diversity of resident bird populations are predicted to be mild following the LIG or retreat of a previous glaciation (Rodríguez-Gómez & Ornelas 2015). Palaeoecological species distribution models for the LGM indicate that the eight hummingbird species did not experience the latitudinal pattern (contraction to southern refugia and northern postglacial range expansion) observed in temperate bird species with more northern distributions (e.g. Milá et al. 2000, 2006, Smith et al. 2011, Pulgarin-R & Burg 2012) except S. platycercus (see also Malpica & Ornelas 2014). The suitable habitats predicted by the LGM models suggest that the responses of these hummingbirds were idiosyncratic in the transition from the LIG to LGM and from LGM to present.
Current distribution model and its projection into the LGM and LIG scenarios indicate that the distribution of suitable habitat for L. amethystinus has remained relatively stable since the LIG, a pattern similar to that found for one of the cloud forest-adapted plant species that they pollinate (Ornelas & González 2014) and contrasting with those inferred for the other hummingbird species inhabiting other geographical regions or other habitats that have been more impacted by Pleistocene glaciations (e.g. northern postglacial range expansion in S. platycercus, Malpica & Ornelas 2014). For the other hummingbird species, regional distribution patterns seem to have been altered during the LGM, including the expansion of suitable conditions during LGM glacial period and contraction during interglacials observed for P. longirostris and A. cyanocephala and expansion under LGM conditions to different geographic areas in C. curvipennis, C. rufus, A. viridifrons and D. eliza. We analysed geographical patterns of genetic diversity and tested whether areas with Quaternary habitat and climate stability sustain increased levels of genetic diversity. We found that genetic diversity decreases with decreased habitat stability in rangerestricted species (Fig. 2), from LIG to present for both measures of genetic diversity and from LGM (CCSM) to present for haplotype diversity (h). None of the relationships changed significantly in widespread species. Therefore, our results seem to indicate that range-restricted species do not maintain stable geographical ranges over longer periods of time as compared to widely distributed species, and that genetic variation within populations (p) decreases as populations are subjected to climatically unstable habitat during the LIG or LGM (CCSM) to present transitions. Note that none of these associations were statistically significant under the MIROC, LGM model, perhaps because of decreased precipitation levels (Braconnot et al. 2007).
The GLMs models of climate stability also showed that the effects varied between widely distributed species and those with restricted areas of distribution. In the full GLM models, haplotype diversity was positively related with VT from LIG to present and VP negatively from CCSM to present conditions in rangerestricted species. However, when these GLM models were conducted by species, the associations between genetic diversity measures and those of climate stability from LIG or LGM to present conditions, the directionality and significance of these relationships seemed idiosyncratic species' responses. As evidenced from the different responses to changes in habitat stability from the LIG to present and from the LGM to present conditions, the eight species studied have not responded equally to changes in climate stability in this topographically complex region. That is, it is possible that habitat differences of the hummingbird species influenced the different patterns observed (Ghalambor et al. 2006, Bell et al. 2010, Qu et al. 2014). Unlike the situation in the Rocky Mountains for migratory populations of S. platycercus, in which the nearby LGM heavily ice-covered area presumably reduced suitable habitat to southern Mexico and Guatemala for sedentary populations (Malpica & Ornelas 2014), the glaciers in Mexico were restricted to high altitudes mainly along the Trans-Mexican Volcanic Belt (the low limit of glaciations being estimated at above 3000 m a.s.l.; Caballero et al. 2010) and, therefore, stable areas of suitable habitat for montane species existed on slopes below the glaciers throughout the Pleistocene. This hypothesis is supported by the wide altitudinal range (1600–3300 m a.s.l.) of sedentary S. platycercus populations restricted to premontane and pine-oak forests and implies a broad thermal tolerance and a considerable ecological plasticity (i.e. post-glacial long-distance migration, Ghalambor et al. 2006). The same explanation is invoked for the apparent in situ diversification of L. amethystinus, a cloud forest-adapted species with a wide altitudinal range (between 600 to 2900 m a.s.l.) reaching adjacent coniferous forests. Palaeoclimate modelling has indicated that the elevational range of cloud forests in the region was significantly displaced downslope and expanded into lowlands during the LGM but not during the LIG (Ramírez-Barahona & Eguiarte 2014). Therefore, it is likely that suitable cloud forest habitat was not significantly reduced and this allowed the persistence of L. amethystinus populations isolated at moderate altitudes with altitudinal migration as observed today. By contrast, the narrower altitudinal ranges (i.e. a more limited ecological tolerance) may have restricted P. longirostris, C. curvipennis, C. rufus, A. cyanocephala, A. viridifrons and D. eliza to habitats at lower elevation and a lower capability of dispersal during glacial periods.
Lastly, we found evidence that precipitation seasonality (BIO15) is the most important predictor of genetic diversity (h, p) in widely distributed species, whereas annual mean temperature (BIO1), temperature seasonality (BIO4), annual mean precipitation (BIO12) or altitude seemed unimportant in determining current genetic diversity. Both genetic measures (h, p) were negatively associated with precipitation seasonality in widely distributed species and temperature seasonality and mean annual precipitation were positively associated in range-restricted species. These results indicate that the decreased genetic diversity in range-restricted hummingbird species is not only linked to decreased Quaternary habitat and climate stability from past to present, as suggested by the negative relationship from LIG or LGM to present conditions in rangerestricted species (Tables 2 and 3), but also to precipitation seasonality (BIO15) which is associated with the current seasonality in the availability of floral resources in the region. Characterizing habitat and climate stability during LIG to LGM transitions linked to genetic diversity changes of migratory species (facing highest seasonality in the availability of floral resources) may represent a template for comparative phylogeography in the region.
Acknowledgements
We thank C. Bárcenas, R.A. Jiménez, A. Malpica, E. Ruiz-Sánchez, A. Vázquez, L.M. García-Feria, and S. Wethington for field and lab assistance; R.T. Brumfield, F. Sheldon and D. Dittmann (LSUMNS), J. Klicka and Sharon Birks (Burke Museum, University of Washington), R. Canales del Castillo (Universidad Autónoma de Nuevo León), J. Klicka (University of Nevada), C. Lara (Universidad Autónoma de Tlaxcala), A. Navarro-Sigüenza and B. Hernández (MZFC, UNAM) and P. Escalante (IB, UNAM) for providing tissue samples essential to this work; Utku Perktaş and Hakan Gür for the invitation to contribute at the special issue about phylogeography and ecological niche modelling in Folia Zoologica; and Roger Guevara and three anonymous reviewers for their valuable criticisms and suggestions. Our fieldwork was conducted with the permission of Mexico' Secretaría de Medio Ambiente y Recursos Naturales, Instituto Nacional de Ecología, Dirección General de Vida Silvestre (permit numbers: D00-02/3269, 02038/07, 01568/08, 02517/09, 07701/11) and the assistance of S. Wethington (Hummingbird Monitoring Network) and F. & T. Engelman (Rocky Mountains National Park) in the U.S.A. (Federal banding permit numbers 22902, 23315, 23315A). Financial support was provided by grants 25922-N, 61710, and 155686 from the Consejo Nacional de Ciencia y Tecnología (CONACyT) and research funds from the Instituto de Ecología, A.C. (20030/10563) awarded to JFO, AOR, FRG, SGL and YLV, are supported by doctoral scholarships (262563, 224634, 250318, and 262561, respectively) from CONACyT.
Literature
Appendices
Supplementary online materials
Fig. S1. Species distribution models for Phaethornis longirostris (a), Campylopterus curvipennis (b), Campylopterus rufus (c) and Amazilia cyanocephala (d) at Last Interglacial (LIG, 140-120 ka), Last Glacial Maximum (LGM, 21 ka; CCSM and MIROC models) and present climatic conditions, respectively. The predicted probability of occurrence is represented as an index of suitability between 0 and 1. Low values (< 0.8, orange to light grey) indicate that the conditions are unsuitable for the species to occur, whereas high values (> 0.8, red) indicate that the conditions are suitable (URL: http://www.ivb.cz/folia/download/ornelas_et_al._fig._s1.pdf).
Fig. S2. Species distribution models for Amazilia viridifrons (e), Lampornis amethystinus (f), Doricha eliza (g) and Selasphorus platycercus (h) at Last Interglacial (LIG, 140-120 ka), Last Glacial Maximum (LGM, 21 ka; CCSM and MIROC models) and present climatic conditions, respectively. The predicted probability of occurrence is represented as an index of suitability between 0 and 1. Low values (< 0.8, orange to light grey) indicate that the conditions are unsuitable for the species to occur, whereas high values (> 0.8, red) indicate that the conditions are suitable (URL: http://www.ivb.cz/folia/download/ornelas_et_al._fig._s2.pdf).
