Impact of Climate Change on Potential Distribution of Chinese White Pine Beetle Dendroctonus armandi in China

Dendroctonus armandi (Coleoptera: Curculionidae: Scolytidae) is a bark beetle native to China and is the most destructive forest pest in the Pinus armandii woodlands of central China. Due to ongoing climate warming, D. armandi outbreaks have become more frequent and severe. Here, we used Maxent to model its current and future potential distribution in China. Minimum temperature of the coldest month and precipitation seasonality are the two major factors constraining the current distribution of D. armandi. Currently, the suitable area of D. armandi falls within the Qinling Mountains and Daba Mountains. The total suitable area is 15.83 × 104 km2. Under future climate scenarios, the total suitable area is projected to increase slightly, while remaining within the Qinling Mountains and Daba Mountains. Among the climate scenarios, the distribution expanded the most under the maximum greenhouse gas emission scenario (representative concentration pathway (RCP) 8.5). Under all assumptions, the highly suitable area is expected to increase over time; the increase will occur in southern Shaanxi, northwest Hubei, and northeast Sichuan Provinces. By the 2050s, the highly suitable area is projected to increase by 0.82 × 104 km2. By the 2050s, the suitable climatic niche for D. armandi will increase along the Qinling Mountains and Daba Mountains, posing a major challenge for forest managers. Our findings provide information that can be used to monitor D. armandi populations, host health, and the impact of climate change, shedding light on the effectiveness of management responses.


Introduction
The Chinese white pine beetle, Dendroctonus armandi Tsai and Li (Coleoptera: Curculionidae: Scolytidae), a bark beetle native to China, is distributed in the Pinus armandii Franch. woodlands of central China [1]. As a pioneer species, this bark beetle invades healthy P. armandii individuals older than 30 years. Subsequently, secondary bark beetles attack the infected or withered hosts, resulting in landscape-level tree mortality [2]. Dendroctonus armandi inhabits and feeds on the phloem of the host [3]. When mature adults search for new host trees, they briefly leave the gallery. It is always the females that break through the first line of defense, by boring through the host's bark. The females then release aggregation pheromones to attract males for mating and colonization. The females subsequently lay eggs in the phloem, and the larvae excavate feeding tunnels in this tissue [4]. In the Qinling Mountains, the number of beetle generations varies with elevation, with ca. 2 per year below 1700 m, 1.5 per year at 1700-2150 m, and 1 per year above 2150 m [5]. Outbreaks of this beetle have occurred every year in recent decades, altering the ecosystem and causing economic losses.
Because insects are poikilotherms and have short generation times, their life history is closely related to temperature [6]. They are therefore particularly sensitive to the warming Dendroctonus armandi presence data was collected in the field, from May to October 2016-2019. The fieldwork was conducted with the assistance of the Department of Forestry Protection. The distance between the sampling sites was >20 km. When the distance between two or more occurrence points was <20 km, we recorded only one occurrence. First, we judged beetle presence by observing resin-flow from the bark, because P. armandii responds to attack by D. armandi by releasing large amounts of distinctly fluid resin. Further, we used pheromone-baited traps to detect beetle presence in areas with population outbreaks.
The traps consisted of six identical black plastic funnels aligned vertically over a cup (model YBQ-LD-004, Geruibiyuan Technology Co., Ltd., Beijing, China), similar to those commonly used to collect bark beetles. As an attractant we used <bait: mixture of terpenes and pheromones> (Sino-Czech Trading Company, Beijing, China) packed in a 15 mL slow-release polyethylene vial with a release rate of 200 mg/day. Each trap was vertically suspended from a manila rope tied between two host trees at the forest edge and positioned more than 1 m from the closest host tree. The lower end of the trap was 1.5 m above ground. The traps were spaced at least 30 m apart. Liquid paraffin was used as a killing agent in the cups. Trapped D. armandi were collected weekly, counted, and then stored in 70% ethanol. Location name, longitude, latitude, and altitude were documented for all occurrence points. In addition, occurrence records were collected from the literature [1][2][3][4][5][8][9][10]; points without accurate location data were excluded. In total, we obtained 132 point-based occurrence records from the literature and fieldwork (Figure 1).  5,[8][9][10]; points without accurate location data were excluded. In total, we obtained 132 point-based occurrence records from the literature and fieldwork (Figure 1). For the host distribution variable, we gathered presence data for the host P. armandii using three methods. First, we registered P. armandii in the field from 2016 to 2018. Raw data were processed as for the D. armandi data. Second, we obtained P. armandii occurrence data from the Global Biodiversity Information Facility (GBIF, http://www.gbif.org (accessed on 27 April 2021)) and the Chinese Virtual Herbarium (CVH, http://www.cvh.org.cn/ (accessed on 27 April 2021)). Although GBIF contains >1000 records, many lack coordinates and altitudes, or are duplicates, and were removed. Third, we consulted the literature for P. armandii occurrence records using a Chinese academic website (https://www.cnki.net/ (accessed on 27 April 2021)). From the all of the above records, we discarded 78 records with incomplete geographic information. In total, 236 P. armandii occurrence records were used to generate its potential distribution via Maxent.

Environmental Variables
Nineteen bioclimatic variables, three terrain variables, and one host variable were used to predict the potential distribution of D. armandi (Table 1). Nineteen bioclimatic variables, with a spatial resolution of 2.5 arc minutes (approximately 4.7 × 4.7 km), were freely downloaded from www.worldclim.org. Terrain variables were obtained from the United States Geological Survey (USGS, http://edcdaac.usgs.gov/gtopo30/gtopo30.html (accessed on 27 April 2021)). We used Arc-GIS v.10.2 (Environmental Systems Research Institute Inc., Redlands, CA, USA) to convert all variable layers into ASCII format for use with Maxent. We accessed the National Fundamental Geographic Information System (http: //www.ngcc.cn (accessed on 27 April 2021)) and downloaded the National Map as the analytical base map. For the host distribution variable, we gathered presence data for the host P. armandii using three methods. First, we registered P. armandii in the field from 2016 to 2018. Raw data were processed as for the D. armandi data. Second, we obtained P. armandii occurrence data from the Global Biodiversity Information Facility (GBIF, http://www.gbif.org (accessed on 27 April 2021)) and the Chinese Virtual Herbarium (CVH, http://www.cvh.org.cn/ (accessed on 27 April 2021)). Although GBIF contains >1000 records, many lack coordinates and altitudes, or are duplicates, and were removed. Third, we consulted the literature for P. armandii occurrence records using a Chinese academic website (https://www.cnki.net/ (accessed on 27 April 2021)). From the all of the above records, we discarded 78 records with incomplete geographic information. In total, 236 P. armandii occurrence records were used to generate its potential distribution via Maxent.

Environmental Variables
Nineteen bioclimatic variables, three terrain variables, and one host variable were used to predict the potential distribution of D. armandi (Table 1). Nineteen bioclimatic variables, with a spatial resolution of 2.5 arc minutes (approximately 4.7 × 4.7 km), were freely downloaded from www.worldclim.org (accessed on 27 April 2021). Terrain variables were obtained from the United States Geological Survey (USGS, http://edcdaac.usgs.gov/ gtopo30/gtopo30.html (accessed on 27 April 2021)). We used Arc-GIS v.10.2 (Environmental Systems Research Institute Inc., Redlands, CA, USA) to convert all variable layers into ASCII format for use with Maxent. We accessed the National Fundamental Geographic Information System (http://www.ngcc.cn (accessed on 27 April 2021)) and downloaded the National Map as the analytical base map. We used a dataset on future climate from the Climate Change, Agriculture and Food Security program (CCAFS, http://ccafs-climate.org/ (accessed on 27 April 2021)) to determine the future distribution of D. armandi under different climate scenarios. Representative concentration pathways (RCPs) were announced in the fifth Intergovernmental Panel on Climate Change report. These include four greenhouse-gas concentration trajectories, representing scenarios in which the total radiative forcing in 2100 is 2.6 W/m 2 , 4.5 W/m 2 , 6 W/m 2 , and 8.5 W/m 2 (in excess of 1750 W/m 2 ) [25]. We used the BCC-CSM 1.1 climate system model, developed by the Beijing Climate Center, China. After more than ten years of development, its performance and functions have been effectively improved [26]. It has been widely used for climate prediction. Of the four available RCPs, we used RCP 2.6, RCP 4.5, and RCP 8.5 to model the future distributions of D. armandi in the 2030s and 2050s.
To avoid overfitting of the models due to multicollinearity within environment variables, we selected variables with higher maximum entropy gain and nonzero regression coefficients of LASSO (least absolute shrinkage and selection) [27,28]. LASSO is a regression analysis method that performs both variable selection and regularization, in order to enhance the prediction accuracy and interpretability of the statistical model. The method regularizes model parameters by shrinking the regression coefficients, reducing some to zero. The feature selection phase occurs after the shrinkage; in this phase, every non-zero value is selected for use in the model. This improves model accuracy because coefficient shrinkage reduces variance and minimizes bias [29]. LASSO was applied using the R package glmnet (http://www.web.stanford.edu/~hastie/Papers/Glmnet_Vignette.pdf (accessed on 27 April 2021)).

Species Distribution Modelling
We used Maxent v.3.3.1 (http://www.cs.princeton.edu/wschapire/MaxEnt (accessed on 27 April 2021)) to predict the potential distribution of D. armandi. Ten-fold crossvalidation was used to train and validate the model. The occurrence dataset was randomly divided into 10 equal-sized subsets. Of these, a single subset was retained as the validation data for testing; the remaining subsets were used as training data. For each subset, 90% of the occurrence data were used to train the single model, and the remaining data were used to test the predictive ability of the model. The mathematical modelling process was provided by the Maxent function in the R package Maxent (https://www.rdocumentation. org/packages/maxnet/versions/0.1.2 (accessed on 27 April 2021)). The cross-validation process was repeated 10 times, with each of the 10 subsets used once as validation data. The 10 results were then averaged to produce a single estimation. The area under the receiver operator characteristic (ROC) curve (AUC) was used to evaluate the prediction performance of each Maxent model. When the ROC curve is at 45 • in the ROC space, the AUC value is close to 0.5, indicating that the model is a random model with an accuracy of 50%; when the AUC > 0.5, the model is more accurate than a random model; when the AUC < 0.5, the model is less accurate than a random model. The closer the AUC is to 1, the better the model performance [30,31] The impact of sampling bias on species distribution modeling must be noted [32]. Sampling intensity often varies between sites and sampling sites are usually biased toward areas which are more accessible. Oversampling in some geographic spaces can cause repetition in the ecological space when building the model; this will affect simulation of the ecological needs of the species. This kind of sampling bias will generally cause overfitting of niche models to species requirements, thus reducing model transferability [33]. In brief, these situations may result in the observed species distribution not reflecting the real distribution. To counter the sampling bias, "pseudo-absence" with the same spatial bias as the presence points is recommended to be introduced into the model. For this purpose, Arc-GIS v. 10.2 was used to generate 1000 background points throughout China with the same bias [34]. As a result, there were a total of 1132 occurrence points with the same bias in this study. Next, we selected points for model calibration using a subsampling regime to reduce sampling bias and spatial autocorrelation [35]. We generated models using all occurrence points and measured spatial autocorrelation among model pseudo-residuals (1-probability of occurrence generated by model) by calculating Moran's I at multiple distance classes [36,37]. Significance was determined using permutation tests. A minimum distance of 275 km was detected, so we selected the occurrence point that close to the centroid of each grid cell. This procedure reduced the number of occurrences to 138 points used for model calibration. The procedure greatly reduced sampling bias and spatial autocorrelation, resulting in evenly distributed occurrence points across space [36].
The output layers were imported into Arc-GIS v.10.2 for further analysis, and the final distribution map was generated. Using occurrence records from different sources may cause some sampling bias. A ten-percentile training presence logistic threshold was adopted to define the minimum probability of suitable habitat for D. armandi [38]. Based on this, the habitat suitability of D. armandi was divided into four grades: unsuitable (0-0.1), poorly suitable (0.11-0.3), moderately suitable (0.31-0.6), highly suitable (0.61-1). We then calculated suitable areas under future climate scenarios by multiplying the number of presence grid cells by their spatial resolution.

Model Performance and Variables Selection
When running the training and test dataset independently, the model had robust evaluation metrics: the five-number summary of AUC based on 10-times cross-validation (0.838, 0.878, 0.894, 0.921, 0.946) shows that the model performance was excellent. Based on the regression coefficients and percentage contributions of the environmental variables used in the model, we selected the following five variables as the most important for determining the habitat suitability of D. armandi: minimum temperature of coldest month (Bio6, 1.4832; 29.9%), precipitation seasonality (Bio15, −0.8965; 19.3%), mean temperature of driest quarter (Bio9, 0.7865; 16.8%), annual mean temperature (Bio1, −0.4614; 11.4%), precipitation of driest quarter (Bio17, −0.2488; 9.8%) ( Table 2). The response curves show the quantitative relationship between environmental variables and the habitat suitability, and the responses of five variables to D. armandi suitability are shown in Figure 2. According to the response curves, annual mean temperature (Bio1) ranges from 5.6 to 17.1 • C, minimum temperature of coldest month (Bio6) ranges from −11 to 2 • C, mean temperature of driest quarter (Bio9) ranges from −5 to 7 • C, precipitation seasonality (Bio15) ranges from 63 to 94, and precipitation of driest quarter (Bio17) ranges from 12 to 74 mm. These variables were considered to have a significant effect on the survival of D. armandi. used in the model, we selected the following five variables as the most important for determining the habitat suitability of D. armandi: minimum temperature of coldest month (Bio6, 1.4832; 29.9%), precipitation seasonality (Bio15, −0.8965; 19.3%), mean temperature of driest quarter (Bio9, 0.7865; 16.8%), annual mean temperature (Bio1, −0.4614; 11.4%), precipitation of driest quarter (Bio17, −0.2488; 9.8%) ( Table 2). The response curves show the quantitative relationship between environmental variables and the habitat suitability, and the responses of five variables to D. armandi suitability are shown in Figure 2. According to the response curves, annual mean temperature (Bio1) ranges from 5.6 to 17.1 °C, minimum temperature of coldest month (Bio6) ranges from −11 to 2 °C, mean temperature of driest quarter (Bio9) ranges from −5 to 7 °C, precipitation seasonality (Bio15) ranges from 63 to 94, and precipitation of driest quarter (Bio17) ranges from 12 to 74 mm. These variables were considered to have a significant effect on the survival of D. armandi.

Current Potential Distribution of D. armandi
The current potential distribution of D. armandi is located in the Qinling Mountains and Daba Mountains, including southern Shaanxi, western Henan, northwest Hubei, north Chongqing, northeast Sichuan, and southeast Gansu Provinces (Figure 3). The highly suitable habitats occur in discontinuous patches, and the moderately suitable habitats are embedded in these patches. The total suitable area amounts to 15.83 × 10 4 km 2 . The total areas for highly, moderately, and poorly suitable habitats are 4.41 × 10 4 km 2 , 4.97 × 10 4 km 2 , and 6.45 × 10 4 km 2 , respectively. Among the six provinces, Shaanxi possesses the largest suitable area, amounting to 6.73 × 10 4 km 2 .

Future Potential Distribution of D. armandi
Under future climate scenarios, the potential distribution area of D. armandi is projected to increase slightly (Figure 4a-f). The response was strongest under the RCP 8.5-

Current Potential Distribution of D. armandi
The current potential distribution of D. armandi is located in the Qinling Mountains and Daba Mountains, including southern Shaanxi, western Henan, northwest Hubei, north Chongqing, northeast Sichuan, and southeast Gansu Provinces (Figure 3). The highly suitable habitats occur in discontinuous patches, and the moderately suitable habitats are embedded in these patches. The total suitable area amounts to 15.83 × 10 4 km 2 . The total areas for highly, moderately, and poorly suitable habitats are 4.41 × 10 4 km 2 , 4.97 × 10 4 km 2 , and 6.45 × 10 4 km 2 , respectively. Among the six provinces, Shaanxi possesses the largest suitable area, amounting to 6.73 × 10 4 km 2 .

Current Potential Distribution of D. armandi
The current potential distribution of D. armandi is located in the Qinling Mountains and Daba Mountains, including southern Shaanxi, western Henan, northwest Hubei, north Chongqing, northeast Sichuan, and southeast Gansu Provinces (Figure 3). The highly suitable habitats occur in discontinuous patches, and the moderately suitable habitats are embedded in these patches. The total suitable area amounts to 15.83 × 10 4 km 2 . The total areas for highly, moderately, and poorly suitable habitats are 4.41 × 10 4 km 2 , 4.97 × 10 4 km 2 , and 6.45 × 10 4 km 2 , respectively. Among the six provinces, Shaanxi possesses the largest suitable area, amounting to 6.73 × 10 4 km 2 .

Future Potential Distribution of D. armandi
Under future climate scenarios, the potential distribution area of D. armandi is projected to increase slightly (Figure 4a-f). The response was strongest under the RCP 8.5-

Future Potential Distribution of D. armandi
Under future climate scenarios, the potential distribution area of D. armandi is projected to increase slightly (Figure 4a-f). The response was strongest under the RCP 8.5-2050s scenario, with a total suitable area of 16.75 × 10 4 km 2 , followed by the RCP 4.5-2050s    Under the RCP 2.6-2030s and -2050s scenarios, the area of highly suitable habitat is expected to increase in southern Shaanxi, and the area of poorly suitable habitat is expected to increase in southeast Gansu Provinces. Under the RCP 4.5-2030s scenario, the highly suitable area in southern Shaanxi is projected to increase continually, whereas the poorly and moderately suitable areas are projected to decrease. Under the RCP 4.5-2050s scenario, no change in the poorly suitable habitat area is projected. At the junction of Shaanxi and Hubei Province, the moderately and highly suitable areas are both projected to increase. Under the RCP 8.5 scenario, the highly suitable area is projected to increase in southern Shaanxi, northwest Hubei, and northeast Sichuan Provinces. By the 2050s, the highly suitable area is projected to increase by 0.82×10 4 km 2 . Overall, in the near future, D. armandi will gain increasingly suitable climatic space in the Qinling Mountains and Daba Mountains.

Discussion
Bark beetles comprise a large subfamily of true weevils but only a small fraction of the more than 6000 bark beetle species globally can cause landscape-scale plant mortality [39]. In central China, D. armandi destroys healthy P. armandii forests. Previous research on D. armandi has focused on integrated control. However, because the beetle spends most of its life cycle protected under the bark, measures to prevent and control outbreaks have not been fully effective [40]. Our current SDM analysis makes it possible to classify and monitor its habitat in a targeted way, based on habitat suitability maps. This both reduces the need for service staff, materials, and funding, and improves the efficiency of monitoring work. Under the RCP 2.6-2030s and -2050s scenarios, the area of highly suitable habitat is expected to increase in southern Shaanxi, and the area of poorly suitable habitat is expected to increase in southeast Gansu Provinces. Under the RCP 4.5-2030s scenario, the highly suitable area in southern Shaanxi is projected to increase continually, whereas the poorly and moderately suitable areas are projected to decrease. Under the RCP 4.5-2050s scenario, no change in the poorly suitable habitat area is projected. At the junction of Shaanxi and Hubei Province, the moderately and highly suitable areas are both projected to increase. Under the RCP 8.5 scenario, the highly suitable area is projected to increase in southern Shaanxi, northwest Hubei, and northeast Sichuan Provinces. By the 2050s, the highly suitable area is projected to increase by 0.82×10 4 km 2 . Overall, in the near future, D. armandi will gain increasingly suitable climatic space in the Qinling Mountains and Daba Mountains.

Discussion
Bark beetles comprise a large subfamily of true weevils but only a small fraction of the more than 6000 bark beetle species globally can cause landscape-scale plant mortality [39]. In central China, D. armandi destroys healthy P. armandii forests. Previous research on D. armandi has focused on integrated control. However, because the beetle spends most of its life cycle protected under the bark, measures to prevent and control outbreaks have not been fully effective [40]. Our current SDM analysis makes it possible to classify and monitor its habitat in a targeted way, based on habitat suitability maps. This both reduces the need for service staff, materials, and funding, and improves the efficiency of monitoring work.
Our findings confirmed that D. armandi is restricted to the Qinling Mountains and Daba Mountains, despite there being large numbers of hosts in northern and southwestern China. Under all the climate scenarios we studied, the highly suitable area is projected to increase somewhat in the near future. In particular, southern Shaanxi, northwest Hubei, and northeast Sichuan Provinces may become centers of distribution for this species. Although the scale of the changes in suitable area varied among the scenarios, all of the scenarios will present difficulties for forest managers in the coming decades. To address this, better monitoring and management are required in these areas. For instance, regular surveys should be conducted to ensure early detection of outbreaks, and vulnerable forest areas should be identified from maps of projected suitable distributions. Pest populations should be carefully monitored, using both visual inspection and trapping systems; this will help to determine when pest-control activity is needed. Robust forest monitoring and reporting systems should be established, to ensure timely warnings of the effects of climate change on pests and hosts, and to measure the effectiveness of management responses. These procedures are an important step in developing management strategies that integrate monitoring systems and projected climate change impacts, when conducting vulnerability and risk assessments.
In this study, temperature and precipitation were found to be the variables that constrain the current distribution of D. armandi; the minimum temperature of coldest month is particularly important. The complex relationship between temperature and physiological processes affects the geographic distribution of a species [41]. In particular, overwinter survival is a dominant factor limiting the distribution of insects. D. armandi relies on supercooling for overwintering survival. The ability to supercool is a key survival indicator for species living at low temperatures. Among its life-cycle stages, the larvae are considered the most tolerant to cold [42]. During the autumn and winter, the larvae acclimate to decreasing temperatures by increasing cryoprotectant production, thereby reducing their supercooling point. Cold tolerance is dynamic and changes throughout the life-cycle. In the coldest month, January, the supercooling point of larvae decreases to the lowest level (−7.5 ± 0.2 • C). However, it has been reported that −6.0 • C causes 100% larval mortality [43]. Observations of field populations suggested that temperatures below −10 • C reduced larvae survival [44]. In addition, for many species of bark beetle, synchronous adult emergence and life-cycle timing are required to kill the host trees [39]. Synchronous adult emergence is an important regulator of insect seasonality and synchrony, and ultimately of the mean fitness of the population. To improve their chances of surviving adverse conditions such as extreme cold or heat, vulnerable life-cycle stages must be synchronized with the appropriate seasons (the phenomenon known as seasonality). Maintaining an appropriate seasonality is critical for bark beetle population growth and outbreak potential [45]. Disrupted synchronization can severely impact insect populations, even leading to widespread mortality [46]. If the rainfall period becomes intermittent or overlong (for instance, in southern China), this may postpone the emergence of D. armandi adults.
Based on results of variable selection, host distribution has little weight as a predictor variable, and we hypothesize that at the macro-scale level the potential distribution of this beetle is limited more by temperature and precipitation variables than by the host distribution; that is, the distribution of D. armandi is not limited by the host's distribution. In fact, in many parts of the Qinling and Daba Mountains where P. armandii trees occur, D. armandi was absent. The host species is distributed far more widely than the beetle, including in northern and southwestern China. Moreover, the beetle has high host selectivity and it only attacks P. armandii older than 30 years; that is, there are enough hosts to provide for the beetles. Our findings are similar to those reported for bark beetles in North America. For example, the northward movements of Dendroctonus frontalis in United States were found to be constrained by the minimum annual temperature isotherm [47]. In western Canada, the distribution of Dendroctonus ponderosae, the mountain pine beetle, was found to be limited by the minimum winter temperature [48]. In the Sierra Madre Occidental, Mexico, the Dendroctonus rhizophagus was found to be limited by the maximum temperature isotherms [49]. However, the host distribution is dynamic in the future and the host tree distribution will respond dynamically to climate change: forest trees may persist via migration, adapt to the new conditions, or become locally extinct [50]. Some predicative results showed the host distribution in Qinling Mountains and Daba Mountains, i.e., the occurrence area of this beetle, is conservative and the change is relatively small [51,52]. Therefore, the response of D. armandi to the host may be more prominent in the local geographic range. In addition, it is worth noting that the gaps between the actual dispersal of P. armandii in the real world and that predicted by the models cannot be ignored. As a large-seeded pine species, in particular, it has a short dispersal distance that is limited to the area around the parent trees, and mainly relies on animal dispersal, which are adverse factors affecting its dispersal. Previous studies have indicated that most of the P. armandii seeds were dispersed less than 20 m [53]. In other words, the long-distance seed dispersal events are almost impossible to achieve based on the current dispersal distance, which would mean zero dispersal on the scale of predictive models. Coupled with geographical barriers, such as mountains and rivers, communications between species are almost blocked, ensuring the short-distance dispersal dominates in a local geographic range. Therefore, the research on dispersal of P. armandii seed needs to be explored. This would allow modeling of how the distribution changes of the beetle and its host tree in a local geographical range will be affected by climate change. To summarize, climate change will affect bark beetle-host interactions in complex or nonlinear ways [54,55].
Maxent models describe association between the occurrence of species populations and environmental variables [22]. Usually, tree mortality induced by bark beetle has been correlated with the climate variables describing conditions during the mortality events [56]. Such analyses are usually retrospective, rather than describing the processes leading to plant mortality. Moreover, because these relationships are likely to be different under future climatic conditions, ecosystem models must include phenology data if they are to explain physiological changes in response to a changing climate. The effects of temperature on insect physiology have been studied for decades, with the primary focus being on how temperature affects development time and survival, and how photoperiod and temperature affect diapause [57,58]. Phenological models are driven by the functional relationships between physiological process and temperature, rather than by statistical relationships [59]. Models describing bark beetle phenology require detailed information on the responses of individual beetles. Globally, there are at least 30 bark beetle species that can cause landscape-scale plant mortality [56]. Of these, sufficient data to model climate-driven phenology is available for only six species. Therefore, to model the phenology of D. armandi, we first need to collect detailed physiological information about it. Using this data, we will be able to use a phenological model to predict how temperature will impact life-cycle timing and ultimately population success in this species, based on a quantitative description of the physiological processes that are impacted by temperature.
There are still some limitations in our study. Firstly, there are the relatively low number of geo-references occurrence points in available databases. Moreover, the occurrence points collected through the field are relatively limited. Although the Maxent model has been shown to perform well with a small sample size, the sample size may still affect the accuracy of the model results [60]. Secondly, only 23 environmental variables were used to model the potential distribution of D. armandi, which did not consider the effects of biotic factors, such as dispersal, geographic barriers, competition, predation, and herbivory, that often also play roles in determining species distributions; this is clearly a critical limitation. Finally, although the individual SMD model in this study showed high prediction accuracy, there are still some limitations in prediction precision. The accuracy and performance of individual SMD varies widely among methods and species [61,62]. Some studies have shown that methods integrating multiple individual models provide robust estimates of potential species' distributions, which provide a means to improve the accuracy of models [63]. All of the factors mentioned above may cause differences between predicted distributions and actual distributions.

Conclusions
We used Maxent to map climate habitat for D. armandi, in terms of its climatic niche, based on 132 occurrence records. The minimum temperature of the coldest month and precipitation seasonality constrain its current distribution, which is within the Qinling