Potential change in the distribution of an abundant and wide-ranging forest salamander in a context of climate change

Climate change already affects species in many ecosystems worldwide. Since climate is an important component of a species’ ecological niche, up-to-date information about climatic niches is needed to model future species distributions in a context of climate change. The eastern red-backed salamander (Plethodon cinereus) is a wide-ranging woodland species and one of the most abundant vertebrates in northeastern North America. Though salamanders contribute to several forest ecosystem functions, little is known about their climatic niche and future distribution. Using a dataset of 400,090 observations from 8302 localities in 5 Canadian provinces and 22 American states, we determined the current climatic niche of P. cinereus and predicted how the species’ distribution could shift in a context of climate change, especially in the northern part of its range. We also aimed to document factors that could affect the species’ distribution. We show that P. cinereus can live in various geographic and climatic conditions and tolerate a wide range of seasonal temperatures. The species’ current potential and future (until 2061– 2080) distributions show a gap of up to 400 km with the northern limit of its current observed distribution. Assuming a mean colonization rate of approximately 100 m per year, we calculated that P. cinereus would need about 4000 years to reach the northern limit of the future distribution range modeled for the 2061–2080 period. The climate-modeled future distribution suggests that the presence of P. cinereus could decrease in the south and increase in the north. This, combined with the potential presence of habitats that are unsuitable for the species’ colonization in the north and with interspecific interactions in the south, could induce a contraction of the species’ range. Regardless of climate warming, the physical environment and natural and anthropic disturbances could also limit the species’ northern post-glaciation migration.


Introduction
In recent years, concerns about the effects of climate warming on species distribution have renewed interest in studying ecological niches and predicting changes in species' distribution ranges (Parmesan and Yohe 2003, Root et al. 2003, Chen et al. 2011. Since climate is an important factor controlling species distribution Daw-son 2003, Thomas 2010), up-to-date information about climatic niches is needed to model future species distributions under climate change. Anticipating changes in species distribution and documenting potential associated factors are key to adapting ecosystem management, particularly for keystone species.
Amphibians, particularly plethodontid sala-manders, are an important faunal component of North American forests (Burton andLikens 1975, Semlitsch et al. 2014). Many studies have revealed the importance of salamanders in food webs and their contribution to nutrient cycling (Walton et al. 2006, Best andWelsh 2014) and carbon sequestration (Best and Welsh 2014). Given their physiological requirements, they are generally very sensitive to habitat disturbance (deMaynadier and Hunter 1995). Recently, some studies have reported declining abundance of some plethodontid species in North America (Milanovich et al. 2010, Caruso and Lips 2013, Kroschel et al. 2014, pointing to climate change as a contributing factor for the observed decline. Other studies have predicted negative effects of climate change on the survival of many salamander species in the northeastern United States (Milanovich et al. 2010, Sutton et al. 2015 as well as on other amphibian species (Parmesan and Yohe 2003). Nevertheless, projected species distributions of amphibians are rarely presented (Early and Sax 2011).
The eastern red-backed salamander (Plethodon cinereus) is one of the most abundant vertebrates in the forests of northeastern North America, where it is often considered as a keystone species (Burton andLikens 1975, Walton 2013). Most of its present distribution range (Moore and Ouellet 2015) is encompassed in a large territory that, until recently on a geological timescale, was glaciated and uninhabitable. Although some recent studies have discussed the possible effect of climate on the distribution of the color phenotypes (Gibbs and Karraker 2006, Moore and Ouellet 2015, Cosentino et al. 2017) of this polymorphic species Ouellet 2014, Ouellet andMoore 2016), little is known about its climatic niche. Previously, we assembled a large dataset for this species (236,109 observations in 1148 localities and dating from 1880 to 2013) to test color morphs in P. cinereus populations as an indicator of climate change (Moore and Ouellet 2015). Given this objective, we excluded many observations that did not contain morph information.
In the present study, we update this data-base to include all observations available for P. cinereus in North America. Our main objective is to establish the current climatic niche (regarding temperature and precipitation) and distribution for P. cinereus, and to predict how the species' distribution could shift in a context of climate change. We concentrated on the northern part of its distribution range because the post-glaciation migration of the species is likely not over, and because of the potential effect of climate warming, including changes in snow precipitation and depth, on the species' northern distribution. We think that the geographical distribution of species that take refuge in the soil during the winter, such as P. cinereus, could suffer from climate change, since more frequent and severe soil frosts in the absence of sufficient snow cover on the forest floor could reduce their distribution, particularly in the northern part of their range.
Since other factors such as land use, biotic interactions, evolutionary change and dispersal ability (Davis et al. 1998, Pearson andDawson 2003) may also determine species distribution, we documented those that could affect the distribution of P. cinereus in a context of climate change. Our analysis, based on a compilation of 400,090 observations in 8302 localities and dating from 1913 to 2013 across the species' range, should help to better predict the future distribution of this important forest species.

Data compilation and current distribution of Plethodon cinereus
An exhaustive literature review was conducted to compile all available data on the presence of P. cinereus in North America. More than 400 people or institutions were contacted to complement the published literature (Moore and Ouellet 2015, Appendix S1) and to access unpublished field records, when available. The HerpNET/VertNet data portal was also accessed to obtain additional records. Most data were collected at the site scale, but some were collected at the township or the county scale. For each locality, coordinates (latitude, longitude) and sampling period were noted. When coordinates were unavailable from the literature or the authors of the inventory, they were determined using Google Earth (v.7.1.1.1580Earth (v.7.1.1. , 2013. Elevation was calculated using a digital elevation model (1 km resolution).
Climatic niche of Plethodon cinereus and its potential current northern distribution A species' climatic niche can be defined as the broad-scale temperature and precipitation conditions where it thrives (Soberón 2007). Climatic variables compiled to evaluate the climatic niche of P. cinereus were minimum, maximum, mean and median mean annual temperatures, number of frost-free days, and total precipitation. Mean climate values for the sampling period were used if salamanders were monitored for more than 1 year. If an inventory was conducted for only 1 year, we also calculated the mean climate values for 5 years before and after the observation date. This allowed us to take into account the likely presence of salamanders prior to and after the observation for a given year, and to normalize extreme values. Climate data for years prior to 1913 were not available. We estimated the average yearly climate conditions of each site using BioSIM software (v.10.2.5.39, 2013) developed by the Canadian Forest Service (Régnière and St-Amant 2007). These estimates are based on climatic variables interpolated from 8 nearby weather stations in Canada and the United States , adjusted for differences in elevation, latitude and longitude.
Snow precipitation from 1950 to 2014 was also compiled using BioSIM to evaluate annual snow accumulation pattern during this period. From 53 (for Prince Edward Island) to 1137 stations (for Québec) were used (mean: 460 stations) to calculate the mean amount of snow precipitation by province or state. Precipitation was assumed to be snow (in mm of rain equivalent) when daily mean temperature was below 0 °C. Regressions were calculated by province or state, taking into account the repeated measurements on a same station.
The current northern limit of P. cinereus' potential distribution was modeled as described later in the text. Moreover, assuming that the total annual precipitation in this area (Jobidon et al. 2015) lies within the range observed elsewhere for this species, we estimated the potential current northern distribution using the current isotherm for the grouped minimum annual mean of daily mean temperatures (−0.9 °C).
The species distribution model (SDM) was developed using the Biomod2 package (Thuiller et al. 2016) in R software (R Core Team 2016). This package models species distribution within an ensemble forecasting framework. Many modeling techniques are offered to simulate the relationships between a given species and its environment. We used various modeling algorithms to relate species occurrence to environment within When only presence data are available, an alternative is often used, such as creating pseudoabsences, in which points without locality records are used instead of absences (Chefaoui andLobo 2008, Barbet-Massin et al. 2012). Various strategies are available to select pseudo-absence points with Biomod2. Chefaoui and Lobo (2008) obtained results that show that the pseudo-absence selection method greatly influences the obtained model and generates different model predictions in the gradient between potential and realized distributions. Nonetheless, comparisons of various SDMs show that presence-absence models tend to perform better than presence-only models (Elith et al. 2006). Pseudo-absence models can be considered as an intermediate methodological approach between presence-only and presenceabsence distribution models. The random selection of pseudo-absences can be a good method because it includes many absences near the envi-ronmental domain of the presences in the modeling process (Chefaoui and Lobo 2008).
First, we considered that the species was present in grid cells with one or multiple observations. Then, we created 10 different sets of pseudo-absence points using the same number of pseudo-absence points as the number of presences. We used the random pseudo-absences in northeastern North America background data to select the pseudo-absences. As there is no evaluation data available, we created 10 random partitionings and divided them into 2 subsets: 70% to calibrate and train the models, and the remaining 30% to test them (Araújo et al. 2005). We tested and evaluated each model (8 modeling algorithms × 10 pseudo-absence repetitions × 10 partitionings, for a total of 800 models) according to the True Skill Statistic (TSS) and the Area Under the Receiver Operating Characteristic curve (AUROC).
After calibrating and evaluating the models, we projected the species' potential distribution over space and time. Consensus current projection is basically the mean of all 800 projections weighted by their TSS evaluation score. To establish the global consensus future projection, we used the mean of all 12 future climates × 4 GHG emissions scenarios × consensus future (weighted 8 modeling algorithms × 10 pseudo-absence repetitions × 10 partitionings). We transformed each probability of occurrence into presence or absence data using the threshold that maximised the model's accuracy according to the TSS score, which Jiménez-Valverde and Lobo (2007) showed to produce the most accurate predictions. We generated the response curves of every model and of the consensus current model using 100 points varying across the range of the variable of interest, while holding the other variables constant at their median value. A plot of these predictions thus reflects the effects of the selected variable. We also added four future modeled distribution of P. cinereus in the Supporting Information (Figs. S1, S2, S3, and S4), based on the four individual RCPs of GHG emissions scenarios described previously.
Moreover, we calculated the overlap be- tween current and future models using similarity "I" and "D" statistics (Schoener 1968, Warren et al. 2008), which range from 0 (niche models have no overlap) to 1 (niche models are identical). We also used distribution change calculations to evaluate areas of stability, expansion, and loss between present and future estimated distributions (see Box 3 in Guisan et al. 2014).

Climatic niche
Climatic data indicates that the species has a relatively wide climatic niche; P. cinereus experiences varying weather conditions across its range and can tolerate a wide range of seasonal temperatures (  Figure 2 represents the mean and standard deviations of evaluation scores for the various modeling algorithms used with Biomod2 (Thuiller et al. 2016). The models' predictive accuracies varied among modeling algorithms, with RF showing the best performance (as shown by its high evaluation scores for both TSS and ROC, Fig. 2). With a permutation procedure, Biomod2 shows the relative importance of the variables in the models (Thuiller et al. 2016). Maximum temperature of the warmest month (BIO5), minimum temperature of the coldest month (BIO6), and precipitation during the coldest quarter (BIO19) are the most important (Table 3). The probability of presence of P. cinereus is influenced by the maximum temperature of warmest month (BIO5) and reaches a maximum probability around 24 °C (Fig. 3). The salamander's probability of presence decrea-ses with the minimum temperature of coldest month (BIO6), reaching a minimum below −20 °C. Finally, its probability of presence is greatest when precipitation of the coldest quarter (BIO19, likely snow in many areas) reaches about 180 mm of rain equivalent. where the eastern red-backed salamander (Plethodon cinereus) was found for 400,090 observations dating from 1913 to 2013. Triangles show the mean annual temperature by province or state. Dashed lines indicate the mean minimal and maximal temperatures (all data) for a period ranging at least 5 years before and after the observation date. Province/state abbreviations and number of sites/populations for each province or state are given in Table 2.

Current, potential, and future geographic distribution of Plethodon cinereus
Data compilation reveals an extent of occurrence (commonly called EOO) for P. cinereus of approximately 1 713 000 km 2 (Fig. 4). The species has been reported at latitudes ranging from 35.3° N (North Carolina) to 50.9° N (Québec), and at longitudes ranging from 59.8° W (Nova Scotia) to 93.6° W (Minnesota). The species is found at elevations ranging from 0 to 1652 m (Mount Rogers, Virginia) (mean = 428 m, median = 299 m), mainly in temperate forests and more rarely in boreal forests.
The modeled current potential distribution of P. cinereus generally fits well with its current observed distribution (Fig. 5), but does not include some northern locations. In this context, we also used the current isotherm for the grouped minimum annual mean of daily mean temperatures (−0.9 °C) to evaluate the current potential northern distribution of P. cinereus. This exercise revealed a wide gap of up to 300 km in Québec and Ontario, between the limits of the current and of the potential northern distributions, that is unoccupied by P. cinereus. In the south, the model predicts that the species, in principle, could be found further south and southwest than its current range, and also that P. cinereus could be found in Newfoundland, Canada, which is not presently the case.
Essentially, the consensus modeled future distribution predicts that the presence of P. cinereus will decrease in the south of its range and increase in the north for the 2061-2080 period (Fig. 6). When taken individually, the 4 RCPs of GHG emissions scenarios were relatively similar to one another (Figs. S1-S4) and to the consensus modeled future distribution (Fig. 6), except RCP8.5, the most pessimistic scenario.
The gap between the modelled future distribution's northern limit and that of the current observed distribution increases to up to 400 km. Similarity statistics show that present and future modelled distributions are different (I = 0.745, D = 0.604, P < 0.001). Calculated areas of stability (0.30), expansion (0.30), and loss (0.40) between present and future modelled distributions reveal that the future species distribution range could be smaller than the current one.   Table 2.   Table 2.

Current distribution and climatic niche of Plethodon cinereus
Data compilation shows that P. cinereus has a broad geographic distribution, as expressed by its wide extent of occurrence (Fig. 4), and that the species can be found at elevations as high as 1652 m. Moreover, climatic data indicates that the species experiences diverse weather conditions, both geographically and seasonally (Table 2; Fig. 1). In a recent study, Quintero and Wiens (2013a) suggested that local seasonal temperature variations are the main factors determining the niche breadth of some plethodontid species. Not surprisingly, and as shown by annual and grouped extreme climatic values in areas where P. cinereus is found, low-temperature extremes were observed in the Canadian provinces of Ontario and Québec, and high-temperature extremes were observed in the American states of Indiana, Virginia and North Carolina, in the southern part of the species' range (Table 2; Fig. 1a). However, the effective winter temperature for this species is not well known, since the salamander generally hibernates underground. Considering the species' relatively poor dispersal ability, its wide climatic niche is probably best explained by behavioral plasticity regarding habitat use, especially foraging strategy. Salamanders burrow to protect themselves against extreme weather conditions such as summer heat and winter cold (Taub 1961). Sunday et al. (2014) mentioned that a species' vulnerability to climate warming and extreme events depends on its behavioral versatility in habitat use and on the energetic consequences of thermal retreats. Storey and Storey (1986) reported that P. cinereus is intolerant to freezing and does not survive below −1.5 °C. In winter, snow cover may isolate and protect soils and salamanders from cold winter temperatures in some areas, provided that snowfall precedes periods of extreme cold. Soil under snow cover is unlikely to get colder than 0 °C (Houle et al. 2002). Taub (1961) observed that hibernation survival in P. cinereus depends on winter severity and snow depth. Interestingly, our analysis revealed that snowfall over the last 65 years decreased significantly (by 5-31%; P ≤ 0.006; data not shown) across the species' range in all eastern Canadian provinces and in the states of Minnesota, Michigan, and Maine. Moreover, our modeling procedure identified precipitation during the coldest quarter (likely snow precipitation in many areas) as one of the most important variables in the models; salamander presence is maximal when the precipitation during the coldest quarter reaches approximately 180 mm of rain equivalent, but decreases sharply when levels drop below this value (Fig. 3). This reinforces the idea that snow precipitation and depth could be a limiting factor in the distribution of P. cinereus. In line with these observations, some studies project that climate change will decrease snow depth and snow cover duration in forests of northeastern North America. This in turn could lead to more frequent and severe soil frosts in these ecosystems (Isard andSchaetzl 1998, Campbell et al. 2005). Such a scenario could affect P. cinereus populations in the northern part of its range and reduce its distribution.

Unfilled current potential distribution
Part of the current potential distribution of P. cinereus in Ontario and Québec is unoccupied by the species (Fig. 5). This situation represents an interesting opportunity for future research on factors that might limit the dispersion of P. cinereus. Since total annual (Fig. 1c) and snow (Fig. S5) precipitation in this area are within the range observed elsewhere for this species (Jobidon et al. 2015), lack of rain precipitation or snow accumulation is unlikely to explain its absence in this area. Low urbanization and a poorly developed road network could also have limited surveys and led to an underestimation of the species' abundance. Moreover, given that P. cinereus is sensitive to forest disturbances (deMaynadier and Hunter 1995), intense logging and recurrent fires in boreal forests in the northern part of the species' range (Stocks et al. 2002, Bouchard and Pothier 2011, Boucher et al. 2014) could have created large areas of unsuitable habitats for this species.
In northern areas such as the Hudson Plains (Fig. 4), sparse vegetation due to physical con-straints (e.g., organic deposits, wetlands, bedrock without or with very thin surficial deposits, Robitaille et al. 2015) may also explain the lack of observations of P. cinereus.

Future distribution and colonization rate
We observe a gap of up to 400 km between the species' consensus modeled future distribution for 2061-2080 and the northern limits of its current observed distribution (Fig. 6). Assuming a mean colonization rate of approximately 100 m per year (see comments below), P. cinereus would need approximately 4000 years to reach the northern limit of its modeled future distribution. We estimated this colonization rate using the maximal distance reached by P. cinereus (~1300 km, Nova Scotia) since the last glaciation (Fig. 6), assuming approximately 15 000 years of range expansion (Holman 1995). This implies a colonization rate of approximately 90 m per year. Using a similar approach, Cabe et al. (2007) estimated the rate of expansion of P. cinereus at 80 m per year. These estimates seem to be conservative, since other studies have reported a dispersal ability of up to 143 m within a single season (Kleeberger and Werner 1982, Woolbright and Martin 2014, Sterrett et al. 2015. The northern limit of the future distribution remains theoretical, since many limiting factors (previously described) could hinder the movement of P. cinereus to the north. A lack of forest connectivity could accentuate the problem. At best, the species could become restricted to areas with less frequent disturbances (e.g., shorelines) and have a more fragmented distribution as a result. A recent study suggested that climate-induced environmental changes can play a dominant role in driving isolation and divergence of Plethodon species (Thesing et al. 2016).
Regardless of these limitations, rates of projected climate change largely exceed those of past climatic niche evolution for vertebrate species (Quintero and Wiens 2013b). A warming trend of up to 2.0 °C in annual temperatures has already been reported within the range of P. cinereus over the last century Karraker 2006, Houle et al. 2007). Nonetheless, there is no indication that P. cinereus populations have been at risk in recent decades (Jaeger 1980, Kroschel et al. 2014. A recent study involving plethodontid salamanders in the southern Appalachian Highlands (southeastern United States) projected that species with more southerly and smaller ranges risk greater projected habitat loss (Milanovich et al. 2010). Sutton et al. (2015) predicted that all the conservation-priority salamander species evaluated in their study could lose 3 to 100% of their current distribution in the northeastern United States by 2050. Given its relatively high climatic tolerance, P. cinereus could be less vulnerable to climate change than southern plethodontid species. Adams et al. (2007) suggest that the ecological tolerance of this species may help it adapt more rapidly to changing environmental conditions. In contrast, others suggested that projected increases in regional drought and temperatures will negatively affect P. cinereus persistence (Muñoz et al. 2016).
In the present study, the climate-modeled future distribution (consensus) suggests that the presence of P. cinereus could decrease in the south of its range and increase in the north for the 2061-2080 period (Fig. 6). These changes could be even more important with RCP8.5, the most pessimistic scenario (Fig. S4). Also, our calculation of the possible future distribution changes reveal a 10% difference between range expansion (0.30) and loss (0.40). This, combined with the potential presence of habitats unsuitable for the species' colonization in the north as well as with possible changes in interspecific interactions in the south (Smith and Pough 1994, Adams and Rohlf 2000, Price and Shields 2002, Arif et al. 2007, Jaeger et al. 2016, could induce a contraction of the species' range. Moreover, hybridization, evolution and interspecific similarities (Thurow 1956, 1957, Highton and Webster 1976, Petranka 1998, Sites et al. 2004, Lehtinen et al. 2016) could complicate future evaluation of the species' southern range.

Conclusion
Our study shows that P. cinereus has a broad climatic niche. It can live in various geographic and climatic conditions and can tolerate a wide range of seasonal temperatures. Therefore, this species could be less vulnerable to climate change than others with narrower climatic niches. The species' range could nonetheless contract in the context of climate change, given the potential presence of habitats unsuitable for the species' colonization in the north and possible changes in interspecific interactions in the south. Moreover, physical environment as well as natural and anthropic disturbances may already limit northern post-glaciation migration of this species, regardless of climate warming.