Abstract
We used ecological niche modelling to study the relative roles of climate and interspecific interactions in defining the parapatric contact between closely related species (Crotalus mitchellii and C. stephensi) and to predict refugia during the last glacial maximum. The modelled suitable habitat for C. stephensi correctly predicts the existing parapatric border between it and C. mitchellii, suggesting that C. stephensi's range at the border is limited by climatic factors. In contrast, the suitable habitat for C. mitchellii does not correctly predict the existing parapatric boundary; rather the suitable habitat of this species extends into the range of C. stephensi, suggesting the latter species, not climatic factors, limit the range of C. mitchellii. Modelled refugia of C. stephensi are much smaller than modern suitable habitat and are partially situated at the current parapatric border, whereas the refugia of C. mitchellii are similar to its current suitable habitat, though also shifted to the south. Ecological niche modelling appears to be a useful tool for studying the interplay between climate and competition in determining boundaries between parapatric species. It also appears to be useful for predicting past suitable habitats of species, because predicted refugia are congruent with independent estimates from molecular phylogeography.
1. Introduction
Species’ geographic ranges are determined by environmental, intraspecific, and interspecific interactions [1]. Climate and competition are two important factors limiting species’ ranges [2–4], where climate is a limiting factor to a species’ fundamental niche and competition is a limiting factor to a species’ realized niche. The concept of the fundamental niche was formalized by Hutchinson [5] as an “-dimensional hypervolume” involving ecological variables that allow a species to exist and persist. The limits of the fundamental niche of a species may be realized if a species is a dominant competitor; otherwise, a species might occupy its realized niche, a subset of the fundamental niche ([5], but see [6]).
One approach to exploring and describing a species’ niche is by ecological niche modelling (ENM). ENM combines species’ occurrence data and environmental variables to predict geographic areas that may provide suitable habitat for a species. This approach generally assumes a Gleasonian model for species’ response to climate change in that species have individualized tolerances and responses to climate fluctuations [7]. Graham et al. [8] provided support for this model of species’ response based on changes in fossil mammal assemblages during the Late Quaternary climate fluctuations. In a theoretical approach to studying species’ ranges, Holt [9] emphasized that species’ traits associated with environmental tolerances evolve by natural selection and species’ ranges may also move, contract, or expand independent of climate fluctuations.
Many software tools have been developed to facilitate data and algorithm processing for ENM (e.g., ANUCLIM, DIVA-GIS, GARP, GRASP, MAXENT, openModeller Desktop, etc.). The software packages generally implement different algorithms for predicting suitable habitat by associations between occurrences and environmental variables, termed the “correlative approach” [10]. To project the calculated suitable habitat onto a map, digital mapping programs are available (e.g., ArcView, DIVA-GIS, Idrisi, etc.). Suitable habitat is determined using the correlative approach by finding geographic areas that are environmentally similar to areas in which organisms are located. Factors affecting a species’ range may be recognized by exploring patterns produced by predicted suitable habitat of one or more species [10]. For example, Sánchez-Cordero et al. [4] used ENM to create maps of potential suitable habitat of individual species within a felid community. They calculated the overlapping pattern of the suitable habitat and identified competition as an important factor in limiting the southern border of the range of one species in that assemblage.
Comparing patterns of multiple species’ suitable habitats may be instructive for evaluating biotic interactions influencing the species’ ranges; further, comparing suitable habitat for a single species over time is instructive to describe how abiotic factors might influence the species’ range size and location. Waltari et al. [11] used ENM and a correlative approach to predict Pleistocene refugia for species at the Last Glacial Maximum (LGM) by projecting the suitable habitat of a species onto a climate reconstruction at the LGM. They compared the predicted suitable habitat of species from ENM with predicted refugia of species from molecular phylogeographic studies. They found that many of the species they compared have high geographical overlap and significant spatial correlations between the two predictions.
The rattlesnakes Crotalus mitchellii and C. stephensi are sister species inhabiting southwestern North America and maintain a parapatric border. The parapatric border occurs along a zone less than a few kilometers wide approximately cutting the south central Mohave Desert on a southeast/northwest transect [12]. These species are presumed to be ecologically similar and were, until recently, recognized as conspecific subspecies, C. m. pyrrhus and C. m. stephensi [13, 14]. It has been demonstrated by morphological differentiation [12] and molecular divergence [15, 16] that the two populations are in fact separate species. Crotalus mitchellii continues to be recognized with subspecies (C. m. pyrrhus occurs at the parapatric border) whereas C. stephensi is now recognized as a monotypic taxon. Using molecular phylogeography, Douglas et al. [15] proposed three allopatric Pleistocene refugia. A refugium was identified near Death Valley for C. stephensi in the northern part of its current range. Two refugia were identified for C. m. pyrrhus to the east and west of the apex of the Sea of Cortez. This phylogeographic study provides a framework to reference ENM predictions of Pleistocene refugia for C. m. pyrrhus and C. stephensi.
We used ENM to study the relative roles of climate and competition in defining parapatric contact between closely related species, and to predict the suitable habitat of these species during past climatic conditions, in this case at the LGM. The use of ENM to address the parapatric boundary between these taxa is useful because actual ecological data for populations at the contact zone are nonexistent.
2. Methods
2.1. Species’ Occurrence Data
We collected and georeferenced species’ occurrence data for Crotalus mitchellii pyrrhus () and C. stephensi () in the Mohave Desert, USA. Data on occurrence come from both textual descriptions of site localities and actual GIS coordinates associated with vouchered specimens from multiple museum collections (see Acknowledgments). We examined all specimens used in this study to verify identifications. This study is focused within the Mohave Desert to explore the parapatric contact between species; therefore, the occurrence data include instances of presence only for the Mohave Desert. The occurrences include all of the species’ range of C. stephensi, because it occurs mostly in the Mohave Desert. However, we only include part of the range of C. m. pyrrhus, because it occurs outside of the Mohave Desert. When multiple species’ occurrences were obtained from one site, only one datum was counted as a presence at that GIS coordinate.
2.2. Ecological Niche Modelling
We used 19 bioclimatic variables from Worldclim, which includes means and variations in temperature and precipitation derived from interpolations of monthly global climate layers with high spatial resolution [17]. The bioclimatic variables were transformed to be more biologically meaningful than raw monthly temperature and precipitation data and encompass different variable types (i.e., yearly trends, seasonality, and extreme environmental conditions). Our analyses were preformed within the North American region, and bioclimatic variables are 2.5′ spatial resolution.
Two algorithms were used in a consensus approach to conservatively develop ENMs [18]: Maxent [19] and GARP [20]. We used the desktop version of Maxent (version 2.3) [19] and open Modeller Desktop [21] to construct ENMs using the Maxent and GARP algorithms, respectively. Both methods use only presence data in the model, because validity of absence data is usually questionable [19, 20]. Maxent calculates a percent contribution for each of the bioclimatic variables used to build species’ suitable habitat models. The bioclimatic variables and their percent contributions were compared between the species qualitatively. Statistical comparisons are not appropriate for these highly correlated bioclimatic variables. All input parameters follow methods in Waltari et al. [11]; although they used a different program to implement the GARP algorithm, Desktop GARP (version 1.1.6). An ENM was obtained from each algorithm and projected on current climate maps, resulting in two suitable habitat models.
We trained models with 80% of the occurrences and tested models with the remaining 20%. We used the area under the receiver operating characteristic curve (AUC) for model validation [22]. LGM climate layers for two climate models, Community Climate System Model (CCSM) [23] and the Model for Interdisciplinary Research on Climate (MIROC) [24], were downloaded from PMIP2 [25]. For projections of species’ suitable habitat at the LGM, a suitable habitat model was obtained from each algorithm and projected on a climate model reconstruction, resulting in four potential refugia reconstructions for each species at the LGM. We constructed suitable habitat model by setting thresholds on output data from the two programs running the ENM algorithms and constructed a model consensus following the conservative approach detailed in Waltari et al. [11]. Suitable habitat maps were made using a combination of DIVA-GIS and ArcGIS Desktop 8.1.
We do not account for potential evolutionary change in species’ suitable habitat from the LGM to modern times. Although these species have been shown to evolve morphologically over short evolutionary time scales [26], their suitable habitat has a much slower rate of evolutionary change [27]. Over the past three glacial-interglacial cycles, the rate at which rattlesnakes have shifted their suitable habitats in response to climate change was two to three orders of magnitude more than the rate at which suitable habitat shifted due to evolutionary change [27]. Therefore, we considered the potential evolutionary change influencing the shift in suitable habitat between the LGM and the modern for C. m. pyrrhus and C. stephensi to be negligible.
We were unable to quantitatively test the refugia hypotheses of the molecular phylogeographic study from Douglas et al. [15] and our study, but we made descriptive inferences about the two hypotheses. The proposed refugia based on the molecular phylogeography were given in a textual format [15]. While we were able to locate the approximate geographic locality, we were unable to infer the approximate size/area of the geographic locality, hindering any possible quantitative analysis of these hypotheses.
3. Results
All ENM model results have a high AUC, validating the models. The AUC for Crotalus stephensi and C. m. pyrrhus is 0.997 and 0.990, respectively. The percent that individual bioclimatic variables contributed to building the suitable habitat model of a species differed between Crotalus stephensi and C. m. pyrrhus (Table 1). Precipitation of the warmest quarter is the most important bioclimatic variable for these desert species, with 28.8% contribution for C. stephensi and 28.1% contribution for C. m. pyrrhus. Mean diurnal range, isothermality, and mean temperature of the wettest quarter are important variables for C. stephensi and not C. m. pyrrhus, with greater than 10% contribution discrepancy between the two species. On the other hand, temperature seasonality and mean temperature of the coldest quarter are important variables for C. m. pyrrhus with a greater than 10% contribution discrepancy between the two species.
The suitable habitat predicted by Maxent was relatively smaller for both species than GARP for all climate projections (data not shown). ENMs were combined for further analyses by defining the suitable habitat of a species as suitable habitat predicted by both algorithms or projected on both climate models [11]. The modern suitable habitat is larger in area for C. stephensi than for C. m. pyrrhus (Figure 1). The refugia for both species at the LGM are similar in size, with C. m. pyrrhus having slightly larger and more continuous potential refugia (Figure 2). The suitable habitat of C. stephensi dramatically increases from glacial to interglacial time. The suitable habitat of C. m. pyrrhus slightly increases from glacial to interglacial time. The modern suitable habitat for C. m. pyrrhus and C. stephensi definitively overlaps; C. stephensi has a suitable habitat with a boundary that approximately predicts the parapatric contact that currently exists between these two species. C. m. pyrrhus fails to occupy potentially suitable habitable within the current range of C. stephensi.
(a)
(b)
(a)
(b)
4. Discussion
Interspecific competition is probably an important factor in maintaining a parapatric contact between Crotalus stephensi and C. mitchellii pyrrhus. An overlap in their suitable habitat indicates geographic space in which both species may have suitable habitat (Figure 1); likewise, a better competitor is identified by the exclusion of one species from occupying all of its suitable habitat, resulting in a parapatric boundary. In this context, C. stephensi is likely the better competitor, effectively excluding C. m. pyrrhus from occupying its suitable habitat. It is not likely that intraspecific dispersal of C. m. pyrrhus limits its range at the parapatric contact zone. Douglas et al. [15] provided evidence to support recent eastward population expansion by C. m. pyrrhus, which implies this species is a good disperser. On the contrary, C. stephensi inhabited the Mohave Desert region since at least the Mid-Pliocene. Based on these observations, one would expect C. stephensi to be the better competitor, having had more time to adapt to the climate of the Mohave Desert than C. m. pyrrhus. Therefore, with respect to the parapatric contact between these species, it seems that C. stephensi is limited by climate conditions and C. m. pyrrhus is limited by C. stephensi.
The modern suitable habitat of C. stephensi extends far north above its observed range within the Mohave Desert (Figure 1(a)). The Mohave Desert is a rain shadow desert, and species well adapted to this environment might be expected to cue on bioclimatic variables associated with the rain shadow effect. The suitable habitat of C. stephensi seems to follow along the east side of the Sierra Nevada and the Cascade Range where the rain shadow is somewhat continuous. Other biologic or geologic factors, besides the biotic and abiotic factors studied here, may exclude C. stephensi from occupying the northern extent of its suitable habitat, such as, geologic barriers, competitors, or intraspecific migration/dispersal rates. Within the Mohave Desert, C. stephensi occupies almost its entire suitable habitat.
The modern suitable habitat of C. m. pyrrhus extends northwestward and southeastward beyond its observed occurrences within the Mojave Desert (Figure 1(b)). The northwest extent of the suitable habitat pushes well into the observed range of C. stephensi and covers approximately half of the observed range of C. stephensi. This part of the suitable habitat is not actually occupied by C. m. pyrrhus. The southeast extent of the suitable habitat reaches to areas around Phoenix, Arizona. C. m. pyrrhus does occur in this region (though not as far east as predicted by the suitable habitat model) and extends southwestward into Baja California. Douglas et al. [15] found two molecular clades within C. m. pyrrhus partitioned into an eastern and western clade. C. m. pyrrhus in the Mohave Desert belongs to the eastern clade, and the suitable habitat for C. m. pyrrhus, from occurrences only within the Mohave Desert, covers the entire range of this eastern clade. The suite of bioclimatic variables that the eastern clade of C. m. pyrrhus responds to are likely different from the ones that the western clade responds to; otherwise, we may have found the suitable habitat of C. m. pyrrhus extending south into Baja California and west into coastal southwestern California.
Refugia predicted by ENM for C. stephensi at the LGM are greatly reduced in area compared to the modern suitable habitat (Figure 2(a)), again suggesting the range of this species is highly dependent on environmental factors. Refugia are scattered throughout and also somewhat to the north of the modern range of C. stephensi and stretch far in length along the sides of the modern range of C. m. pyrrhus. There also seems to be a significant area of refugia to the west of the Sierra Nevada. Refugia predicted for C. m. pyrrhus at the LGM only slightly decrease in area compared to its modern suitable habitat (Figure 2(b)), again suggesting that the range of this species was not dramatically affected by extreme changes in climate conditions between the LGM and now. Douglas et al. [15] proposed refugia for C. stephensi near Death Valley and refugia for C. m. pyrrhus at two allopatric localities to the east and west of the apex of the Sea of Cortez, accounting for the east and west molecular clades within C. m. pyrrhus. From ENM, predictions of refugia for C. stephensi encompass Death Valley and predictions of refugia for C. m. pyrrhus do occur on the east side of the apex of the Sea of Cortez. These refugia predictions are consistent with the hypotheses of Douglas et al. [15].
Peterson et al. [28] found conserved ecological niches between sister taxon pairs in birds, mammals, and butterflies. They suggested that vicariant speciation events do not induce niche divergence and over longer geologic time scales niche evolution proceeds slowly. In contrast, the sister taxon pair we studied have greatly diverged in environmental variables that characterize the niche since their hypothesized vicariant speciation event. A molecular clock estimate dates the split between C. m. pyrrhus and C. stephensi to approximately 3.7 Ma [15]. The development of the Bouse Embayment [29], which spread west through the Imperial Valley [30] and reached the Transverse Range during the Mid-Pliocene [31], is hypothesized as the vicariant geological event that separated C. m. pyrrhus and C. stephensi [15]. One species, C. m. pyrrhus, was considered by Douglas et al. [15] to be a generalist due to recent range expansion and large geographic area that it currently occupies. The other species, C. stephensi, was considered a specialist because most of its range is restricted to the Mohave and southern parts of the Great Basin deserts. Subsequent to these species’ divergence event, evolution of the niche proceeded to generate species-specific habitat requirements (Table 1) as well as differential species-specific factors (competition and climate) limiting their distributions within the Mohave Desert.
In conclusion, using ENM seems to be useful for studying the interactions between climate and competition in shaping parapatric contact between two species. In this study, the existing boundary of one species was predicted by ENM, but ENM also suggests the boundary of the second species is limited by the first: Crotalus stephensi appears to exclude C. mitchellii pyrrhus from suitable habitat through competition. ENM is also useful for predicting refugia of species and, supporting molecular phylogeographic hypotheses, it provides compatible independent estimates of refugia at the LGM.
Acknowledgments
The authors thank Rob P. Guralnick for assistance with North American paleoclimate layers. They are grateful to the following institutions: BYU, CAS, MBM, MSB, MVZ, SDNHM, SNHM, UAZ, UTEP, UTA, and UT, for providing access to specimens and space for examination. They acknowledge the Laboratoire des Sciences du Climat et de l’Environnement (http://pmip2.lsce.ipsl.fr/) for collecting and archiving data. The PMIP2 Data Archive is supported by CEA, CNRS, MOTIF, and the Programme National d’Etude de la Dynamique du Climat.