Mapping Invasion Potential of the Pest from Central Asia, Trypophloeus klimeschi (Coleoptera: Curculionidae: Scolytinae), in the Shelter Forests of Northwest China

Simple Summary Trypophloeus klimeschi Eggers (Coleoptera: Curculionidae: Scolytinae) causes substantial mortality to Populus alba var. pyramidalis individuals in Xinjiang, China. Currently, the number of host trees killed by this bark beetle is increasing. Climate change has exacerbated this problem, causing its range to expand. Here, using Maxent, we simulated its distributions in the 2030s and 2050s under representative concentration pathways (RCPs) 2.6, 4.5 and 8.5. The distribution expanded the most under RCP 8.5. It is large enough to pose substantial challenges for forest managers across northern China. Our study contributes to the construction and protection of shelter forests in Northern China because we present novel evidence of the potential impacts of climate change on forestry. Abstract Temperature and precipitation are the two main factors constraining the current distribution of Trypophloeus klimeschi. Currently, T. klimeschi is mainly distributed in South Xinjiang, where it occurs between the southern edge of the Tianshan Mountains and northern edge of the Tarim Basin. In addition, Dunhuang in northern Gansu also provide suitable habitats for this bark beetle. Two other potential areas for this species are in or near the cities of Alaer and Korla. Under future climate scenarios, its total suitable area is projected to increase markedly over time. Among the climate scenarios, the distribution expanded the most under the maximum greenhouse gas emission scenario (representative concentration pathway (RCP) 8.5). Jiuquan in Gansu is projected to become a suitable area in the 2030s. Subsequently, T. klimeschi is expected to enter western Inner Mongolia along the Hexi Corridor in the 2050s. In southeastern Xinjiang, however, the suitable area in northern Ruoqiang and most areas of Korla may decrease. By the 2050s, it is large enough to pose substantial challenges for forest managers across northern China. Our findings provide information that can be used to monitor T. klimeschi populations, host health, and the impact of climate change, shedding light on the effectiveness of management responses.


Introduction
Trypophloeus klimeschi Eggers (Coleoptera: Curculionidae: Scolytinae) is a newly recorded species in China [1]. T. klimeschi is native to the Kyrgyz Republic, which borders the Xinjiang Uygur Autonomous Region [2]. In 2003, this beetle was first found in Aksu, Xinjiang [1,3]. Subsequently, the beetle spread rapidly to adjacent areas and is now found in Dunhuang, Gansu Province [4]. T. klimeschi first attacks the branch shoots of Populus alba var. pyramidalis Bunge and then gradually spreads to the main trunk. The beetle completes its life cycle in the phloem under the bark, except during a brief dispersal when the adults search for new host trees [4,5]. The injured branches turn yellow and wither, maintenance [7]. All afforestation projects are flourishing with the long-term support and protection of the Chinese government. The host trees are widely distributed in northern China, especially in Xinjiang. To better understand the climatic environment, the distribution profile of the P. alba var. pyramidalis is as follows (Table 1) [7,26].

Climate Belt Environmental Overview Distribution Records
Northern warm temperate deciduous broad-leaved forest The annual average temperature is 7-12 • C, the average temperature of the coldest month is −10∼−3 • C, the average temperature of the warmest month is 18∼27 • C, and the annual precipitation is 500-700 mm.

Temperate grassland
The annual average temperature is −3-−9 • C, the accumulated temperature of ≥ 10 • C is 1600-3200 • C, and the average temperature of coldest month is −7-−29 • C. The annual precipitation is 150-500 mm, mostly below 350 mm, mainly in summer. This climate belt is distributed mainly in Xinjiang. The annual average temperature is 4-9 • C in northern Xinjiang and 7-14 • C in southern Xinjiang. The average temperature in January is −20-15 • C in northern Xinjiang and −10 • C-5 • C in southern Xinjiang. Most of the annual rainfall is below 50-100 mm, and the least is only 10-20 mm.

Occurrence Collection
The fieldwork was conducted with the assistance of the Department of Forestry Protection from May to October 2018-2019. According to the infected area, the county forestry departments were regarded as the survey units. Along the edge of the Tarim Basin, the fieldwork first started in Hotan and finally reached Dunhuang. The detailed geographical locations of occurrence area are shown in Figure 2. 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. Due to the cryptic nature of the bark beetle, traps were used to monitor beetle presence. The trap in this study consisted of six identical black plastic funnels aligned vertically over a cup, similar to those commonly used to collect bark beetles (Supplementary Figure S1). Generally, the attractant was hung in the middle and lower part of trap. Then, each trap (W × H = 18 × 50 cm, YBQ-LD-004, Geruibiyuan Technology Co., Ltd., Beijing, China) 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 away from the ground. The traps were spaced at least 30 m apart. The components of the attractant included nonanal, 2-methylbutanal, decanal, 2-hydroxybenzaldehyde, (Z)−3-hexen-1-ol benzoate, methyl, benzoate, methyl salicylate, and geraniol, based on the behavioral responses of T. klimeschi to an olfactory stimulants trial [6]. Then, eight compounds were dissolved in solvent hexane. Finally, the attractant was added into a 15-mL slow-release polyethylene vial (Sino-Czech Trading Company, Beijing, China) at a release speed of 200 mg/day. Liquid paraffin was used as a killing agent in the cups. Trapped T. klimeschi 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 data were collected by consulting the relevant literature. Those points without accurate location data were excluded. For example, Kyrgyz Republic had a few records, but lacked coordinates. As a result, there were only seven points obtained through the literature, and others were from the field survey. Finally, we obtained 89 point-based occurrence records ( Figure 1). We gathered data regarding the presence of P. alba var. pyramidalis via three methods. First, we recorded P. alba var. pyramidalis in the field from 2016 to 2018. Similarly, raw data were processed as for the T. klimeschi data. Second, we obtained P. alba var. pyramidalis presence data from the GBIF (http://www.gbif.org (accessed on 11 March 2021)) and the Chinese Virtual Herbarium (CVH, http://www.cvh.org.cn/ (accessed on 11 March 2021)). Those records without coordinates and altitudes, including duplicates, were deleted. Third, we consulted the literature for P. alba var. pyramidalis records using a Chinese academic website (https://www.cnki.net/ (accessed on 11 March 2021)). As a result, we discarded 31 occurrence records with incomplete geographic information. In total, 324 P. alba var. pyramidalis 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 T. klimeschi (Table 2). Nineteen bioclimatic variables, with a spatial resolution of 2.5 arc minutes, were freely downloaded from www.worldclim.org (accessed on 11 March 2021). Terrain variables were obtained from the United States Geological Survey (USGS, http://edcdaac.usgs.gov/gtopo30/gtopo30.html (accessed on 11 March 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 11 March 2021)) and downloaded the National Map as the analytical base map.  [28]. There has been a 41% increase in total radiative forcing since 1990. If the current trend in greenhouse gas concentrations continues, mean annual global temperature is expected to increase between 1.8 • C and 4 • C during the 21st century [29]. From the perspective of the global environment, all scenarios are likely to occur. 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 [30]. It has been widely used for climate prediction. In this study, RCP 2.6, RCP 4.5, and RCP 8.5 were selected to model the future distributions of T. klimeschi in the 2030s and 2050s.
To eliminate over-fitting of the models due to multicollinearity within environmental variables, variables with higher maximum entropy gain and non-zero regression coefficients according to the least absolute shrinkage and selection (LASSO) regularization method were selected [31,32]. LASSO is a regression analysis method that performs both the selection of variables 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 [33]. The feature selection phase occurs after 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. LASSO was applied using the R package glmnet (http://www.web. stanford.edu/~hastie/Papers/Glmnet_Vignette.pdf (accessed on 11 March 2021)).

Species Distribution Modeling
We used Maxent v.3.3.1 (http://www.cs.princeton.edu/wschapire/MaxEnt (accessed on 11 March 2021)) to predict the potential distribution of T. klimeschi. Ten-fold crossvalidation was performed to train and validate the model. The occurrence dataset was randomly divided into 10 equal-sized subsets. Of these, a single subset was retained for model-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 cross-validation process was repeated 10 times, with each of the 10 subsets used once as validation data. The 10 results were then be averaged to produce a single estimation [34]. All mathematical modeling process were based on the maxnet function in the maxnet R package (https://www.rdocumentation. org/packages/maxnet/versions/0.1.2 (accessed on 11 March 2021)). The area under the receiver operator characteristic (ROC) curve (AUC) was used to evaluate the predictive 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%; at AUC > 0.5, the model is more accurate than a random model; at AUC < 0.5, the model is less accurate than a random model. The closer AUC is to 1, the better the model performance [35,36]. All R codes used for analysis can be obtained from GitHub (https://github.com/RayLing88/Ninghang_SDM (accessed on 11 March 2021)).
The impact of sampling bias on species distribution modeling must be noted [37]. Sampling intensity varies between sites. Oversampling in some geographic spaces can cause repetition in the ecological space when building the model; this affects the simulation of the ecological needs of the species. This kind of sampling bias generally causes overfitting of niche models to species requirements, thus reducing model transferability [38]. To counter sampling bias, we randomly generated 1000 pseudo-absence points throughout China, using Arc-GIS v. 10.2 [39].
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 tenth-percentile training presence logistic threshold was adopted to define the minimum probability of a suitable habitat for T. klimeschi [40]. Based on this, habitat suitability of T. klimeschi was divided into four levels: unsuitable (0-0.1), poorly suitable (0.11-0.3), moderately suitable (0.31-0.6), and highly suitable (0.61-1). We then calculated the suitable areas under the future climate scenarios by multiplying the number of presence grid cells by their spatial resolution.

The Centroid Changes Using SDM Tool-Box
We used the SDM tool (http://www.sdmtoolbox.org (accessed on 11 March 2021)) to calculate distributional changes between two binary SDMs (i.e., the current and future SDMs). Centroid analysis summarizes the core distributional shifts in many species' ranges, and reduces each species' distribution to a single central point (a centroid); it then creates a vector file depicting the magnitude and direction of the change over time [41].

Model Performance and Variables' Response Curves
The model prediction showed a good performance, as the AUC values from the training and test datasets were 0.898 and 0.876, respectively. Based on the regression coefficients and percentage contributions of the environmental variables used in the model, the top five variables-mean temperature of coldest quarter (Bio11; −0.8217, 36.3%), precipitation of wettest month (Bio13; 0.5534, 25.1%), mean temperature of warmest quarter (Bio10; −0.4035, 17.7%), mean temperature of driest quarter (Bio9; −0.1864, 9.6%), and mean diurnal range (Bio2; −0.0569, 7.7%)-were selected as the important factors determining the habitat suitability of T. klimeschi (Table 3). Response curves show the quantitative relationship between environmental variables and the habitat suitability. According to the response curves (Figure 3), we obtained the thresholds (existence probability > 0.3) for the five variables. Mean diurnal range ranged (Bio2) from 11.7 • C to 16.4 • C, mean temperature of driest quarter (Bio9) ranged from −12.3 • C to 2.6 • C, mean temperature of warmest quarter (Bio10) ranged from 12.6 • C to 26.7 • C, mean temperature of coldest quarter (Bio11) ranged from −14.8 • C to −4.2 • C, and precipitation of wettest month (Bio13) ranged from 5 to 42 mm.

Current Potential Distribution of T. Klimeschi
The current potential distribution of T. klimeschi is in South Xinjiang and northern Gansu, amounting to 25.58 × 10 4 km 2 ( Figure 4). The total areas for highly, moderately and poorly suitable areas are 2.56 × 10 4 km 2 , 6.91 × 10 4 km 2 and 15.35 × 10 4 km 2 , respectively. As shown in Figure 2, the suitable area in South Xinjiang is distributed between the southern edge of the Tianshan Mountains and the northern edge of the Tarim Basin, amounting to 20.36 × 10 4 km 2 . In northern Gansu, this beetle is distributed in Dunhuang. Two other potential areas for this species are in or near the cities of Alaer and Korla. The highly suitable habitats occur in discontinuous patches, and the moderately suitable habitats are embedded in these patches.

Future Potential Distribution of T. Klimeschi
Under future climate scenarios, the potential distribution of T. klimeschi is projected to increase markedly over time (Figure 5a-f, Figure 6). The response was strongest under the RCP 8.5-2050s climate scenario, with the area increasing to 36.7 × 10 4 km 2 , followed by the RCP 4.5-2050s climate scenario. Compared with the current distribution, the highly suitable area in Alaer still exists in this future scenario. However, the highly suitable area in Korla will change to a moderately suitable area. Under all assumptions, Jiuquan in Gansu is projected to become a suitable area for T. klimeschi in the 2030s. Then, T. klimeschi is expected to enter western Inner Mongolia along the Hexi Corridor in the 2050s. In southeastern Xinjiang, however, the suitable area in northern Ruoqiang and most of Korla is projected to shrink.  Under the RCP 2.6-2030s scenario, the suitable area in Dunhuang is projected to expand. By the 2050s, the suitable area is projected to increase continually in Jiuquan. However, the moderately and highly suitable areas in northern Ruoqiang, and most of the area of Korla are projected to decrease persistently. Under the RCP 4.5-2030s scenario, T. klimeschi will continue to expand along the Hexi Corridor. By the 2050s, the population in Jiuquan is projected to spread to the junction of Inner Mongolia and Gansu Province. Under the RCP 8.5-2030s scenario, T. klimeschi is expected to expand to most regions of western Inner Mongolia. By the 2050s, T. klimeschi will gain the largest climatically suitable space. Overall, in the near future, T. klimeschi will gain an increasingly suitable climatic niche in northwest China.

Discussion
In northwest China, many P. alba var. pyramidalis stands are being destroyed by T. klimeschi. Previous research on this beetle has focused on integrated control. However, due to its effective concealment, measures to prevent and control outbreaks have not been fully effective. 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.
The results showed that T. klimeschi occurs between the southern edge of the Tianshan Mountains and the northern edge of the Tarim Basin. In other words, in south Xinjiang, the southern edge of the Tianshan Mountains and the northern edge of the Tarim Basin are the northernmost and southernmost boundaries of the geographical distribution of T. klimeschi, respectively. In addition, Dunhuang in northern Gansu also provides a suitable climate that supports the survival of this beetle. Furthermore, it is restricted to this area, despite there being large numbers of hosts in northern China, because the climate of this area is suitable for it. Under all of the climate scenarios that we studied, the suitable area in northern Gansu is projected to expand continuously. Along the Hexi Corridor, T. klimeschi could spread into western Inner Mongolia. Jiuquan will develop into the next center of the spread. Although the scale of the changes in suitable areas varied among the scenarios, all of the scenarios suggest 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 selected as the variables that constrain the current distribution of T. klimeschi; the mean temperature of the coldest quarter is particularly important. The complex relationships between temperature and physiological processes affect species' geographic distributions [42]. In particular, overwinter survival is a dominant factor limiting the distribution of insects. T. klimeschi relies on larvae for overwintering survival [6]. 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. It has been reported that −13 • C caused more than 90% larvae mortality [5]. Observations of field populations suggest that temperatures below −15 • C reduced larvae survival [5,6]. For many species of bark beetle, synchronous adult emergence and life-cycle timing are required to kill the host trees [43]. 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). The complex relationships between temperature and the physiological processes involved in phenology affect bark beetle population dynamics and distribution in many ways. Mechanisms promoting seasonality are critical to bark beetle population growth and outbreak potential [44]. Trypophloeus klimeschi adult infestations begin in mid-May and reach a peak from late May to mid-June, completing two generations per year [6]. During adult activity, the extreme variation of diurnal temperature and scarce rainfall in summer (for instance, in the interior of Tarim Basin) are likely to affect its seasonality. Disrupted synchronization can severely impact the growth and development of the insect, even leading to death [45]. The climatic conditions between the southern edge of the Tianshan Mountains and northern edge of the Tarim Basin, which support the presence of this bark beetle, are unique to these areas. These specific temperature and precipitation conditions constrain this species to this region.
Trypophloeus klimeschi is absent in most regions of Xinjiang, where P. alba var. pyramidalis grows. The host trees in the vast northwest of China are distributed far more widely than the beetle. The northward or southward expansion of the suitable habitat, especially in Xinjiang, did not occur in our projections. The results indicated that the distribution of the beetle is not limited to the host's distribution. Considering the mean temperature of the coldest quarter as the only limiting factor, the southward expansion of T. klimeschi is likely to occur. However, the Tarim Basin in the south Xinjiang has the largest Taklimakan Desert in China, where there is no water, and this limits the survival of the host trees. When it comes to the northern limit of the beetle distribution, the Tianshan Mountains are the natural geographical barrier prevent the beetle from spreading northward. Moreover, the mean temperature of the coldest month in northern Xinjiang (−20 • C) might not allow the overwintering of this beetle according to its cold resistance. Hence, intertwined factors restrict the southward or northward spread of the beetle. Our findings are similar to those reported for bark beetle species in North America. For example, the northward movement of Dendroctonus frontalis in the United States was constrained by the minimum annual temperature isotherm [46]. In western Canada, the distribution of Dendroctonus ponderosae, the mountain pine beetle, was constrained by minimum winter temperature [47]. In the Sierra Madre Occidental, Mexico, the Dendroctonus rhizophagus distribution was limited by the maximum temperature isotherm [48]. Furthermore, the host tree distribution will respond dynamically to climate change; forest trees may persist via migration, adapt to the new conditions, or go locally extinct [49]. All three scenarios may occur in the host's response to climate change. However, due to the fact that P. alba var. pyramidalis is the most important afforestation tree species in Northwest China, we believe that its distribution will continue to expand. The range of T. klimeschi may change to track these environmental changes. Although our results suggested alterations in the potential distribution of T. klimeschi, they did not provide absolute predictions. So far, we have not found that this bark beetle feeds on other host trees. We only considered the influence of P. alba var. pyramidalis distribution on T. klimeschi distribution. Therefore, there are still various uncertainties. For example, if the host species of the beetle become more diverse in the future climate, the current projections may prove to be conservative. However, climate change will affect bark beetle-host interactions in complex and nonlinear ways [50,51]. Therefore, future studies should consider the impacts of climate change on the host trees when studying the effects of climate change on the distribution of this bark beetle.
There are some limitations in our study. First, our occurrence data were relatively small. Although the Maxent model has shown advantages in terms of predictive performance for use with a small sample size, the sample size may still affect the accuracy of the model results [52]. Second, only 23 variables were used to model the potential distribution of T. klimeschi. The model is based on an ideal niche without considering the effects of biotic factors such as dispersal, geographic barriers, competition, predation, and herbivory, which often also play roles in determining species distributions, which is clearly a critical limitation. Finally, despite the individual SMD model in this study showing a high prediction accuracy, there are still some limitations in the precision of prediction. The accuracy and performance of individual SMDs vary widely among methods and species [53,54]. Some studies have showed that methods integrating multiple individual models provide robust estimates of potential species' distributions, providing a way to improve the accuracy of models [54,55]. All of the above factors may cause differences between the predicted distributions and actual distributions.
Maxent models describe the association between the occurrence of species populations and environmental variables [23]. We found that bark-beetle-induced tree mortality was correlated with the climate variables describing conditions during the mortality events. Such analyses are usually retrospective, rather than describing the processes leading to plant mortality [56]. 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 processes 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. Of these, sufficient data to model climatedriven phenology are available for only six species, including D. frontalis, D. ponderosa, and D. rufipennis [43]. Therefore, in order to model the phenology of T. klimeschi, we first need to collect detailed physiological information for it. Using these 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.

Conclusions
Temperature and precipitation constrain the current distribution of T. klimeschi. With climate change, the eastward expansion of its suitable area will occur. The impacts of climate change will also increase over time. By the 2050s, the total area suitable for this species under the RCP 8.5 scenario will reach 36.7 × 10 4 km 2 . More and more shelter forests will be invaded by this beetle. Therefore, pest monitoring and control measures should be taken to prevent it from spreading further eastward.