Climatic dissimilarity associated with phylogenetic breaks

Abstract

Shared phylogenetic breaks often are associated with clear geographic barriers but some common phylogeographic breaks may lack obvious underlying mechanisms. A phylogenetic break involving multiple taxa was found in the Baja California Peninsula that was associated with a past sea barrier. However, geological evidence is lacking for this barrier’s past existence, and despite its current absence, the genetic breaks have persisted. This work explores the relationships between the current climatic niches for matrilineages of 11 vertebrate species as a possible explanation for the current geographic partitioning of matrilineages. We evaluated the climatic occupancy of each matrilineage through ecological niche models, background similarity, niche overlap, niche divergence, and Mantel tests. We found disparities in the climatic occupancy between north and south matrilineage of each taxon. Northern matrilineages are associated with lower temperatures and winter rains, while southern matrilineages reside in areas with higher temperatures and summer rains.

Discontinuidades filogenéticas a menudo están asociadas con barreras geográficas, pero barreras filogeográficas comunes pueden estar faltas de obvios mecanismos subyacentes. En la península de Baja California se ha encontrado un patrón de discontinuidades filogenéticas en múltiples taxones que está asociada con una barrera marina que existió en el pasado. Sin embargo, a pesar de que ni está actualmente presente la barrera, ni se ha encontrado evidencia de su existencia, las discontinuidades genéticas parecen persistentes. Este trabajo explora las relaciones entre los nichos climáticos actualmentes existentes para linajes maternos de 11 especies de vertebrados y las condiciones climáticas donde se encuentra cada linaje, como una posible explicación para el actual aislamiento entre linajes. Evaluamos las relaciones climatológicas de cada linaje materno mediante modelos de nicho ecológico, superposición de nichos, divergencia entre nichos, pruebas de similitud de fondo y pruebas de Mantel. Encontramos disparidades en las preferencias climáticas entre linajes maternos norte y sur de cada taxón. Los linajes del norte están asociados con temperaturas más bajas y lluvias invernales, mientras que los del sur prosperan con temperaturas más altas y lluvias de verano.

Baja California Peninsula, ecological separation, matrilineages, realized ecological niche

Linajes maternos, nicho ecológico, Península de Baja California, separación ecológica

Historical geographic processes such as population division, range expansion, and long-range colonization influence genetic variation (Templeton et al. 1995). Generally, researchers assume that genetic breaks are the result of geographic barriers to dispersal, cryptic species boundaries, or recent contact between historically allopatric populations (Neigel and Avise 1986; Nei and Takahata 1993; Wakeley and Hey 1997). However, in some cases, the geographic positions of phylogenetic breaks are inconsistent with known geographic barriers (Soltis et al. 1997; Gibbs et al. 2000; Bond et al. 2001; Puorto et al. 2001) or biogeographic limits (Burton 1998; Irwin et al. 2001a), yet nevertheless coincide with physical and behavioral differences between populations (Irwin et al. 2001a, b). Consideration therefore should be given to the possibility that genetic breaks not always are the result of geographic factors (Irwin 2002).

The Baja California Peninsula (BCP) of northwestern Mexico has a narrow shape (mean 80 km width × 1,220 km length) and nearly north-south alignment (Riddle et al. 2000a) demonstrated the presence of several biogeographic and climatic provinces in the BCP (see also Álvarez-Castañeda et al. 1995; Riddle and Hafner 2006; Garcillán et al. 2010; González-Abraham et al. 2010). The central part of the peninsula, known as a climate transition zone, is an area that marks the border between the south and north climatic provinces of the peninsula (Garcillán et al. 2010; González-Abraham et al. 2010). In that area, a recurrent genetic break has been recorded in many species, including plants (Nason et al. 2002; Clark-Tapia and Molina-Freaner 2003; Garrick et al. 2009), invertebrates (Crews and Hedin 2006; Garrick et al. 2013), and vertebrates (reptiles, birds, and mammals; Upton and Murphy 1997; Zink et al. 1997, 2001; Douglas et al. 2006, 2007; Álvarez-Castañeda and Lorenzo 2017).

The genetic breaks in the mid-peninsular region have been explained by the presence of an ancient seaway (Upton and Murphy 1997; Lindell et al. 2006, 2008). The seaway hypothesis is based on differences in fossils found in north and south portions of the peninsula (Helenes and Carreño 1999) and on topographic analysis (Dolby et al. 2015). However, the existence of a seaway is debated because there is no geological evidence to support it (Brusca et al. 2005; Lindell et al. 2006; Garrick et al. 2009; Dolby et al. 2015). This seaway apparently retreated about one million years ago (mya), as estimated based on inter-population divergences (Upton and Murphy 1997; Holt et al. 2000; Hafner and Riddle 2005; Lindell et al. 2005; Munguía-Vega 2011). One million years is enough time for the genetic break to disappear, particularly in species with high dispersion capacities such as birds (Irwin 2002). However, there are other factors that can maintain phylogenetic divergence once the geographic barrier has disappeared, including reproductive isolation among lineages and ecological or morphological divergences between groups (Irwin 2002). For instance, Giarla et al. (2018) demonstrated that ecological distances better explain genetic distances than do geographic distances among populations of an Asian island rodent species.

Maintenance of a genetic discontinuity without geographical barriers would not be rejected by the following hypotheses: (i) currently, there are no geographic barriers in the mid-peninsular region that would restrain gene flow in dispersive taxa; (ii) matrilineages of different species with genetic breaks show a contact zone (Álvarez-Castañeda 2007; Álvarez-Castañeda et al. 2009; Álvarez-Castañeda et al. 2010); and (iii) the same geographic area constitutes the distributional limit of many species, including low and medium-vagility species (e.g., gophers and birds, respectively; Álvarez-Castañeda and Patton 1999; Rios and Álvarez-Castañeda 2010). These findings suggest that the genetic divergence between BCP matrilineages is maintained by factors other than a geographic barrier.

All these conditions give rise to the following questions. Do matrilineages inhabit different climatic conditions? Are matrilineages of species with genetic breaks associated with external barriers to hybridization, such as the disparities in their realized ecological niches? To address these questions, we tested the following hypotheses: (i) the realized ecological niches of northern and southern vertebrate lineages are different; (ii) realized niches of each matrilineage have remained geographically separated between areas with different climatic conditions; and (iii) the current physical characteristics of the habitat in the area of contact between matrilineages do not restrain the potential reconnection between matrilineages. The objective of this work was to investigate an alternative hypothesis focused on realized ecological niches and their relationship with the environment to explain the current maintenance of genetic divergence in different vertebrates of the BCP.

Materials and Methods

We focused on groups of BCP fauna with available occurrence data that have been the subject of phylogeographic analyses (Supplementary Data SD1). We selected 11 species: 6 mammals (Ammospermophilus leucurus and Otospermophilus beecheyi [Rodentia: Sciuridae]; Dipodomys merriami and Chaetodipus fallax [Rodentia: Heteromyidae]; Thomomys nigricans [Rodentia: Geomyidae]; and Lepus californicus [Lagomorpha: Leporidae]; Whorley et al. 2004; Álvarez-Castañeda 2007; Rios and Álvarez-Castañeda 2007; Trujano-Álvarez and Álvarez-Castañeda 2007; Álvarez-Castañeda et al. 2009, 2010; Rios and Álvarez-Castañeda 2010; Álvarez-Castañeda and Cortés-Calva 2011; Trujano-Álvarez and Álvarez-Castañeda 2013; Álvarez-Castañeda and Lorenzo 2017); 3 reptiles (Crotalus ruber [Squamata: Viperidae], and Uta stansburiana and Urosaurus nigricaudus [Squamata: Iguania]; Upton and Murphy 1997; Aguirre et al. 1999; Lindell et al. 2008); and 2 birds (Campylorhynchus brunneicapillus [Passeriformes: Troglodytidae] and Auriparus flaviceps [Passeriformes: Remizidae]; Zink et al. 2001; Table 1; Supplementary Data SD2). For each species, we divided occurrence data according to matrilineages for model calibration. Our hypothesis focuses on the BCP sensu stricto, that is, insular records were not considered.

Table 1.

Number of presence records for each lineage of vertebrates with genetic breaks in the Baja California Peninsula and values associated with these temperature and precipitation records.

Vertebrates with a genetic break in the middle of the Baja California Peninsula	Number of records		Temperature (°C)						Precipitation (mm)
			Northa		Southa		Differenceb		Northa		Southa		Differenceb
	Northa	Southa	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max
Mammals
Ammospermophilus leucurus	106	74	13.0	23.7	17.1	24.3	4.1	0.6	60.0	395.0	66.0	518.0	6.0	123.0
Otospermophilus beecheyi	28	16	7.2	19.1	16.5	22.5	9.3	3.4	115.0	769.0	119.8	330.0	4.8	−439.0
Dipodomys merriami	228	114	12.1	23.7	17.1	23.9	5.0	0.2	55.0	447.0	45.0	499.0	−10.0	52.0
Chaetodipus fallax	135	30	12.9	22.4	19.2	21.4	6.3	−1.0	85.0	406.0	61.0	123.0	−24.0	−283.0
Thomomys nigricans	106	150	8.0	22.9	13.7	24.3	5.7	1.4	53.0	735.0	66.0	712.0	13.0	−23.0
Lepus californicus	61	52	9.8	23.7	18.2	24.2	8.4	0.5	60.0	635.0	63.0	430.0	3.0	−205.0
Birds
Campylorhynchus brunneicapillus	54	94	13.5	22.8	17.1	24.2	3.6	1.4	53.0	358.0	65.0	518.0	12.0	160.0
Auriparus flaviceps	29	53	14.6	22.8	20.0	24.2	5.4	1.4	53.0	327.0	61.0	422.0	8.0	95.0
Reptiles
Urosaurus nigricaudus	58	150	8.6	22.8	14.6	24.3	6.0	1.5	55.0	707.0	62.0	649.0	7.0	−58.0
Uta stansburiana	316	131	8.6	23.6	20.4	24.3	11.8	0.7	58.0	126.0	44.0	451.0	−3.0	325.0
Crotalus ruber	124	92	15.1	23.6	14.5	24.4	−0.6	0.8	61.0	368.0	47.0	659.0	−14.0	291.0

^aMatrilineages.

^bDifferences: southern values minus northern values.

First, we explored the climatic variation in the area each vertebrate matrilineage currently occupies. Second, we characterized and modeled the realized ecological niches for each matrilineage under two climatic scenarios (warm: present; cold: Last Glacial Maximum [LGM]: ~21,000 years ago]). Third, we evaluated the potential correlation between geographic and ecological conditions of the realized niches for each matrilineage in the mid-peninsular area where the genetic breaks of multiple taxa are found.

Climatic variation between matrilineages.

We obtained occurrence data for reptiles and birds from the VertNet database, supplemented with data from the literature (Upton and Murphy 1997; Aguirre et al. 1999; Zink et al. 2001; Lindell et al. 2008); records for mammal species were gathered from the mammal collection of the Centro de Investigaciones Biológicas del Noroeste (CIB) and from the literature (Supplementary Data SD2).

Occurrence records were screened according to various criteria to avoid potential biases resulting from spatial sampling (Boria et al. 2014). First, we removed duplicate records; second, we removed poorly georeferenced records (mismatches between actual coordinates and recorded location data, or insufficient information for confirming that the georeferencing was correct); finally, we established 1 km as the minimum acceptable distance between records of each species using the “ecospat” package in R (Table 1; Di Cola et al. 2017; R Core Team 2020). These filters allow us to maximize the information obtained on ecological niches of the matrilineages, reducing the sampling effects of the geographic space (Varela et al. 2013; Boria et al. 2014).

We obtained values for two climate variables: mean annual temperature (°C) and mean annual precipitation (millimeters [mm]), from the WorldClim database (Hijmans et al. 2005) for each locality. These two variables were selected because they are the climate variables usually employed to measure the climate at the Earth’s surface (Hartmann 2016). We calculated minimum, maximum, and median values for each variable for each matrilineage.

Ecological niche modeling.

We created ecological niche models (ENMs) for each matrilineage. ENMs aim to identify areas of suitable climatic space given a set of known localities within a defined model area (Guisan and Zimmermann 2000; Loyn et al. 2004; Soberón and Nakamura 2009). The calibration area was unique for each matrilineage, considering its current known distribution according to publications where the genetic divergence of each species was reported (Anderson and Raza 2010; Barve et al. 2011; Supplementary Data SD3). However, not all the potential environmental combinations that a species can tolerate (i.e., the fundamental niche = N_F) are available always in a local environment. Therefore, the set of favorable conditions for a species that are available in the environment is known as realized niche (N_R) and is what will be analyzed with ENMs (Hutchinson 1957; Peterson et al. 2011; Soberón and Peterson 2011).

ENMs were based on the 19 bioclimatic variables in the WorldClim data set, at a spatial resolution of 30 arc-seconds (~1 km). These variables represent summaries of means and variation in temperature and precipitation (Hijmans et al. 2005). For each matrilineage, we selected one set of variables from the WorldClim data set. We selected each set according to the following steps. First, we evaluated the correlations between the 19 variables with a Pearson’s correlation test. For each pair of variables with R² >0.8, one of the two was selected (Dormann et al. 2013), namely the one most relevant to the biology of the species and with a greater contribution to the model, as measured by the jackknife test carried out in MaxEnt (Phillips et al. 2006). The Pearson’s correlation test was carried out in R (R Core Team 2020).

We built niche models using MaxEnt 3.4.1 presence-background data technique (Phillips et al. 2006), which estimates the ecological niche of a species by finding the distribution of maximum entropy subject to the constraints of observed presences in the environmental space (i.e., the distribution that is most spread out, or closest to uniform; Phillips et al. 2006). We built separate models for each matrilineage of each species. For model calibration, records for each matrilineage were split into 75% for calibration and 25% for model evaluation, using the geographically structured method, dividing the calibration area into two sections by latitudinal gradient (Radosavljevic and Anderson 2014). MaxEnt needs a collection of pixels from the predefined study area with its associated climatic conditions, which is typically called the background (Elith et al. 2010). We used a random sample of 50,000 pixels through the study area delimited for each matrilineage (Supplementary Data SD3).

Aiming to evaluate the performance and control the complexity of the models used, these models were constructed using different types of feature classes (FCs) and regularization multiplier (RM) values. The first function refers to the algorithm applied to variables by MaxEnt; the second imposes a function penalty to avoid overfitting the models (Phillips et al. 2006; Phillips and Dudík 2008). We used 29 different combinations of FCs (linear [“L”], quadratic [“Q”], product [“P”], threshold [“T”], hinge [“H”], LQ, LP, LT, LH, QP, QT, QH, PT, PH, TH, LQP, LQT, LQH, LPT, LPH, QPT, QPH, QTH, PTH, LQPT, LQPH, LQTH, LPTH, and LQPTH), along with RM values ranging from 0.1 to 1 (with 0.1-point intervals) and from 2 to 5 (with 1-point intervals), and two sets of previously selected variables. All tuning steps were implemented using the “Kuenm” package in R (Cobos et al. 2019). Model performance was evaluated using yield (omission rate [OR] ≤ 10%) and model complexity (Akaike information criterion [AIC]). We selected those models that fulfilled the following: (i) OR ≤ 10% and (ii) delta AIC ≤ 2 (as a complexity threshold in the models, with values closer to 0 indicating less complexity; Warren and Seifert 2011). We calibrated the final models for each matrilineage using the MaxEnt parameterizations selected (Supplementary Data SD4). The final models were converted into binaries considering a 10th training presence (Araújo and Peterson 2012), which accepts up to 10% omission of the records in model training (Araújo and Peterson 2012; Escalante et al. 2013). R code is available in the Github repository (https://github.com/Leticia9/ENM-BCP).

Transfers to the cold-weather scenario.

To infer how the availability of past environments might have influenced phylogeographic patterns, we transferred the ENMs to the climatic conditions of the LGM. We used three Global Circulation Models (GCMs): Community Climate System Model (CCSM), Model for Interdisciplinary Research on Climate (MIROC), and Max Planck Institute for Meteorology (MPI). We used different GCMs because each provides different climatic parameterizations for the past and, consequently, alternative hypotheses of potential past distributions (Varela et al. 2015). To carry out the respective interpretations for transfers, we calculated the Mobility-Oriented Parity (MOP) to identify areas with non-analog conditions between calibration and transfer regions. MOP estimates Euclidean distances between calibration and transfer regions in multivariate environmental spaces; subsequently, MOP identifies and excludes those values that fall outside the environmental range of the calibration region; this test was implemented in the “ntbox” package in R (Elith et al. 2010; Owens et al. 2013; Osorio-Olvera et al. 2020).

Quantification of ecological niches between matrilineages.

Comparison tests were undertaken of the realized ecological niche between matrilineages in the geographic space (G-space) and ecological space (E-space). To evaluate the environments where each matrilineage is located, we carried out “Background Similarity Tests” (BSTs) in geographic dimensions (Warren et al. 2008, 2009). This test explores whether ENMs from populations with partially or totally non-overlapping distributions are more different from each other than expected by chance (H₀ = similarity). That is, it asks whether the ENM of a population predicts the presence of the other better than expected by chance (Warren et al. 2008). It is based on the generation of random points through the accessible area of a species or population in numbers equal to the amount of record data available for that matrilineage (Warren et al. 2009). The overlap between each model pair is calculated by pooling empirical occurrence points and randomizing their identities to produce two new samples with the same numbers of observations as empirical data (Warren et al. 2009). We measured the overlap between predictions of habitat suitability from the models using the Schoener’s D index (Schoener 1968) and the I index (modified Hellinger distance; Warren et al. 2008). The null hypothesis (H₀ = similarity) is rejected if D or I values fall below the 5th percentile in the random-replicate distribution of similarity values (Glor and Warren 2010). Tests were undertaken between matrilineages of the species studied using the “ENMTools” package in R (Warren et al. 2010; R Core Team 2020). We carried out 100 replicates to build the null distribution, and sample sizes for background points matched the number of points available for the other matrilineages. As background, we used the first five components of a principal component analysis of the 19 climatic variables of WorldClim, which account for 95% of the total variance. A Broken Stick test indicated the first three principal components (PCs) were significant, accounting for just under 87% of the cumulative variance; we used the first five PCs as a means to encompass 95% of the variance. R code is available in the Github repository (https://github.com/Leticia9/ENM-BCP).

To evaluate the equivalence and similarity between matrilineages, we carried out the “Niche Overlap Test” (NOT) and “Niche Divergence Test” (NDT) on the ecological dimensions (Brown and Carnaval 2019). The NOT estimates the similarity between the occupied niches of a species; it considers the total accessible environmental space within the geographic range of the species. When the NOT indicates significant differences in the total environmental space occupied by each of the two species, there is support for the hypothesis that they currently occupy different niches. However, this finding alone is insufficient to state whether niches have diverged or niche similarity (or lack thereof) is due to other causes. The NDT estimates the portion of the accessible environmental space that is shared by the two populations; it allows us to estimate the degree of similarity in the occupied niches of two species sharing a common environmental background. If the NDT yields a significant result, it indicates that the niches of two species sharing a common accessible environmental space are not equivalent, supporting the hypothesis that their fundamental niches are the result of divergent evolution. Both tests characterize the E-space as two axes of a PC analysis of environmental variables across the distribution range of both species. The difference between these approaches is that the NOT quantifies the total area of the two species, whereas the NDT quantifies only the portion of the niche that is shared by both species. Tests were carried out between matrilineages of the species studied using the “Humboldt” package in R (Brown and Carnaval 2019; R Core Team 2020). These tests calculate the degree of similarity between the occupied niches of two species by calculating the modified Hellinger distance (I) and Schoener’s D and compare these indices with those obtained when the occurrences of the two species are determined after resampling. The R code is available in the Github repository (https://github.com/Leticia9/ENM-BCP).

Relationship between physical and ecological conditions in the mid-peninsular region.

We explored potential correlations between landscape and ecological characteristics for each matrilineage in the area of genetic breaks of the vertebrates studied, by applying a Mantel test between resistance distance and ecological conductance measures (Mantel 1967). We first built distance matrices with physical resistance measures based on altitude and ecological conductance measures; ecological distance matrices were constructed for each matrilineage (Supplementary Data SD5).

Matrices were constructed using the Circuitscape software (Shah and McRae 2008). Circuitscape uses raster layers to measure friction and the random walk algorithm to estimate how the landscape might either restrain or facilitate the movement between nodes. For the landscape resistance matrix, we used the elevation raster (http://srtm.csi.cgiar.org/) as the friction layer and presence records in the area of genetic breaks for northern and southern matrilineages of each species as nodes (Supplementary Data SD5). For the ecological conductance matrix, we used the niche models of each matrilineage as the conductance layer and presence records in the area of genetic breaks for northern and southern matrilineage of each species studied as nodes (Supplementary Data SD5).

With the data from distance matrices, the Mantel test was run using the “vegan” package in R (Oksanen et al. 2019) through the Pearson’s correlation with 999 permutations. The Mantel test compared physical resistance distance with ecological conductance distance for each matrilineage (the R code is available in the Github repository: https://github.com/Leticia9/ENM-BCP). In the Mantel test, H₀ assumes that the response variable “Y” is not linearly correlated with the variable “X”. The objective, therefore, is to evaluate whether the association (positive or negative) is more robust than expected at random (Mantel 1967; Olano and Luzuriaga 2008).

Results

Climatic variation between matrilineages.

We obtained 2,201 occurrence data points for 11 vertebrate species from collections (CIB) and databases available online (VertNet; Supplementary Data SD2); the number of records per species and clade are provided in Table 1. The minimum and maximum values of temperature associated with presence records differed between matrilineages. Northern matrilineages displayed lower maximum and minimum temperature values relative to southern matrilineages (Table 1). The greatest difference between minimum temperature values was found for U. stansburiana, and the greatest difference between maximum temperature values between matrilineages corresponded to O. beecheyi.

For precipitation, the northern matrilineages of U. stansburiana, Cr. ruber, D. merriami, and C. fallax showed minimum precipitation values higher than minimum precipitation values associated with the presence of their southern matrilineages. On the other hand, the northern matrilineages of U. nigricaudus, O. beecheyi, T. nigricans, and L. californicus were associated with higher maximum precipitation values than their southern counterparts (Table 1; Supplementary Data SD6).

Ecological niche modeling.

We created 812 candidate models for each matrilineage of each species studied. The optimal model settings selected for each matrilineage were statistically better than the default settings. Niche models under current warm-weather conditions attained good performance: all models obtained less than 10% OR and a delta AIC equal to zero (see Supplementary Data SD4). The variables with the greatest contribution to mammalian niche models were different. For birds, precipitation variables made the greatest contribution for northern matrilineages, while temperature variables were most important for southern matrilineages (Supplementary Data SD4). The opposite occurred for reptiles, that is, temperature variables were most important for northern matrilineages, and precipitation variables for southern matrilineages (except for U. stansburiana; Supplementary Data SD4).

ENMs for mammals were consistent with the known distribution of both matrilineages. Spatially, the niche models of O. beecheyi, C. fallax, and T. nigricans showed a separation between the areas of climate suitability of their respective matrilineages (Fig. 1; see Supplementary Data SD7). Transfers to cold-weather scenarios followed the same pattern as for reptiles and birds. Northern matrilineages were able to find large areas of climatic suitability; this was not the case for southern matrilineages, however, because areas of climatic suitability were restricted to the southern tip of the peninsula (Fig. 1; Supplementary Data SD7). In the particular case of D. merriami, no sites of climatic suitability were observed for the southern matrilineage under two transfer scenarios (CCSM and MPI). On the other hand, the southern matrilineage of C. fallax shows large areas of climatic suitability. However, the MOP analysis indicates a low environmental similarity between calibration and transfer areas, so these areas may represent an extrapolation of the model (Supplementary Data SD7 and SD8).

Fig. 1.

Representative examples of ecological niche models (ENMs) for three vertebrate taxa (mammals, reptiles, and birds), and their respective transfers to the cold-weather scenarios (Last Glacial Maximum, ~21,000 years ago), with the three Global Circulation Models used (Community Climate System Model [CCSM], Model for Interdisciplinary Research on Climate [MIROC], Max Planck Institute for Meteorology [MPI]), under which sea levels would have been lowered, hence the larger surface areas shown for the Baja California Peninsula (BCP). These examples represent the pattern of divergences between the realized ecological niches of northern and southern matrilineages of species with genetic breaks in the BCP. Lines represent areas with non-analogous environments (na.env) between calibration and transfer areas; vertical for northern matrilineages and horizontal for southern matrilineages.

ENMs in the current warm scenario for birds (C. brunneicapillus and A. flaviceps) were consistent with the known distribution of each matrilineage (Fig. 1; Supplementary Data SD7). The cold (past) weather transfers showed suitable area for the northern matrilineages was much larger, in geographic extension, than for the southern matrilineages (Fig. 1; Supplementary Data SD8). In the quantification analysis of the MOP environmental similarity, a greater environmental similarity was observed under the MPI scenario relative to the calibration scenario for both northern and southern matrilineages. As for reptiles, the southern matrilineages showed little environmental similarity (Supplementary Data SD8).

Under warm (current) weather conditions, niche models for northern matrilineages of reptiles matched the known distribution (Fig. 1; see Supplementary Data SD7). For southern matrilineages of reptiles (U. stansburiana, U. nigricaudus, and C. ruber), the models under current warm conditions suggest areas of climatic suitability to the south of the peninsula and east of their respective distributions (Fig. 1; Supplementary Data SD7).

Transfers of the ENMs of reptiles showed large areas of the BCP where northern matrilineages were able to find climatically suitable areas under the cold-weather scenario (LGM; Fig. 1; see Supplementary Data SD7). The transfers of southern matrilineages suggest that climatic suitability for these matrilineages receded toward the southern part of the peninsula. The three transfer scenarios showed the same expansion pattern of areas of climatic suitability for the realized niche of northern matrilineages along with a reduction of areas of climatic suitability for the realized niche of southern matrilineages (Fig. 1; Supplementary Data SD7).

The MOP analysis showed that of the three GCMs tested for northern reptiles, MPI was the scenario with the greatest environmental similarity, that is, showing equivalent climates between calibration and transfer areas, while MIROC was the least similar scenario (Supplementary Data SD8). For the southern matrilineages of reptiles, the three GCMs showed large areas of low environmental similarity between the warm (calibration) and cold (transfer) scenarios, particularly in the northern part of the peninsula (Supplementary Data SD8).

In general, the MOP analysis (Supplementary Data SD8) contrasting the environmental characteristics between calibration (warm-weather conditions) and transfer (cold-weather conditions) areas yielded some areas with non-analogous climatic conditions between both areas. The GCM with the highest number of pixels where extrapolation is required, evidencing non-analogous conditions between calibration and transfer areas, was MIROC. In contrast, the GCM with the lowest number of pixels with extrapolation in transfers was the MPI, which included six models with some pixels showing extrapolation in the central part of BCP (Supplementary Data SD7 and SD8). It should be highlighted that those pixels showing extrapolation were interpreted with caution in our discussion. However, the pattern derived from our analyses was found for the three GCMs, where the northern matrilineages possibly find favorable weather conditions in other areas besides their current distribution. Thus, our interpretation is unaffected regardless of whether the transfers of GCM–MIROC are considered or not.

Quantification of ecological niches between matrilineages.

The BSTs indicate that the environments where each matrilineage is found are not more similar than expected by chance (P < 0.05; Table 2). That is, each matrilineage is found in particular environments that differ from those of its counterpart.

Table 2.

Comparison of niches in the ecological space: Niche Overlap Test (NOT) and Niche Divergence Test (NDT) and geographic space: Background Similarity Test (BST) of northern and southern matrilineages of vertebrates inhabiting the Baja California Peninsula.

	BST		NOT		NDT
Species (northern–southern matrilineages)	D	I	D	I	D	I
Mammals
Ammospermophilus leucurus	0.397	0.677	0.048	0.063	0.192a	0.295a
Otospermophilus beecheyi	0.490a	0.806a	0.025	0.098a	0.039a	0.066a
Dipodomys merriami	0.410	0.694	0.150	0.188	0.129a	0.211a
Chaetodipus fallax	0.207b	0.516b	0.121	0.241	0.141a	0.245a
Thomomys nigricans	0.450	0.749	0.107	0.140	0.155a	0.234a
Lepus californicus	0.298	0.623	0.019	0.029	0.021a	0.026a
Birds
Campylorhynchus brunneicapillus	0.592c	0.862c	0.076	0.144	0.132a	0.237a
Auriparus flaviceps	0.459	0.775	0.028	0.047	0.042a	0.066a
Reptiles
Urosaurus nigricaudus	0.355	0.668	0.010	0.016	0.022a	0.035a
Uta stansburiana	0.337	0.591	0.119	0.132	0.132a	0.200a
Crotalus ruber	0.341	0.613	0.032	0.041	0.046a	0.059a

P < 0.05. D = Schoener’s index; I = modified Hellinger distance. BST: H₀ = similarity; NOT: H₀ = equivalent occupied niches; NDT: H₀ = fundamental niches are not divergent.

^a  H  ₀ not rejected.

^bLower values for D and I. ^cHigher values for D and I.

The bird C. brunneicapillus obtained the highest values for D and I indices (with values ranging from 0 to 1, where 1 represents identical ENMs); this suggests that the differentiation level between its matrilineages is the lowest of the whole set of species studied. Conversely, the mammal C. fallax obtained the lowest values of D and I, suggesting a greater differentiation between the environments that each matrilineage inhabits, compared to the other species (Table 2). According to the NOTs and NDTs, the current niches of the northern and southern matrilineages are not equivalent. However, these results do not support the divergent evolution of niches but reflect different environments to which each matrilineage has access (Table 2).

Relationship between physical and ecological conditions in the mid-peninsular region.

The Mantel test between physical resistance distance and ecological conductance distance yielded differential results between matrilineages. No statistical significance was observed for the northern matrilineages (Table 3), suggesting no correlation between the measures of physical landscape resistance and ecological conductance for the northern matrilineages. This finding supports the idea that altitude in the area of contact between matrilineages does not affect the spatial expression of the ecological niche of northern matrilineages. As regards the correlation tests between resistance distance and ecological conductance for the southern matrilineages, H₀ (no correlation between variables), was rejected in six (54%) tests (Table 3); this result suggests for this matrilineages that landscape conditions (altitude) have a correlation with the spatial expression of the ecological niche of southern matrilineages.

Table 3.

Mantel test to evaluate the correlations between resistance distance and ecological conductance measures of species with genetic breaks in the mid portion of the Baja California Peninsula.

Species	Northern		Southern
	RM	P value	RM	P value
Mammals
Ammospermophilus leucurus	0.150	0.193	0.230	0.148
Otospermophilus beecheyi	0.198	0.288	−0.041	0.524
Dipodomys merriami	0.014	0.498	0.287	0.020
Chaetodipus fallax	−0.137	0.997	−0.110	0.981
Thomomys nigricans	−0.137	0.704	0.221	0.102
Lepus californicus	−0.059	0.464	0.336	0.038
Birds
Campylorhynchus brunneicapillus	0.387	0.069	0.548	0.032
Auriparus flaviceps	−0.105	0.709	0.789	0.001
Reptiles
Urosaurus nigricaudus	−0.083	0.546	0.131	0.111
Uta stansburiana	0.020	0.408	0.350	0.012
Crotalus ruber	0.338	0.039	0.739	0.001

H  ₀ = no correlation between variables; RM = Mantel test statistic (r). Bold indicates significant value (P ≤ 0.05) and H₀ rejected.

Discussion

Climatic variation between matrilineages.

The results from the comparison of the environments where each lineage is found suggest that individual matrilineages have particular requirements in terms of the realized niche. The variables analyzed indicate that the greatest environmental disparity between lineages involves temperature (Table 1). Temperature appears to be more important for reptiles than endotherms (Huey 1982). The maximum difference was observed in the minimum temperatures for the northern and southern lineages of U. stansburiana (11.8°C; Table 1). This species displays different morphotypes according to the climatic conditions where it thrives (Sinervo et al. 2000; Zamudio and Sinervo 2000; Sinervo and Zamudio 2001).

The north and south regions of the BCP show distinct rainfall patterns: winter rains in the north (71–82%) versus summer rains in the south (75–90%; García and Mosiño 1968; Salinas-Zavala et al. 1998). These differences in precipitation pattern are evident in the phenology of plants, which mainly bloom in spring in the north (Delgadillo 1998) and mainly late summer in the south (León de la Luz et al. 1996). The difference in bloom is in relation to the rainy season and has an effect on seed availability: seeds are an important food source for many vertebrate species.

Whorley and Kenagy (2007) noted that the antelope squirrel (A. leucurus) shows different reproductive cycles throughout its range associated with seasonality of precipitation and ecosystem composition and diversity. Although no studies are available on the reproductive cycle of the antelope squirrel in the BCP, differences in rainfall patterns between the north and south of the peninsula may affect its reproductive cycles. Differences in reproductive timing also could limit gene flow.

The marked climatic differentiation also affects the distribution range of mammal species, as in the case of D. merriami (Álvarez-Castañeda et al. 2009), T. nigricans (Álvarez-Castañeda and Patton 2004; Trujano-Álvarez and Álvarez-Castañeda 2007), A. leucurus (Whorley et al. 2004), C. rudinoris and C. arenarius (Riddle et al. 2000a), and other species (Goldman and Moore 1945; Álvarez and LaChica 1974; Wiggins 1980). Based on on our results (Table 1), the range of precipitation in the distributional area of northern matrilineages is lower than in the south.

Ecological niche of matrilineages.

Our results: ENMs in warm-climate scenario, transfers to cold-weather conditions, and comparisons of climatic suitability (Fig. 1; Table 2; Supplementary Data SD7) support the idea that there are differences in the realized niches of matrilineages of species with genetic divergence in the middle of the BCP. This study shows that the current niches are not equivalent; however, this may be due to the fact that each population has access to different environments. When species encounter either different ecological conditions or abrupt changes in climate, individuals can adapt to local environments, which may result in subpopulations with slightly different tolerance ranges (Grismer and Greene 2002; Labra et al. 2009; Soberón and Peterson 2011). Some of the species analyzed show adjacent geographic distributions (U. nigricaudus, C. ruber, and C. brunneicapillus), while others show overlapped distributions (U. stansburiana, A. leucurus, D. merriami, and L. californicus). However, the ecological conditions in the areas where they live are dissimilar (Table 2).

It has been shown that environmental variations can lead to genetic exchange among individuals of a population inhabiting similar environments and reduce the exchange among populations living in different environments (Wang 2013; Sexton et al. 2014). Frequently, there is no close correspondence between geographic and ecological spaces, implying that species in geographic proximity are not necessarily ecologically similar (Soberón et al. 2017). Thus, phylogenetically affine species with nearby distributions also may show ecological divergence (Blair et al. 2013).

In past transfers (cold-weather scenario), the disparity of realized niches between the northern and southern matrilineages is more noticeable. Transfers to the LGM indicated that northern matrilineages were able to find areas of climatic suitability in much of the peninsula, while models showed only small areas with suitable climate conditions for southern matrilineages. During the LGM, climate conditions were colder in the peninsula (Ferusquía-Villafranca and Torres-Roldán 1980). These climate conditions favored northern matrilineages, as these were associated with colder temperatures than southern matrilineages (Table 1).

The southern BCP has been considered a Pleistocene refuge for invertebrates (González-Trujillo et al. 2016; Valdivia-Carrillo et al. 2017). In particular, the eastern slope of the southern mountain ranges of the peninsula (“La Giganta” and “La Laguna”) may have functioned as a barrier preventing the passage of wind from the Pacific Ocean, thus contributing to maintaining warm conditions on the eastern slope (González-Medrano 2012), and thereby serving as a refuge for southern matrilineages.

The Global Circulation Models (GCMs) used for transferring niche models to cold-weather conditions showed divergences. MIROC was the GCM with the lowest ecological similarity between current warm and past cold conditions. In contrast, the MPI model showed the highest environmental similarity between the two scenarios contemplated. These results are similar to those reported by Graham et al. (2014), where transfers under the MIROC scenario were more dramatic versus the CCSM scenario. These results are due to the fact that MIROC simulates warmer and more humid conditions (Otto-Bliesner et al. 2006).

Relationship between physical and ecological conditions in the mid-peninsular region.

Authors of phylogeographic studies rightly argue that species with limited dispersal capabilities are more likely to show phylogeographic breaks caused by geographic barriers relative to species with greater dispersal capabilities (e.g., Avise 1994; Bond et al. 2001; Irwin 2002). However, some studies show that genetic isolation also can be explained by ecological data (Irwin 2002). The strength of a phylogeographic break caused by a barrier preventing gene flow tends to decrease rapidly over time once this barrier disappears (Irwin 2002). The phylogeographic break persists for a long time only for species with very limited dispersal capabilities. These low dispersal conditions are likely to produce phylogeographic breaks unrelated to a geographic barrier (Irwin 2002).

The mid-peninsular seaway, which a number of authors have used to explain genetic divergence in multiple species, is hypothesized to have closed about a million years ago (Upton and Murphy 1997; Holt et al. 2000; Lindell et al. 2005; Munguía-Vega 2011). This time has been sufficient for northern and southern matrilineages of species to come in contact and for the genetic break to disappear. However, this has not happened for many vertebrates of mainly low vagility (Riddle et al. 2000b; Álvarez-Castañeda et al. 2008, 2009; Trujano-Álvarez and Álvarez-Castañeda 2013). Other mammals (C. arenarius, C. rudinoris, C. spinatus, Canis latrans, Neotoma bryanti, Peromyscus fraterculus, P. maniculatus, Urocyon cinereorgenteus, Vulpes macrotis; STA-C, personal observation) have been studied and do not show genetic divergence. The Mantel test showed no relationship between resistance distance and ecological conductance measures for northern matrilineages and the existence of a relationship for six southern matrilineages. These findings suggest that in the area of contact between matrilineages, southern matrilineages show more resistance to movement throughout the area of contact than do northern matrilineages.

To recapitulate, this study found that matrilineages of species within which are found discrete breaks of genetic divergence have ecological niches that are associated with distinct climatic conditions for distinct populations; these ecological distinctions would act as an ecological barrier. Northern matrilineages are associated with the winter rainfall pattern and lower temperatures, contrasting with southern matrilineages, which are associated with the summer rainfall pattern and higher temperatures. It was evident that under different global climatic conditions (warm versus cold), the spatial expression of ecological niches remains distinct.

Currently, altitude in the area of contact between matrilineages does not constrain reconnection between matrilineages. Thus, the association with different rainy and flowering seasons suggests the possibility of differential resource availability, likely leading to diachronic reproduction. Reproductive isolation and ecological barriers may be the primary mechanisms leading to genetic divergence between vertebrate matrilineages in the BCP.

This pattern also can occur in other regions of the planet with geographic characteristics similar to that of the peninsula (i.e., a nearly north-south alignment). It therefore is advisable to consider the characteristics of the landscape and analyze the ecological conditions of the species in areas where genetic breaks are found in spite of the apparent absence of physical barriers.

Supplementary Material

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1.—Geographic location of the Baja California Peninsula, climate group (from Köopen 1918), and genetic divergences in vertebrate groups.

Supplementary Data SD2.—Geographic coordinates (in decimal degrees) for each matrilineage of vertebrate species used to generate ecological niche models.

Supplementary Data SD3.—Calibration area for ecological niche models and spatial presence records for north and south matrilineages of vertebrates with genetic breaks in the Baja California Peninsula. Dark gray, northern matrilineages; light gray, southern matrilineages.

Supplementary Data SD4.—Set of variables and parameterization selected for each matrilineages used for calibration of ecological niche models in Maxent. Northern matrilineages (N); southern matrilineages (S); regularization multiplier (RM); area under the curve (AUC); omission rate (OR); Akaike information criterion (AIC).

Supplementary Data SD5.—Details of the relationship between resistance distance and ecological conductance in the mid-peninsular region for vertebrates with genetic divergence in the Baja California Peninsula.

Supplementary Data SD6.—Boxplot diagram of temperature and precipitation values for northern and southern matrilineages of vertebrate species with genetic breaks in the central part of the Baja California Peninsula. Key to species: Ammospermophilus leucurus (Ammos), Otospermophilus beecheyi (Otos), Dipodomys merriami (Dipo), Chaetodipus fallax (Chaet), Thomomys nigricans (Thom), Lepus californicus (Lepu), Campylorhynchus brunneicapillus (Campy), Auriparus flaviceps (Auri), Urosaurus nigricaudus (Uro), Uta stansburiana (Uta), and Crotalus ruber (Crot). The key is followed by the matrilineage to which it belongs: northern (N) and southern (S).

Supplementary Data SD7.—Ecological niche models for vertebrate species (reptiles, birds and mammals) and their respective transfers to the cold-weather scenario (LGM, ~21,000 years ago), with the three Global Circulation Models used (CCSM, MIROC, MPI). Light gray represents models for northern matrilineages. Dark gray represents models for southern matrilineages. Middle gray represents the overlap between northern and southern models.

Supplementary Data SD8.—Analysis of environmental similarity through “Mobility-Oriented Parity” between calibration and transfer areas of ecological niche models for vertebrates with intra-peninsular genetic breaks. L = Low environmental similarity. H = High environmental similarity.

Acknowledgments

We thank the Collection of Mammals of the Centro de Investigaciones de Baja California Sur (CIBNOR) for the availability of mammal records of the BCP. Thanks also to Sean P. Maher, Pedro Peña, Andrés Lira-Noriega, José Anadón, and Townsend Peterson for their valuable comments on the manuscript. María Elena Sánchez-Salazar edited the English manuscript. LC-S is a recipient of a student fellowship from the Consejo Nacional de Ciencia y Tecnología of Mexico (CONACYT 337421).