Would future climate warming cause zoonotic diseases to spread over long distances?

Source of data on the distribution of the three-toed jerboa

The distribution data of the three-toed jerboa was obtained from the Global Biodiversity Information Facility (http://www.gbif.org), the National Specimen Information Infrastructure (http://www.nsii.org.cn), a review of existing literature, and field data accumulated by our research group over the past 21 years. A total of 646 distribution points of Dipus sagitta were obtained. The buffer module of the ArcGIS software was used to eliminate the distributed data points, ensuring that there was only one distribution point in each grid to avoid the over-fitting phenomenon caused by too many distribution points. The spatial resolution of environmental data was 2.5 min (about 4.5 km), and the buffer diameter was set to 10 km. When the distance between two distribution points was less than 10 km, the buffer area overlapped, and only one distribution point was reserved (Wang, Xu & Li, 2017). After this process, a total of 216 effective distribution points were reserved (Fig. 1). The longitude and latitude coordinates of the sample were then converted into .csv format, which was used to construct the MaxEnt Model.

Note: For the distribution data of Dipus sagitta, please refer to the supporting materials.

Environmental data sources

The 19 bioclimatic variables collected in this study (Table 1) were all derived from the World Climate Database (http://www.worldclim.org/). The coordinate system was WGS84, the grid size was 25 km², and the data spatial resolution was 2.5 min. The time range of contemporary climate data was from 1960 to 2020. The Last Glacial Maximum (LGM, about 21,000 years ago), mid-Holocene (about 6,000 years ago), and future climate scenarios used the CCSM4 universal climate system model developed by the National Center for Atmospheric Research (NCAR). For future climate data, two greenhouse gas emission scenarios, RCP4.5 and RCP8.5 (RCP = Representative Concentration Pathways), from the Fifth Assessment Report of the IPCC were selected to represent a low and high impact of rising greenhouse gas concentrations on the future climate, respectively.

Filtering and processing of environment variables

The LGM, mid-Holocene, current, and future climate scenarios all included 19 bioclimatic variables (Table 1). Since some of the bioclimatic variables are highly correlated, in order to avoid overfitting of the model caused by the multicollinearity of environmental variables (Lebedev et al., 2018; Michael, 2003), DIVA-GIS 7.5 software (http://www.diva-gis.org) was used to extract information on the 19 climate variables at 216 distribution points. R software was used to carry out a Pearson correlation analysis (Pearson et al., 2006). Environmental variables with correlation coefficient r < 0.8 were retained, and environmental variables with correlation coefficient r > 0.8 were selected to have a larger contribution rate in the initial model test (Fig. 2). After final screening, seven environmental variables were obtained: mean annual temperature (Bio01), mean diurnal range (Bio02), temperature seasonality (Bio04), mean temperature of wettest quarter (Bio08), annual precipitation (Bio12), coefficient of variation of precipitation seasonality (Bio15), and precipitation of coldest quarter (Bio19).

Model building

Data from the 216 Dipus sagitta distribution points and seven climate variables in different periods were introduced into MaxEnt 3.3.3 to predict the potential distribution of Dipus sagitta under climate change conditions. A total of 75% of the distribution points of Dipus sagitta were randomly selected as training sets for model construction, and the remaining 25% were used as test sets for model verification. Cross-validation was adopted, meaning the species data were randomly divided into 10 parts, with one part randomly selected as the test set each time, and the remaining nine parts used as training sets, repeated ten times. By default, the maximum number of iterations was 500 and the maximum number of background points was 10,000.

A key parameter setting in MaxEnt is feature class, which has five feature types: linear feature (L), quadratic feature (Q), fragmented feature (H), product feature (P), and threshold feature (T). With more distribution points, MaxEnt uses more of these feature types by default (Elith, Kearney & Phillips, 2010; Wang et al., 2017; Kong, Li & Zou, 2019), with little impact on overall performance. Because the distribution points of the selected species in this study were >80, the default feature selection of MaxEnt was selected.

This data was used to plot environmental variable response curves to evaluate the impact of each climate variable on the model prediction results. Receiver operating characteristic (ROC) curves were created and the jack-knife method was used to assess the importance, or contribution, of each environmental variable to distribution gain. The jack-knife method calculates the training scores when simulating “with only one variable,” “without variables,” and “with all variables,” respectively. When the “with only one variable” had a higher score, it indicated that the environmental factor being used as the variable had a higher predictive ability and thus a greater contribution to species distribution. When the training scoring ability of the “with only one variable” model decreased, it indicated that the variable had more unique information and was more important to species distribution (Phillips & Miroslav, 2008). Using this method, the dominant factors affecting the species distribution were determined, and the remaining parameters were kept at their default settings. The model then output the results in the form of logistic output, from which the ASC II (American Standard Code for Information Interchange) file of the average of 10 operations were selected. Raster values giving the distribution probabilities of the species were used as logical values reflecting the extent to which each grid in the target area met the actual ecological niche of the species. These results were then input into ArcGIS 10.2 software (Yang et al., 2013). The automatic classification method in the reclassification tool was used to divide the suitability of the distribution areas into four grades: non-suitable area, somewhat suitable area, moderately suitable area, and highly suitable area.

Geospatial analysis

ArcGIS 10.2 was used to count the different grades of suitable areas in different periods, and SDM toolbox 2.4 tool was used to calculate the potential distribution areas of Dipus sagitta in different periods (Brown, 2014). Using the “Reclass” function in ArcGIS 10.2, the grid values of suitable (including somewhat, moderate, and highly suitable areas) and non-suitable areas of Dipus sagitta were modified to 1 and 0 respectively, then SDM toolbox was added, and the subdirectory of “MaxEnt Tools” in the “SDM Tools” module was selected. The “Distribution Changes Between Binary SDMs” tool was used to calculate changes in potential distribution areas in each period (LGM maximum to mid-Holocene, mid-Holocene to current, current to RCP 4.5-2070, and current to RCP 8.5-2070) and obtain the distribution expansion area, stable area, and contraction area.

Model accuracy evaluation

Receiver operating characteristics (ROC) take every value of the prediction result as a possible judgment threshold and calculates the corresponding sensitivity and specificity, and then evaluates the accuracy of the model. Because the area under the ROC curve (AUC) is not affected by the judgment threshold, it is recognized as the best index to evaluate the prediction accuracy of the model. The value range of AUC is [0, 1]. The larger the value, the farther it is away from the random distribution, indicating a more accurate prediction effect. The evaluation criteria are: 0.7–0.8 is more accurate, 0.8–0.9 is accurate, and 0.9–1.0 is very accurate (Phillips, Anderson & Schapire, 2006).

Because the subject working curve generated by MaxEnt was fuzzy, to improve the clarity, the origin picture numerical drawing toolbox was used to extract the result of the MaxEnt subject working curve. This extracted result was then imported into MATLAB and the plot function was used to redraw the ROC curve. This study also chose a partial receiver operating character (PROC) measurement method to evaluate the prediction performance of the model (http://shiny.conabio.gob.mx:3838/nichetoolb2/). Finally, the model was constructed based on the distribution data of all species. The PROC method uses AUC ratio to evaluate the model. An AUC ratio >1 indicates that the model has better relative random prediction results, while AUC ratio ≤ 1 indicates that the model has worse relative random prediction results (Fan et al., 2019). PAUC was calculated by NicheToolbox (http://shiny.conabio.gob.mx:3838/nichetoolb2/) with 1,000 iterations, E = 0.05.

Laboratory animal ethics

The animal study was reviewed and approved by the Research Ethics Review Committee of Inner Mongolian Agricultural University (NND2023081).

MaxEnt model accuracy verification

The MaxEnt model was used to simulate the potential habitat of Dipus sagitta, and the ROC was used as a measure of the prediction accuracy of the model. The evaluation results of ROC show that the average AUC (referring to the area surrounded by ROC and abscissa axis) of the prediction model of suitable living area of Dipus sagitta was greater than 0.9 (Fig. 3), the average AUC value of the evaluation results of the model PROC performance test was 0.906 (P < 0.001) (Fig. 4), and the AUC ratio was 1.813, which is much greater than 1, indicating that the evaluation results of the prediction model and random model were significantly different, and the performance of the prediction model was significantly stronger than the random model. This suggests that the mean value of the MaxEnt model can be used to reflect the distribution of species with high accuracy.

The importance of environmental variables affecting the distribution of the three-toed jerboa

The influence of seven environmental variables on the potential future distribution area of the three-toed jerboa was assessed using the jack-knife method (Table 2). The top three variables by percentage contribution were: precipitation of the coldest quarter (Bio19, 38.77%), for which the most suitable range was 0–22.79 mm; temperature seasonality (standard deviation; Bio04, 26.15%), for which the most suitable range was 86.72–250; and mean annual temperature (Bio01, 23.76%), for which the most suitable range was 1.37–19.73 °C. The cumulative contribution rate of these three factors was 88.68%. The top three permutation importance values were: the standard deviation of temperature (Bio04, 32.01%), mean annual temperature (Bio01, 30.21%), and precipitation of the coldest quarter (Bio19, 19.69%). The cumulative importance value of these three factors was 81.92%.

The simulation results of the MaxEnt model showed that, among the seven environmental variables in this study, both precipitation and temperature had a certain degree of influence on the distribution area of the three-toed jerboa. The ranking of contribution rates showed that the precipitation factor was more important, but replacement importance value results showed that the temperature factor was more important. These results suggest that the two main factors affecting the contemporary geographic distribution of the three-toed jerboa are precipitation and temperature. Precipitation of the coldest quarter was the environmental factor with the highest contribution rate. The simulation results showed that the most favorable range of precipitation of the coldest quarter for the survival of the three-towed jerboa was 0–22.79 mm, which is in line with this animal being highly adapted to arid conditions (IUCN, 2016).

Potential suitable habitats of the three-toed jerboa in different climatic backgrounds

The potential distribution area of the three-toed jerboa under the current climate background was simulated by the MaxEnt model (Fig. 5). This area covered almost all the sample points, indicating that the simulated distribution area of the three-toed jerboa was in close agreement with its actual distribution. The distribution of highly suitable areas was relatively concentrated (brown area), with a total area of 3.25 ×10⁶ km². This region mainly included China’s Xinjiang region, central and western Inner Mongolia, Shaanxi, Shanxi, and Ningxia. In Mongolia, the highly suitable areas were mainly distributed in the western region. There were also sporadic distributions in Kyrgyzstan. The moderately suitable areas (dark green area) were mainly concentrated in southern Kazakhstan, western Uzbekistan, eastern Turkmenistan, and along the edge of the highly suitable areas, with a total moderately suitable area of 3.51 ×10⁶ km². The somewhat suitable areas (yellow area) were scattered throughout northern Kazakhstan, southern Turkmenistan, southern Uzbekistan, and the borders between Russia, China, and Mongolia. The simulated distribution area also showed areas suitable for the survival of the three-toed jerboa in the Midwest region of North America, but there is no actual distribution of the three-toed jerboa in that region. The main species of rat in the Midwest region of North America is the kangaroo rat (Dipodomys spectabilis), which is very similar to the three-toed jerboa in both morphological and physiological characteristics (Li et al., 2007).

The suitable distribution area of the three-toed jerboa varied greatly across the Last Glacial Maximum, the mid-Holocene, current climate conditions, and future climate scenarios up to the 2070s (Table 3). The highly suitable area for the three-toed jerboa in the Last Glacial Maximum was the smallest, at 1.72 ×10⁶ km², and then it continued to spread and increase, reaching 3.25 ×10⁶ km² in the current climate period, an increase of 88.95%. By the 2070s, under the RCP8.5 scenario, highly suitable areas could be reduced to 2.6 ×10⁶ km². During the Last Glacial Maximum, the moderately suitable area of the three-toed jerboa was small, 2.03 ×10⁶ km², after which it began to increase. The moderately suitable area in the mid-Holocene expanded to 4.17 ×10⁶ km², an increase of 105.49%. From the mid-Holocene to the 2070s under RCP8.5, the area of suitable habitat for the three-toed jerboa shrank slightly, reducing by 0.42 ×10⁶ km² to 3.75 ×10⁶ km², which was still higher than in the Last Glacial Maximum. The somewhat suitable area for the three-toed jerboa in the mid-Holocene was the largest at 9.3 ×10⁶ km², an increase of 99.14%, or 4.63 ×10⁶ km² compared with the Last Glacial Maximum. After this, the somewhat suitable area for three-toed jerboa shrank gradually, and could reduce to 7.26 ×10⁶ km² in the 2070s under RCP8.5.

Spatial changes in the total suitable habitat of the three-toed jerboa under different climate scenarios

To better analyze changes in the suitable area for the three-toed jerboa, the three gradients of suitable area (somewhat, moderately, and highly) were considered together as the total suitable area. The current three-toed jerboa distribution layer was superimposed on the Last Glacial Maximum and the mid-Holocene prediction layer, and the differences in distribution were analyzed (Table 4, Figs. 6, 7). The results showed that from the Last Glacial Maximum to the mid-Holocene, the spatial pattern change rate of the three-toed jerboa was 93.91%, and the expansion rate reached 111.67%. The total suitable habitat has increased on a large scale. Kazakhstan, China, Turkmenistan, Afghanistan, Iran, and Mongolia all experienced varying degrees of expansion. The change rate of the spatial distribution pattern of the three-toed jerboa from the mid-Holocene to the present was 10.19%, with a shrinkage rate of 19.11%. The shrinking areas were mainly concentrated in Afghanistan and Iran, and there were also sporadic shrinkages in central and western China. The expansion rate was 5.49%, concentrated primarily at the junction of Afghanistan, Iran, and Turkmenistan, with some expansion in the southern area of Afghanistan, indicating that the three-toed jerboa began to migrate southward.

The current three-toed jerboa distribution layer was superimposed on the forecast layer under the two carbon emission scenarios for the 2070s, RCP4.5 and RCP8.5, and the change trend was analyzed (Table 4, Figs. 8, 9). Under these future climate change scenarios, the suitable habitat of the three-toed jerboa generally shrinks in the east, west, and south and expands to the north. The areas of expansion are mainly located at the junction of Mongolia and Russia, the junction of eastern Inner Mongolia and Russia, and the northern border of Kazakhstan, while the shrinking distribution areas are mainly concentrated in the southern border of Inner Mongolia, China, eastern Jilin Province, western Kazakhstan, Uzbekistan, Kyrgyzstan, and southern Turkmenistan. Under the RCP4.5 emission scenario, the shrinkage rate is 1.61%, and under the RCP8.5 emission scenario, the shrinkage rate is 1.86%. These results indicate that in order to adapt to high-concentration CO₂ emissions, in future climate warming scenarios, the distribution area of the three-toed jerboa will migrate and expand to higher latitudes.