Changes in the Geographic Distribution of the Diana Fritillary (Speyeria diana: Nymphalidae) under Forecasted Predictions of Climate Change

Climate change is predicted to alter the geographic distribution of a wide variety of taxa, including butterfly species. Research has focused primarily on high latitude species in North America, with no known studies examining responses of taxa in the southeastern United States. The Diana fritillary (Speyeria diana) has experienced a recent range retraction in that region, disappearing from lowland sites and now persisting in two phylogenetically distinct high elevation populations. These findings are consistent with the predicted effects of a warming climate on numerous taxa, including other butterfly species in North America and Europe. We used ecological niche modeling to predict future changes to the distribution of S. diana under several climate models. To evaluate how climate change might influence the geographic distribution of this butterfly, we developed ecological niche models using Maxent. We used two global circulation models, the community climate system model (CCSM) and the model for interdisciplinary research on climate (MIROC), under low and high emissions scenarios to predict the future distribution of S. diana. Models were evaluated using the receiver operating characteristics area under curve (AUC) test and the true skill statistics (TSS) (mean AUC = 0.91 ± 0.0028 SE, TSS = 0.87 ± 0.0032 SE for representative concentration pathway (RCP) = 4.5; and mean AUC = 0.87 ± 0.0031 SE, TSS = 0.84 ± 0.0032 SE for RCP = 8.5), which both indicate that the models we produced were significantly better than random (0.5). The four modeled climate scenarios resulted in an average loss of 91% of suitable habitat for S. diana by 2050. Populations in the southern Appalachian Mountains were predicted to suffer the most severe fragmentation and reduction in suitable habitat, threatening an important source of genetic diversity for the species. The geographic and genetic isolation of populations in the west suggest that those populations are equally as vulnerable to decline in the future, warranting ongoing conservation of those populations as well. Our results suggest that the Diana fritillary is under threat of decline by 2050 across its entire distribution from climate change, and is likely to be negatively affected by other human-induced factors as well.


Introduction
Understanding how species distributions might shift with the changing climate is a critical component of managing and protecting future biodiversity. Hundreds of species in the United States and elsewhere have responded to the warming climate by shifting to higher latitudes or elevations [1][2][3][4]. Such range shifts have been documented in a number of taxa [5][6][7], including alpine plants [8], marine invertebrates [9], marine fish [10], mosquitoes [11], birds [12,13], and butterflies [1,[14][15][16][17][18]. A number of species distribution models have been developed to predict the impacts of climate change on species distributions, including bioclimate envelope models, which are useful first estimates of the potential effects of climate change on altering species' ranges [19]. Bioclimate envelope models work by identifying the climatic bounds within which a species currently occurs, and then delineating how those climatic bounds will shift under various future climate projections [20][21][22][23].
Most often, researchers are limited to presence-only occurrence data, requiring the use of indirect methods to infer a species' climatic requirements [8,24,25]. One of the best performing models using presence-only data is maximum entropy modeling, or Maxent [26], which performs well even with low sample sizes typical of rare species [19,27,28]. Maxent works by comparing climate data from occurrence sites with those from a random sample of sites from the larger landscape to minimize the relative entropy of statistical models' fit to each data set. Species distribution models such as Maxent have been criticized for being overly simplistic, because they do not incorporate external biotic factors such as species interactions [20,27,29]. However, such bioclimate envelope models have been used to project with reasonable accuracy whether species ranges will increase or decrease under a changing climate [19,[30][31][32], which was the primary objective of this study.
Speyeria diana (Nymphalidae) (Cramer 1777) is a butterfly species endemic to the southeastern United States and is currently threatened across portions of its range. This species is of particular conservation interest because it has experienced a range collapse in recent decades resulting in an 800-km geographic and genetic disjunction between western populations in the Ouachita and Ozark Mountains and populations in the southern Appalachian Mountains, and has shifted to a higher elevation at an estimated rate of 18 m per decade [33]. This range contraction is consistent with the predicted effects of a warming climate, and might represent the first such documented case in the southeastern United States, though the region has experienced other environmental changes in recent decades as well [33]. Previous research using coalescent-based population divergence models dated the earliest splitting of the western population from the east at least 20,000 years ago, during the last glacial maximum [34]. In addition, recent geometric morphometric evidence from the wings of S. diana further support this long-term spatial and genetic isolation [35]. In light of these pieces of evidence, we used Maxent to model the future distribution of S. diana under several future climatic scenarios, in order to forecast how the range of the butterfly might shift under predicted conditions. Forecasts of large range reductions (over 50%), or small overlaps between current and future ranges (less than 50%), would suggest high vulnerability to climate change. Range reductions of any size in the western distribution would likely threaten those populations that are genetically isolated and adapted to relatively low dispersal, with the negative effects of genetic drift [34,35].

Study Species
The Diana fritillary, Speyeria diana, is a large and sexually dimorphic nymphalid butterfly, endemic to the southeastern United States. Adult males emerge in late May to early June, with females flying several weeks to a month later [36]. Once mated, each female can lay thousands of eggs singly on ground litter during the months of August and September in the vicinity of Viola spp., the larval host plant for all Speyeria [37]. After hatching, first instar larvae immediately burrow deep into the leaf litter layer of the forest floor, where they overwinter [38]. In spring, larvae feed on the foliage of freshly emerging violets. Adult Diana butterflies are often found along forest edges or dirt roads containing tall, conspicuous nectar sources such as milkweeds, butterfly bushes, or other large summer and fall composites [39][40][41][42]. While males begin to die off in late July, females may persist in large numbers, although somewhat cryptically, through October [42].

Distributional Dataset
We searched for all known records of S. diana, from publications, catalogued and uncatalogued specimens in public and private collections in the United States and Europe, online databases, contemporary field surveys by scientists and amateurs, and our own field surveys. We obtained distributional data from 1323 pinned S. diana specimens from 33 natural history museum collections in the United States and Europe (Table 1). Four hundred thirty-five additional records  were provided by the Butterfly and Moth Information Network and the participants who contribute to its BAMONA project. Our literature survey produced 153 records (1818-2011) across 54 U.S. counties (Table 2). We also collected 469 S. diana butterflies in our own field surveys (Table 3). Our dataset essentially represents a complete dataset of all publicly available records for the species, and is as comprehensive as for any taxon in the region [33]. For this reason, our dataset should be especially informative in creating an accurate bioclimate envelope for the species, as collection bias is a major consideration with ecological niche modeling [43,44].

Species Distributional Modeling
We developed species distribution models using the popular machine-learning algorithm for ecological modeling, Maxent [26]. Maxent estimates a species' probability distribution that has maximum entropy (closest to uniform), subject to a set of constraints based on the sampling of presence-only data [45]. Because of the difficulty and impracticality of obtaining accurate absence data, presence-only data are most often used in species distribution modeling. In order to offset the lack of absence data, Maxent uses a background sample to compare the distribution of presence data along environmental gradients with the distribution of background points randomly drawn from the study area [46][47][48]. Locality data and the randomly sampled background points are combined with climatic data to predict the probability of the species' occurrence within each raster grid cell. We used environmental climate data from WorldClim [49] at 30 arc-second resolution or approximately 1 km 2 grid cells. Bioclimate variables and elevation layers were each clipped to the extent of North America using ESRI (Environmental Systems Research Institute) ArcMap 10.0, and data extracted to S. diana sample localities. Additionally, we collected the same types of locality data for three other species of North American butterflies (Speyeria cybele, Speyeria idalia, Battus philenor), which served as 5628 random background points for our models. We utilized these background data to minimize spatial bias in our modeling, as data represented by similar butterfly species can be used as pseudo-absence data with the same collection bias as our occurrence data, improving the accuracy of the model [50,51].
Climatic variables included 19 derived bioclimatic variables that describe annual and seasonal variation in temperature and precipitation, as well as elevation, averaged for 1950-2000 (Table 4). One concern when modeling species distributions is the strong correlation that occurs between multiple climate variables, which can significantly influence model predictions of species distributions [52]. To test for co-linearity, we performed spatial autocorrelation statistics between all pairs of the 19 bioclimate variables using ESRI ArcMap 10.0. We then selected the most biologically meaningful variable for each group of two or more variables with Pearson correlation coefficients higher than 0.7 (Table 4). This allowed us to reduce the number of bioclimate variables to the nine potentially most important ones, which were: Minimum Temperature of Coldest Month, Mean Temperature of Driest Quarter, Precipitation of Wettest Month, Precipitation of Driest Month, Precipitation of Driest Quarter, Isothermality, Mean Diurnal Range (Mean of monthly (maximum temperature-minimum temperature)), Temperature Annual Range, and Annual Precipitation, along with elevation (Table 4). These variables are typically considered to be important determinants of butterfly distributions, as they relate to life history traits. Butterflies are highly sensitive to weather and climate, particularly changes in temperature and rainfall [53]. For example, mean temperature of the coldest month is related to the overwintering survival of first instar larvae, growing degree days above 5 • C are regarded as a surrogate for the developmental threshold of the larvae, water balance corresponds to the moisture availability for the larval host and adult nectar plants, and the mean temperature of late summer ensures proper adult emergence and mating [54][55][56][57][58][59]. Temperature changes affect all aspects of butterfly life history, from their distribution and abundance [14,54], to their realized fecundity [60,61]. Changes in rainfall levels can influence butterfly larvae indirectly through changes in host plant quality, and generally rainfall is considered to be beneficial because it enhances host plant growth [62].
One concern when modeling species distributions is whether the occurrence records are spatially biased with respect to site accessibility (e.g., towns, roads, trails) [63]. To address this concern, we applied a spatial filter to remove all sampling points that were within 5 km of each other using ESRI ArcMap 10.0. The spatial filter resulted in 254 unique presence points for S. diana that were used in the final model. We first modeled the distribution of these 254 occurrences in present-day climate, and then projected the fitted species distribution under two future climate scenarios for the period 2040-2069 (hereafter referred to as 2050). Future climate scenarios were taken from two global circulation models (GCMs) obtained from www.worldclim.org; the community climate system model (CCSM) [64] and the model for interdisciplinary research on climate (MIROC) [65,66]. These GCMs differ in the reconstruction of several climatic variables and are well known to produce different outcomes for butterfly species [67,68]. For example, in hind-casting Mediterranean butterflies, the CCSM model projects narrower distributions at the last glacial maximum than does MIROC [65,66]. For each of these two GCMs, we considered two different representative concentration pathways (RCPs) [69][70][71][72][73], which are cumulative measures of human emissions of greenhouse gases from all sources expressed in Watts per square meter. These pathways were developed for the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [67] and correspond to a total anthropogenic radiative forcing of RCP = 4.5 W/m 2 (low) and RCP = 8.5 W/m 2 (high) [72,73].
We used Maxent's default parameters [26,50] and a ten-fold cross-validation approach to further reduce bias with respect to locality data. This method divides presence data into ten equal partitions, with nine used to train the model, and the tenth used to test it. These partitions generate ten maps (one map per run), with each raster grid cell containing a value representing the probability of occurrence. These values were used to designate habitat suitability ranging from 0 (unsuitable habitat) to 1 (highly suitable habitat) ( Figure 1). We averaged the resulting maps for the current climate, and for the two GCMs under RCP = 4.5 and RCP = 8.5. This method resulted in the production of a "low" and "high" average prediction for S. diana species distribution in 2050, represented with habitat suitability maps. We measured the goodness of fit for the models using the area under the curve (AUC) of a receiver-operating characteristic (ROC) plot [74]. We used criteria of Swets [75] and considered AUC values higher than 0.7 representative of model predictions significantly better than random values of 0.5 or less [26,27,74]. Because AUC has been recognized as a somewhat questionable measure of accuracy, especially when used with background data instead of true absences [74,76], we also calculated the TSS (true skill statistics), a threshold-dependent evaluation metric [76,77]. The relative importance of each variable's contribution was assessed by sequential variable removal by Jackknife [26].

Results
Species distributional modeling resulted in "excellent" model fits for Speyeria diana, with a mean AUC = 0.91 ± 0.0028 SE, TSS = 0.87 ± 0.0032 SE for RCP = 4.5; and a mean AUC = 0.87 ± 0.0031 SE, TSS = 0.84 ± 0.0032 SE for RCP = 8.5 (Table 1)  Both climate models indicate that the loss of core distributional area is modest, with an average of 91.3% of present distributional areas retained. The most drastic reduction in habitat is apparent across the southern Appalachian Mountains (Figure 2).

Results
Species distributional modeling resulted in "excellent" model fits for Speyeria diana, with a mean AUC = 0.91 ± 0.0028 SE, TSS = 0.87 ± 0.0032 SE for RCP = 4.5; and a mean AUC = 0.87 ± 0.0031 SE, TSS = 0.84 ± 0.0032 SE for RCP = 8.5 (Table 1)  Both climate models indicate that the loss of core distributional area is modest, with an average of 91.3% of present distributional areas retained. The most drastic reduction in habitat is apparent across the southern Appalachian Mountains (Figure 2).

Discussion
Our ecological niche models predicted that the amount of suitable habitat for Speyeria diana will decline substantially by the year 2050 across its entire distribution. Both CCSM and MIROC climate models predicted severe habitat loss and fragmentation in the southern Appalachian Mountains by 2050, with some range expansion predicted into higher latitudes in both eastern and western populations. High elevation habitat will be an important refuge for the species across the entire distribution, as the range of S. diana is already shifting to higher elevations at an estimated rate of 18

Discussion
Our ecological niche models predicted that the amount of suitable habitat for Speyeria diana will decline substantially by the year 2050 across its entire distribution. Both CCSM and MIROC climate models predicted severe habitat loss and fragmentation in the southern Appalachian Mountains by 2050, with some range expansion predicted into higher latitudes in both eastern and western populations. High elevation habitat will be an important refuge for the species across the entire distribution, as the range of S. diana is already shifting to higher elevations at an estimated rate of 18 m per decade [33]. Recent evidence further suggests that some S. diana populations may already be adapting to high elevations, as S. diana female forewings from high elevation populations were found to be narrower than low elevation populations, indicating that these females may be more mobile than those from low elevations with wider forewings [35].
Unlike populations in the eastern distribution, the wing shape of western populations of S. diana appears to be better adapted for lower dispersal, which is in alignment with findings that western populations of S. diana are both spatially and genetically isolated [35]. Our models predicted that the southern edge of the highly suitable habitat in the west will recede by 2050; However, as was found in the southern Appalachian Mountains, the suitable habitat was predicted to expand in the higher elevations of the Ozark and Ouachita mountains of Arkansas. The genetic isolation of western populations may ultimately prevent them from adapting to higher elevations as successfully as populations in the eastern distribution of the species. If this is the case, lower elevation populations will be even more vulnerable to climate change than our models predict.
We would like to note that all ecological niche models should be used and interpreted with caution because of various sources of bias and error that result in inaccurate predictions [78]. Some have questioned the applicability of bioclimatic modeling at regional scales because of the somewhat coarse resolution [79]. However, we are confident that the size of our study area, and our uniquely extensive dataset, provide sufficient data to forecast climate-driven range shifts in S. diana with accuracy. Both global circulation models (CCCM and MIROC) were very closely aligned in their outcomes, indicating strong agreement between them. Climate is well understood to play a primary role in shaping the distributions of species [80], and we are confident in our overall findings that the suitable habitat for S. diana will decline and become increasingly fragmented by 2050.

Conclusions
These results highlight the importance of maintaining connectivity of the suitable habitat for S. diana, especially in the eastern populations that appear most vulnerable to increased fragmentation and loss of suitable habitat. These populations in the eastern distribution of S. diana harbor important genetic diversity that may become lost through genetic drift if these populations become small and isolated. The Ozark and Ouachita Mountains of Arkansas and Missouri appear to be least vulnerable to loss of suitable habitat from climate change, and therefore will be important for the future conservation of S. diana after 2050. As a result of the geographic and genetic isolation of the western populations, conservation of suitable habitat in the west is equally as important as in the east. Our climate models show that the 800-km disjunction across the center of the range of S. diana is not due to complete absence of suitable habitat, but more probably a result of the extensive habitat fragmentation regionally across the Ohio River Valley from agricultural land use change, and other human related factors that were not included in our models. We conclude that maintaining well-connected low and high elevation habitats across the entire distribution of S. diana, both now and into the future, will be necessary for this species, even under conservative forecasts of climate change.