Challenges and Opportunities for Terrapene carolina carolina Under Di ﬀ erent Climate Scenarios

: An unprecedented rate of global climate change as a result of human impacts has a ﬀ ected both endotherms and ectotherms. This is of special concern for ectotherms, such as reptiles, as these species are su ﬀ ering from large population declines and lack the dispersal ability of other taxa. There are many protected areas across the United States; however, these areas are fragmented, which hinders dispersal. We examined species distribution and dispersal capabilities for Terrapene carolina carolina , a relatively narrow range, low dispersal, and vulnerable species. We created climatic suitability models to predict changes in suitable habitat and identiﬁed important predictor variables. We modeled three time periods using MaxEnt and hypothesized that there would be an increase in northern habitat. We found that most of the suitable habitat changed at the northern end of the range and that mean temperature of driest quarter had the most inﬂuence on future predictions. Overall there were relatively moderate changes in suitable habitat, but where these changes occur a ﬀ ects accessibility. As an example, we examined these local scale movements within Oak Openings Region and found that individuals are capable of dispersing to new suitable habitats; however, other physical barriers will hinder movements. In conclusion, there is a critical need to protect this vulnerable reptilian species and our results suggest that T. c. carolina will expand their distribution northward. We suggest that land managers increase connectivity among protected areas to facilitate dispersal, but future studies should incorporate other dynamic ecological factors at ﬁner spatial scale.


Introduction
Global biodiversity is under major threat by a multitude of factors, such as nutrient loading, carbon dioxide enrichment, and invasive species; similar to these other threats, land-use and climate change pose considerable challenges [1]. Climate change distresses both endotherms and ectotherms by affecting various population aspects including density, distribution shifts, phenology, morphology, and genetics [2]. Many species are unable to adapt fast enough to the rapid rate of rising global average temperature. Although the impacts on endotherms is concerning, there have been relatively few studies examining the decline of less charismatic, but highly vulnerable and climatically sensitive ectothermic species [3][4][5][6][7]. This is especially concerning for reptiles, such as turtles, because these species lack the dispersal ability of other taxa such as birds, mammals, and many fishes [8] and population declines are expected to be expedited [5]. Additionally, ectotherm physiology plays a strong role where it is influenced by thermal conditions, which in turn are often highly affected by the local scale physical characteristics of a habitat [9]. In a changing climate, thermal heterogeneity is crucial for providing both sunny, warm patches and shady, cool patches to move in or out of as needed. However, maintaining optimal thermal physiology can be difficult in thermally variable environments, We hypothesized that there would be an increase in suitable habitat in the northern portion of the range as a result of increased warming and precipitation. Therefore, we predicted that there would be a northward shift in distribution. We predicted that habitat suitability would increase over time with greater gains for higher greenhouse gas concentration scenarios. We hypothesized that mean temperature of warmest quarter and precipitation of driest quarter would drive T. c. carolina responses as a result of increased vulnerability to very hot and dry conditions.

CSMs (Regional Scale)
CSMs were created for the range of T. c. carolina by combining (a) species occurrence data with (b) bioclimatic layers to identify past, current, and future distributions of habitat using MaxEnt. MaxEnt applies maximum entropy modeling, a machine-learning technique, which includes a set of environmental layers with georeferenced occurrence localities. Then the model creates a probability distribution for the species where each grid cell contains a predicted suitability condition. Climate suitability was modeled at different time periods using six bioclimate variable layers. We used presence-only data in MaxEnt version 3.3.3k [59,60] to build CSMs for each scenario for the T. c. carolina distribution range. Presence-only records were randomly subsampled and separated withholding 30% of the records for test data and using 70% as training data, with ten replicates.
(a) Occurrence. Species occurrence data for T. c. carolina were obtained from the Global Biodiversity Information Facility (GBIF) [61] and a current radio telemetry study [62]. We downloaded 2579 records (decimal degree) from GBIF across this species range in the United States that were gathered as presence-only records from 1970 to 2019. We used a subset of presence records from the current radio telemetry study in Oak Openings Region. These records were condensed by applying a 75 m apart rule to reduce spatial autocorrelation. The methods for tracking individuals were carried out in accordance to approved guidelines and permitted by Ohio Department of Natural Resources (20-016), Bowling Green State University's Institutional Animal Care and Use Committee (1001429-9), and with permission from Metroparks Toledo.
(b) Bioclimatic layers. Bioclimatic variables represent extremes of limiting environmental factors. Nineteen bioclimatic layers were downloaded from CHELSA (climatologies at high resolution for the Earth's land surface areas) [63,64] for three conditions: LGM (22,000 years ago), current , and future 2050. The future scenarios were taken from the Coupled Model Intercomparison Project Phase 5 Model Community Climate System Model version 4 for four scenarios (Representative Concentration Pathway (RCP) 2.6, RCP 4.5, RCP 6.0, and RCP 8.5). All layers were downloaded in 30 arc-second resolution and clipped to Eastern United States region using the Environmental Systems Research Institute (ESRI) ArcCatalog version 10.2.2 [65]. Climatic layers for the current model were examined using band collection statistics tool in ESRI ArcGIS version 10.2.2 [66] to identify highly correlated variables, i.e., r > 0.80 (Table S1). We selected variables to minimize correlations; the final set included bioclimatic (BIO) 8 (mean temperature of wettest quarter); BIO9 (mean temperature of driest quarter); BIO10 (mean temperature of warmest quarter); BIO16 (precipitation of wettest quarter); BIO17 (precipitation of driest quarter); and BIO18 (precipitation of warmest quarter) layers (Table 1). Each quarter has a temporal resolution of three months of the year.

CSMs Regional Evaluation
We examined the influence of each environmental predictor variable using percent contribution and permutation importance for each climate scenario model. Percent contribution values depend on the pathway taken by the MaxEnt code for the optimal solution and these values may vary depending on the algorithm taken, while permutation importance depends on the final model and is derived from the contribution of each variable from random permutation among the training points [67]. The response curves for the average model were examined to explore how the environmental variable affected the MaxEnt model prediction and the jack-knife results for the average model were used to identify the significance of individual environmental variables.
We evaluated our CSMs by building alternative climatic distribution maps for comparison. These models were created from a bioclimatic profile (i.e., the range of values at each occurrence point) of current conditions (Table 2). We reclassified the climatic variables as suitable (range from minimum to maximum) and unsuitable (values below the minimum and above the maximum) for LGM, current and future scenarios. Using raster calculator, we added each reclassified binary variable together and created a hotspot distribution map. Low overlap ranged from 0-2 layers, medium overlap ranged from 3-4 layers, and high overlap ranged from 5-6 layers. These maps represent suitable climate conditions based on current thermal usage by T. c. carolina.

Occurrence Model (Local Scale)
To place the regional climate results in a local context, we explored a case study within Oak Openings Region in northwestern Ohio. We developed five separate occurrence models for T. c. carolina in northwestern Ohio by combining (a) species occurrence data with (b) bioclimatic data, and (c) local environmental data using MaxEnt. Presence-only records were randomly subsampled and separated withholding 30% of the records for test data and using 70% as training data, with ten replicates. We created five comparison distribution models for Current and Future scenarios. Each model included all six bioclimatic variables, percentage of forest, percentage of early successional, and normalized difference vegetation index (NDVI). Only the six bioclimatic variables varied across models as each set were based on current and the four future emission scenarios.
(a) Occurrence. Species occurrence data for T. c. carolina were obtained from the same subset of presence records from a current radio telemetry study [62]. We used 47 records (decimal degree) within Oak Openings Region that were gathered as presence-only records from 2016 to 2019. These records were at least 75 m apart to reduce spatial autocorrelation. As stated earlier, the methods for tracking individuals were approved and permitted.
(b) Bioclimatic layers. We used the same six selected bioclimatic variables from our regional scale models for current and future conditions. All layers were clipped to Oak Openings Region using ESRI ArcGIS version 10.2.2 [66].
(c) Local environment data. We created two additional environmental variables, i.e., land cover, and NDVI. We ran a supervised land cover classification using training data. We acquired one image from 19 June 2016 Landsat 8 for Path 20, Row 31, and we clipped the scene to only contain our entire study area [68]. Training sites were delineated in ENVI 5.5 [69] across the study site for four classes, i.e., forests, early successional, human-modified, and water. We used a cropland mask [70] and then ran the supervised land cover classification. We merged the cropland layer with the classification map and created a final 5-class land cover map for Oak Openings Region. We converted the categorical land cover map to continuous data in FRAGSTATS version 4.2.1 [71]. To identify the percentage of landscape for forest and early successional land cover, we ran a moving window analysis with a 120 m buffer. A NDVI was created using the 19 June, 2016, Landsat 8 image for our study area in ArcGIS version 10.2.2 [66].
We examined model accuracy using the area under receiver-operating-characteristic curve (AUC) values and converted each average model into a binary suitable and unsuitable distribution map using the logistic threshold for the maximum test sensitivity plus specificity (MSS). Additionally, we examined the influence of each environmental predictor variable using percent contribution and permutation importance for each climate scenario model.

Connectivity Model (Local Scale)
Using a connectivity model, we examined whether individuals would be capable of dispersal along with changing climatic conditions. We used fine-scale movement data from an ongoing tracking study [62] to create three dispersal distance scenarios (minimum: 10 m/d; average: 35 m/d; and maximum: 200 m/d) to determine if T. c. carolina would be able to feasibly migrate within their life time (conservative estimate of 40 years). At this northern portion of the distribution range, we assumed that individuals would be able to travel for 252 days or 36 weeks per year as their active period occurs from mid-March to mid-November, with no movement during overwintering. Individuals in the southern portion of the distribution range would have a longer activity period and we expect they would be able to travel greater distances. Additionally, we would expect variation in movement with changing temperatures; however, we present a conservative estimate with the three dispersal scenarios.
Climate conditions are not the only barriers that T. c. carolina face when attempting to disperse, other physical barriers such as roads and streams reduce landscape connectivity. To estimate whether individuals could disperse to predicted novel areas of suitable climate, we built dispersal cost layers for (a) occupancy models, and five physical factors including (b) roads and streams, (c) elevation, and (d) land cover. We combined each cost layer together for four scenarios: Occupancy Model 1 and five physical factors; Occupancy Model 2 and five physical factors; Occupancy Model 3 and five physical factors; and five physical factors only, then we calculated the least cost path for each scenario for comparison. We created cost raster layers in ESRI ArcGIS 10.2.2 [66] for each occupancy model and each of the five physical layers [6]. To determine actual cost of travel for each physical layer, we used our field data to inform our decision.
(a) Occupancy models. We utilized the three occupancy models that were built in Section 2.3. We created a separate cost layer for each occupancy model in ESRI ArcGIS 10.2.2 [66]. Each average continuous probability map was converted to a binary suitable and unsuitable map using the MSS threshold. We then assigned a cost of 1 to values above and a cost of 2 for values below the MSS threshold.
(b) Road and streams. We acquired road and stream shapefiles from U.S. Census Bureau [72] and clipped the layers to our study area. We then calculated the Euclidean distance for road and stream layers to create a distance from continuous layer. We created an environmental profile layer using the presence-absence occurrence points, like our bioclimatic profile (Table 2), and used the minimum to maximum range values to assign a cost of 1, and a cost of 2 to any values above or below that range.
(c) Elevation. We obtained a digital elevation image courtesy of U.S. Geological Survey Earth Explorer server [73] and we derived slope from this digital elevation model. We created an environmental profile layer and used the minimum to maximum range values to assign a cost of 1, and a cost of 2 to any values above or below that range.
(d) Land cover. We utilized an existing land cover map specifically created for the study area [74]. We estimated how much each land cover would cost based on T. c. carolina habitat use. Upland prairie, upland savanna, sand barrens, upland deciduous forest, and upland coniferous forest were assigned the smallest travel cost of 1, followed by wet prairie, swamp forest, floodplain forest, and wet shrubland a cost of 2, then turf/pasture, Eurasian meadow, and cropland a cost of 3, then residential/mixed, and dense urban a cost of 4, and, finally, perennial pond the highest cost of 5.
Our cost raster layers were simplified for practicality with such large-scale data. Fine-scale models could include greater detail for each cost raster if warranted. Each raster was reclassified into a binary cost raster, except land cover which had 5 cost values, and then we calculated the sum for all the values per pixel to create a final comprehensive dispersal cost raster for the least cost path analysis. We used the cost path tool in ArcGIS version 10.2.2 [66] to derive the path of least resistance from the starting point to the destination point. It uses a cost surface raster, i.e., difficulty of travel, a cost distance analysis, i.e., inputs a source or starting point and calculates the accumulative cost to travel to each cell from the source, and a backlink raster, i.e., the direction from each cell to its lowest cost neighbor. The least cost pathway is calculated by backtracking from the destination point using the directions created by the backlink raster and the total cost accumulated from the distance raster to create a pathway to the source. Here, we estimated the best single pathway for each scenario from the starting point in Oak Openings Preserve to the destination point in Wildwood Metroparks, a straight-line distance of 20.6 km. We calculated the length of the least cost paths in km to determine how far T. c. carolina would need to travel at a minimum.

CSMs (Regional Scale)
At a large scale, climate changes are likely to have relatively modest impacts on suitable habitat in the future compared to changes in the past. For each scenario the average model (i.e., average of 10 replicates) performed better than random with AUC values greater than 0.8 (ranging from 0.808 to 0.823) ( Table 3). Using the MSS threshold value obtained for each scenario, suitable habitat ranged from 24.5%-30.3% across scenarios (Table 3). Habitat suitability continuous probability maps varied among each time period (Figure 1). All models except 2050 RCP 2.6 were largely affected by the environmental variable: mean temperature of the driest quarter based on the percent contribution (Table 4). In contrast, 2050 RCP 2.6 distributions were largely affected by mean temperature of the warmest quarter using the percent contribution (Table 4). We found that suitable habitat decreased over time, 5.8% from LGM to current conditions (Table 3). We found predicted increases in suitable habitat from current to future predictions: 2050 RCP 2.6 (4.6%), 2050 RCP 6.0 (4.4%), 2050 RCP 4.5 (4.7%) and 2050 RCP 8.5 (3.7%) ( Table 3).  We examined the most influential variable using response curves ( Figure S1) and our test data jack-knife results for all our models. We found that for LGM, precipitation of driest quarter had the highest gain when used in isolation and mean temperature of warmest quarter decreased the gain the most when omitted. Additionally, for all other models, we found that mean temperature of driest quarter had the highest gain when used in isolation and precipitation of driest quarter decreased the gain the most when omitted.
Over time, 45%-54% of the habitat remained unchanged in suitability across each scenario. Total change varied from 19% to 29% for T. c. carolina under the different climate scenarios (Figure 2, Table 5). The greatest habitat gain occurred between current conditions and 2050 RCP 4.5 with an increase of 39% in suitable habitat; while the greatest loss of suitable habitat occurred between LGM and current, with a decrease of 42%. All four future scenarios ranged from 10% to 17% loss of suitable climatic habitat and ranged from 36% to 39% gain of suitable climatic habitat. Overall, suitable habitat was lost from LGM to current, while future scenarios gained suitable habitat.  Our alternative CSMs supported our hypothesis that suitable habitat for this species would have a northward expansion. This held true for changes from LGM to current model, while suitable habitat both expanded northward and contracted in the south for future models ( Figure S2). The predicted suitable habitat losses across low to high categories ranged from 4.3% to 31.7% across all scenarios, while predicted suitable habitat gains ranged from 0.0% to 31.4% ( Table 6). The greatest habitat gain for future change occurred for 2050 RCP 8.5 medium category with an increase of 31.4% in suitable habitat, while the greatest loss of suitable habitat for future change occurred 2050 RCP 8.5 high category with a decrease of 31.7% (Table 6). Table 6. Predicted changes in the percent of climatically suitable habitat for T. c. carolina using hotspot maps based on the number of environmental predictor layers that overlap with one another. We classified low overlap (0-2 layers), medium overlap (3)(4), and high overlap (5-6). Abbreviations: Last Glacial Maximum (LGM) and Representative Concentration Pathway (RCP).

# of Layer Overlap
LGM

Occurrence Model (Local Scale)
At the local scale, changes in suitability were readily detectable as a result of climate changes. We found that Model 1, 78% of the habitat was unsuitable and 22% was suitable (Table 7, Figure 3a). The model was largely affected by percentage of forest, precipitation of driest quarter, and precipitation of wettest quarter based on the percent contribution (Table 8). We found that the probability of presence increased as percentage of forest and precipitation of driest quarter increased ( Figure S3). We examined our jack-knife results and found that percentage of forest had the highest gain when used in isolation and precipitation of wettest quarter decreased the gain the most when omitted. In contrast, all future models had less suitable habitat than Model 1, ranging from 14.4% to 18.1% (Table 7, Figure 3b-e), and except for Model 2, all future models were largely affected by percentage of forest based on percent contribution and as this variable increases so does probability of presence (Table 8). Model 2 was largely affected by mean temperature of warmest quarter, followed by percentage of forest, and mean temperature of wettest quarter (Table 8). Additionally, our jack-knife results showed that mean temperature of warmest quarter had the highest gain when used in isolation and precipitation of driest quarter decreased the gain the most when omitted. Both Model 3, Model 4, and Model 5, respectively, were secondly largely affected by mean temperature of driest quarter, mean temperature of warmest quarter, and mean temperature of wettest quarter. Our jack-knife results varied across each future scenario; Model 3 had the highest gain for mean temperature of driest quarter, while Model 4 and Model 5 it was mean temperature of warmest quarter. Conversely, the variables that decreased the gain the most when omitted for Model 3 and Model 4, were mean temperature of driest quarter, while for Model 5 it was mean temperature of wettest quarter.

Connectivity Model (Local Scale)
While suitable habitat may be available both now and in the future, the question is whether it is likely to be accessible. We estimated that T. c. carolina would be able to disperse 101 km, 353 km, and 2016 km, respectively, per dispersal distance scenario (i.e., minimum, average, maximum) within their lifetime. These dispersal distances were based on feasible movement distances and do not consider physical barriers that would hinder movement. However, we examined constraints on movements in our local study area. We found that the least cost pathway from two protected areas varied across the five scenarios. All models had similar required distances to travel, ranging from 22.95 km to 23.84 km. Both Model 1 and Model 5 had similar pathways, while Model's 2 through 4 were more similar (Figure 4). We found that near the upper portion of the pathways, all models converge for the same least cost path. Individuals would have to cross 44 to 49 roads and avoid the Toledo Express Airport that lies in the center between the two parks to get from point A to point B. Even with the minimum dispersal distance, individuals can feasibly migrate despite a multitude of environmental barriers.

Discussion
MaxEnt modeling is a common method for modeling future changes in habitat suitability. Our prediction that this species would have expanding climatic habitat northward was supported with both our MaxEnt and alternative distribution models. Other studies have shown northern expansion for several taxa, such as Ambrosia artemisiifolia [75], Lepidoterans [76], Dolichophis caspius [6], and several mesopredators: Lontra canadensis, Mephitis, Canis latrans [77]. For our study, these changes are based on temperature and precipitation, which are important to reptilian movements. Overall, with warming climate, reptiles will face similar problems such as overheating and will require greater fine-scale heterogeneity which includes more shaded areas provided by more vegetation [78]. Forest canopy cover can provide critical structural complexity in its underlying microclimates [79] and facilitate fine-scale thermoregulation; however, local movements are limited by other physical barriers such as streams, roads, and human-modified landscapes. Therefore, habitat fragmentation may lead to the isolation of dispersing individuals within habitats that become suboptimal in terms of temperature and precipitation [80]. Utilizing land surface temperature maps can help identify critical areas of thermal limits for different species and land managers can manage the habitat to increase or decrease temperatures (e.g., manage amount of canopy cover). As with many problems that species face, there will be winners and losers as in Popescu et al. (2013) [81], which found that reptiles show mixed responses to climate change across species. Our study suggests that T. c. carolina may benefit from climate change at a landscape scale with increased suitable habitat in the future, assuming that habitat is detectable and accessible, which will be dependent on local scale characteristics.
Our results for LGM, current and future regional models supported mean temperature of driest quarter as the most influential bioclimatic variable driving the model, except for RCP 2.6 scenario where mean temperature of warmest quarter was the most influential. This likely occurred because precipitation patterns may change more drastically at higher temperatures, i.e., RCP 2.6 to RCP 8.5 scenarios. Many climate models for North America have shown an increase in precipitation over time [82]. This may occur from warmer winter and spring temperatures that lead to more snowmelt and rain-on-snow events that create severe flooding [83]. We found that greatest change occurred between LGM and current conditions; however, the temporal difference between these two scenarios is much greater than between current and future scenarios. Therefore, we caution making too many conclusions about the change from past to present. Global climate changes are occurring at an unprecedented rate and many of these changes are alterations in temperature and precipitation. Surprisingly, we saw little changes among the four future scenarios. Although our results suggest that T. c. carolina will have the greatest response to these changes within the driest and warmest quarters of the year, as we predicted based on physiological constraints. We found small percent contributions from precipitation of the wettest and warmest quarters (< 5.7%) across all models. This suggests high vulnerability during the driest and warmest temperatures of the year and declines may be expedited if these quarters become longer or temporally shift. Surprisingly at the local scale, we found an opposite trend for the current model where precipitation of wettest and driest quarters mattered the most and mean temperature was less important. This suggests that currently, temperature is vital for regional suitability, while precipitation matters locally. However, when we examined future models, mean temperature was most important, but varied which quarter mattered the most. At the local scale, T. c. carolina utilize a wide range of habitats for thermoregulation; however, when their body temperatures are too high, they seek cover in flooded or wet areas [84], spending days in these areas especially during high temperatures and droughts [85]. Therefore, precipitation and where water is located during the driest quarter would result in greater movement displacement. Physiologically, this makes sense because water stress exacerbates temperature regulation issues, therefore shaded and floodplain areas will become disproportionately more important during the driest season for these vulnerable species. Other studies have utilized different sets of these bioclimatic predictor variables and found that other variables were most important. For example, the most influential variable for D. caspius was mean temperature of coldest quarter [6]; for Pseudopus apodus it was temperature seasonality and annual precipitation [86]. For several coral snake species, they were most affected by annual mean temperature, precipitation of wettest month, and precipitation of warmest quarter [87]. Since each study used a different set of bioclimatic factors, it is difficult to effectively compare variable importance. However, each species will be affected differently based on their ecological requirements, physiology, and location of their range. This underscores the need to translate regional climate change effects to the local scale context.
Under current conditions, our MaxEnt model predicted that T. c. carolina has a concentrated distribution within the middle of the eastern United States which is consistent with known range distributions. We expected that all future scenarios at a regional scale would result in the gain of suitable habitat and our results supported our hypothesis. We found that the proportion of suitable or unsuitable habitat showed moderate changes, i.e., 19% to 27%, it is critical to note, though, that where these changes occur will have large impacts on this species. Most of the changes, both gains and losses, occur on the distribution edge and relatively few changes within the center of the distribution. Although identifying these large-scale changes can help pinpoint critical priority areas for conservation efforts, we need finer-scale data to provide recommendations for on the ground management. In the future, suitable habitat may become inaccessible and isolated, therefore these models should distinguish between suitable and accessible locations [88]. Identifying these areas, requires downscaling from regional to local scale maps. Therefore, we compared local scale models across current and future scenarios. We found that suitable habitat decreases over time and is limited (i.e.,~15%). For our case study, percentage of forest was the driving factor for all models except Model 2 which suggests that the availability of shade will drive individual responses. Depending on which emission scenario is used will vary which quarter (i.e., wettest, driest, or warmest) will drive turtle responses. Once again highlighting the need for local scale climatic datasets, especially when modeling future suitable habitat. However, in all cases, suitable habitat is limited and fragmented, which makes dispersal challenging.
Warming climate may make forest habitat more suitable; however, it will also make open areas less suitable and costly to cross. This suggests that canopy cover is important for dispersal, although NDVI is a good model for canopy cover, it cannot substitute for measuring canopy cover at a fine scale. Our results, though, support the need to incorporate land cover or other remote sensed data with climate factors, as percentage of habitat and NDVI contribute to habitat suitability. We recognize that projecting changes into the future is challenging and deriving models with more models requires a lot of assumptions. Therefore, a mixture of different techniques that address both changes in spatial and temporal scale are needed to tackle this challenge. From a climatic perspective, T. c. carolina may not need to move, especially at the northern part of the distribution, therefore other biological processes (i.e., increased disturbance, fragmentation, competition, predation) will drive local scale movements.
We examined this problem using both our CSMs and finer-scale environmental data in Oak Openings Region to evaluate whether individuals could feasibly disperse despite climatic changes. Based on field data an individual could feasibly disperse to a new area within its lifetime to accommodate a changing landscape. However, this assumes that there are relatively few to no barriers and constant movement each day towards potentially new suitable habitat, as well as an ability to detect these new habitats. For example, in Oak Openings Region, we found that individuals would need to travel a total of 22.95 to 23.84 km to get from one protected area to another (an additional 2 to 3 km from the straight-line distance) using only physical factors on the landscape. It would take individuals approximately 5.5 months to 9.5 years based on our minimum to maximum dispersal scenarios and seasonal activity period to travel this distance. Incorporating climate suitability produced a northern route and different scenarios did not vary much. We did find that the end portion of all models were the same and we suggest that local managers focus on protecting this route for potential dispersal. Although this distance ranges from 23% to 24% of the minimum dispersal distance scenario, it is unlikely that most of the individuals will take the least costly pathway, always following a straight-line path. Individuals would have to carefully traverse a multitude of suboptimal habitat including roads, rivers, crops, and developed areas. Additionally, individuals would have to traverse a variety of thermal gradients that may become too intolerable in the future. Individuals would be able to travel the required distance; however, when considering other structural features, many individuals would likely perish trying to disperse such large distances. T. c. carolina can persist in developed areas and can traverse croplands and flooded areas; however, they are susceptible to road mortality. Deaths on roads often change the population demographic and cause skewed adult sex ratios [89]. Additionally, this species is more susceptible to road mortality than other turtle species because they tend to close their shell when threatened and remain closed for longer periods of time [90]. Increased road traffic will increase the probability of mortality and assuming no further habitat changes, will lead to a threshold effect that causes widespread local extinction [91]. Although we found that there will be more suitable habitat in the future, it may not do any good with continued habitat fragmentation. We found that this species will face many difficult barriers and potentially have increased mortality when dispersing to future suitable habitat even with expanding suitable habitat, it may not be accessible.
Although T. c. carolina may not be as heavily impacted by climatic changes in terms of suitable area at the landscape-scale, it is important to provide a model for local scale context. T. c. carolina will remain under threat from a variety of sources, e.g., habitat destruction, invasive species, environmental pollution, disease, unsustainable use [5]. Therefore, models should incorporate other biologically relevant variables (e.g., roads, land cover, elevation, streams/lakes, etc.) to represent the synergistic effects of multiple threats. Such a multifaceted approach is critical because the response of local populations will be dependent on regional weather patterns and local structural characteristics. In addition, climate change interacts with the landscape by altering the configuration and composition of land cover. For an example, lakes are sensitive indicators of climate alterations such as fluctuating water levels [92] and timing changes in ice formation and thawing [93]. As with many other studies, land-use patterns are often a larger driving factor that influences populations more so than climate alone [94,95]. Our climatic habitat models illustrate that climatic change may be beneficial for this species; however, there are local scale challenges when other factors are examined that will affect how the species responds to these changes.

Conclusions
There is a critical need to protect and manage a flagship reptilian species such as T. c. carolina. We developed CSMs for this vulnerable species and evaluated the climatic variables that are influential in determining the geographic distribution and we predicted changes in potentially suitable habitat between different time periods. Our results predict moderate suitable habitat expansion. Additionally, we have identified areas that are susceptible to loss or gains in habitat within T. c. carolina distribution that can be monitored or managed. One of the most important aspects of our approach is to consider climate change impacts across a range of temporal and spatial scales from regional climate models to local occupancy models. Understanding the local context of large-scale changes is necessary to effectively manage for potential impacts on native flora and fauna and can highlight critical priorities. In addition to identifying where suitable habitat is and will be, we suggest that land managers should work towards increasing connectivity among suitable areas to facilitate dispersal. However, to examine these local scales, we need finer scale climatic data in combination with a variety of remotely sensed and locally collected data to inform our models. Finally, we suggest that future studies incorporate other dynamic ecological factors that influence distribution shifts.