Geographic characteristics of sable (Martes zibellina) distribution over time in Northeast China

Abstract Understanding historical context can help clarify the ecological and biogeographic characteristics of species population changes. The sable (Martes zibellina) population has decreased dramatically in Northeast China since the l950s, and understanding the changes in its distribution over time is necessary to support conservation efforts. To achieve this goal, we integrated ecological niche modeling and historical records of sables to estimate the magnitude of change in their distribution over time. Our results revealed a 51.71% reduction in their distribution in 2000–2016 compared with the potential distribution in the 1950s. This reduction was related to climate change (Pearson's correlation: Bio1, −.962, p < .01; Bio2, −.962, p < .01; Bio5, .817, p < .05; Bio6, .847, p < .05) and human population size (−.956, p < .01). The sable population tended to migrate in different directions and elevations over time in different areas due to climate change: In the Greater Khingan Mountains, they moved northward and to lower elevations; in the Lesser Khingan Mountains, they moved northward; and in the Changbai Mountains, they move southward and to higher elevations. Active conservation strategies should be considered in locations where sable populations have migrated or may migrate to.


| INTRODUCTION
The sable (Martes zibellina), a species of marten (family: Mustelidae), is a medium-sized carnivore distributed widely in Russia, Mongolia, China, North Korea, and Japan (Proulx et al., 2005). It inhabits forests and is associated with coniferous taiga forests. In China, it is distributed in Northern Xinjiang and Northeast China (Monakhov, 2001b(Monakhov, , 2005Proulx et al., 2005). Due to commercial hunting and deforestation, the sable population has decreased dramatically in Northeast China since the l950s (1950-1959Kashtanov et al., 2015;Monakhov, 2005Monakhov, , 2011a and there is an increasing need for its conservation.
Understanding changes in its distribution over time is necessary to support basic conservation efforts.
We integrated ecological niche modeling with historical sable records to estimate the magnitude of change in their distribution over time. To achieve this, we (1) created a historical dataset of sable distribution from multiple resources (e.g., gazetteers and journal articles); (2) reconstructed climate data to represent the environmental variables; (3) estimated the potential habitat for each period using MaxEnt; and (4) evaluated changes in their distribution range and possible impact factors. Our results offer an opportunity to understand the ecological and biogeographic characteristics of species population decline, thereby improving the predictive power of conservation management for sables and other species in similar ecology niches.

| Study area
Northeast China includes Heilongjiang Province, Jilin Province, and the northeast corner of Inner Mongolia (N40°5′-53°17′, E115°30′-135°06′; Yang et al., 2016;Zhang et al., 2016). The climate is a continental monsoon climate with a negative water balance. The annual precipitation is 400-1,000 mm, and the annual average temperature is 1-4°C, with a north-south temperature gradient of 25°C. The region includes coniferous forest, broadleaf mixed forest, secondary forest, woodland shrub and marshy grass areas, and contains more than 2,500 plant species, with a tree cover of ~42.9% (Xiaofeng et al., 2011;Yang & Xu, 2003). Coniferous forest is mainly located in the Greater Khingan Mountains, where it is characterized by Larix gmelinii.
Broadleaf mixed forests are located in the Lesser Khingan Mountains, Changbai Mountains, and Wanda Mountains and are dominated by L. gmelinii, Pinus koraiensis, and Betula platyphylla (Yang & Xu, 2003).
The forested area supports numerous wildlife species, including M. zibellina and its abundant prey (e.g., moles and hares).
In this study, we obtained the boundary of Northeast China from the Thematic Database for Human-Earth System (http://www.data. ac.cn/index.asp). The area covers 938,000 km 2 with an elevation range of 2-2,576 m. The ecological data were obtained in a raster structure with a cell size of 1 km 2 (Figure 1).

| Data
Based on previous research Zhang et al., 2016), we collected historical sable data from five sources: gazetteers, fauna records, nature reserve scientific surveys, scientific research, and news (see the Table S1 in Appendix S1). Conflicting records with unsubstantiated datasets, such as those lacking relevant or detailed descriptions, were excluded from the analysis. We recorded the locality information from all of the extracted records. In total, 1,028 records were collected. To eliminate potential biases caused by clustered occurrences, we removed duplicate records within the same cell. Ultimately, records were selected and divided into different time periods: 902 records in 1950s (1950-1959), 678 in 1960s (1960-1969), 489 in 1970s (1970-1979), 363 in 1980s (1980-1989), 147 in 1990s (1990-1999), 84 in 2000s (2000-2016).
Sables may be affected by human influences (Monakhov, 2001a(Monakhov, , 2001b(Monakhov, , 2011a(Monakhov, , 2011bZhou et al., 2011), but the historical data of human activities (e.g., land cover, vegetation cover from the 1950s) are difficult to be compiled. The change of human population size can reflect the scale of urbanization, the deforestation, and the reclamation of farmland. Therefore, we collected human population size data from the first to the sixth national population censuses and several gazetteers (see the Fig. S2 in Appendix S3).

| Environmental variable selection
Daily climate data were obtained from the surface daily climate dataset (code: SURF_CLI_CHN_MUL_DAY, ver. 3.0, http://data.cma.cn), including mean temperature, maximum temperature, minimum temperature, and precipitation (recorded daily at 20:00). Datasets from 1951 to 2010, from 72 meteorological stations, were downloaded. Environmental variables for 12 months per 10 years were calculated with ANUSPLIN ver. 4.36. This package facilitates the transparent analysis and interpolation of complex multivariate data using thin plate smoothing splines, which were suitable for this research. For each period, bioclimatic variables were derived from the monthly temperature and rainfall values to generate more biologically meaningful variables with the dismo package in R (R Development Core Team, Vienna, Austria). From this, we obtained 19 bioclimatic variables. Topography data were obtained from the SRTM 90 m Digital Elevation Database (http://srtm.csi.cgiar.org/).
The slope and aspect data were calculated from the topography data. In total, 22 variables were obtained for further research.
Overfitting can result from highly correlated variables; therefore, we built a Pearson correlation matrix for 22 variables. We drew a network diagram (see the Table S3 in Appendix S2) using a raster structure with the igraph package in R to facilitate variable selection (Csardi & Nepusz, 2006;Hijmans et al., 2016). In the network diagram, highly correlated variables (|r| ≥ .8) were linked by lines and clustered together. We assumed that the climate variables could constrain the distribution during months when the species face a severe ecological environment in this region. Under the assumption that extreme climate index variables (e.g., BIO5, BIO6, BIO14, etc.) represented limiting factors for the species, we prioritized those variables among the correlated cluster. Each period had a unique set of variables for the ecological niche model (Figure 1; see the Table S3 in Appendix S2).

| Species distribution modeling
We estimated the potential sable distribution in Northwest China over time using MaxEnt ver. 3.3.3k (Elith et al., 2011;Phillips et al., 2006).
We randomly selected 80% of the occurrence records as calibrations to generate the model and used the other 20% as test data to evaluate the model. The model was developed with MaxEnt ver. 3.3.3k with the following parameters: number of random selected background points = 10,000; replicated run type = subsample; replicates = 20; auto feature; convergence threshold = 0.00001; and output format = logistic. The final model was the average of these replicates. To estimate its accuracy, we used the area under the curve (AUC) of the receiver operating characteristic to assess model fit (Feng et al., 2015;Phillips et al., 2006).
In order to transform our models from environmental suitability into presence-absence distributions, we used the threshold (P) calculated with the sensitivity-specificity sum maximization approach F I G U R E 1 Approach for modeling species potential distributions (Liu et al., 2005). We transformed logistic model output to a presence-absence map for each time periods using the mean threshold of maximum training sensitivity plus specificity logistic threshold as the cutoff. Then, all the outputs were divided into two groups: Outputs above the threshold (>P) were grouped as "present," while all other values (<P) were considered to be "absent."

| Spatial analysis
It is important to quantify distribution change spatial patterns to gain a better understanding of threatening factors (Luo et al., 2015;Pacifici et al., 2015;Urban, Zarnetske, & Skelly, 2013). The predicted ranges for each period were determined as the total number of respective re- Human influence and climate are likely to be important factors affecting sable distribution (Kashtanov et al., 2015;Miyoshi & Higashi, 2004;Monakhov, 2001bMonakhov, , 2005Monakhov, , 2011a

| Threatening factors
Global warming and human influences can force animals to alter their habitats, leading to habitat loss and fragmentation (Luo et al., 2015;Rands et al., 2010;Willis, Gillson, & Knapp, 2007;). Compared with the temperature characteristics in the 1950s, the annual mean temperature, maximum temperature of the warmest month, and minimum temperature of the coldest month Occurrence point data were collected from five sources and verified (see Section 2.2). Records were obtained from occurrence points with the tool to spatially rarefy occurrence data. All results were derived from the distribution models. Model accuracy was measured from the AUC curve (mean ± standard deviation [SD]). The ranges were determined as the total of the number of related pixels.  (Table 3, specific information showed in Appendix S3).

| DISCUSSION
Our investigation using historical records to quantify the temporal and spatial dynamics of range change supports the use of multiple datasets (e.g., gazetteers, fauna records, nature reserve scientific surveys, scientific research, and news) to contribute novel insights into our understanding of the dynamics of species responses to human pressure and climate change. In addition, they can be used to track range changes across longer timescales than those usually addressed in ecology or conservation biology. Our analyses combined multiple historical datasets with historical climate and species distribution models, which has not been attempted in previous research.

| Model limitations
This study has several limitations. First, this research cannot provide completely accurate insight into the changes in sable distribution range over time. Due to our use of historical data, the results may contain sample biases resulting from data accessibility, spatial and temporal variables, and nonstandardized sampling (Turvey et al., 2015;Yang et al., 2016;Zhang et al., 2016). The predictors came from the surface daily climate dataset of 72 meteorological stations also contained sample biases resulting from spatial cluster (most of the meteorological stations next to the city) and data accuracy (the data precision increased after the 1990s

| Biogeographic characteristics over time
Based on our results, range change was particularly associated with two factors: climate and human influence. The range change was divergent in the three regions examined due to their complex terrains.
Due to the lower human population density in the Greater Khingan Mountains since the 1950s, the sable population movement toward north mainly resulted from the poleward shifts under global warming.
In the Lesser Khingan Mountains, the influence of the urbanization and climate change through time forced the sable population to shift the habitat to the north which tend to be colder and have the lower human population density. According to the landcover type in 2013 (Fig. S3 in Appendix S3), urban agglomeration and farmland can be the important threat to the sable population in the Lesser Khingan Mountains and Changbai Mountains. Therefore, our results partly demonstrate that these three regional sable populations were influenced by climate change, as well as by other factors (e.g., human population size), which is consistent with the results of Urban and Chen (Chen, Hill, Ohlemuller, Roy, & Thomas, 2011;Urban et al., 2013). This trend of range change would be continued in the near future. Moreover, we suggest that the trend of the human population size through time should be considered in species distribution modeling research.

| Conservation implications
In responding to climate change, there is a particular need for climate  Schloss et al., 2012;Willis et al., 2007). Many studies on this issue have suggested that climate influences local species abundance, community structure and biodiversity, phenology, and species range, but only a few have addressed their research using a historical perspective (Ban et al., 2014;Beaugrand, Edwards, Raybaud, Goberville, & Kirby, 2015;Davies et al., 2014;Mathur & Padalia, 2006;McClenachan et al., 2012;Rick & Lockwood, 2012;Saikia, Kalita, & Saikia, 2009). We suggest that conservation strategies based on historical information may provide a more realistic perspective. The sable population in Northeast China has divided into three regional popula- The population in the Changbai Mountains may be facing the greatest challenge, as their migration to higher elevations in response to climate change limits space for their distribution. In addition, this region has experienced strong human influences for many years (Figure 2; Fig. S2 in Appendix S3). As the vulnerability of this population, due to climate change, may be higher than the other populations, active conservation strategies should be considered, including increasing connectivity, assisted migration, ex situ conservation, and reintroduction.