Predicting the Potential Global Geographical Distribution of Two Icerya Species under Climate Change

: Climate change is predicted to alter the geographic distribution of a wide variety of taxa, including insects. Icerya aegyptiaca (Douglas) and I. purchasi Maskell are two polyphagous and invasive pests in the genus Icerya Signoret (Hemiptera: Monophlebidae) and cause serious damage to many landscape and economic trees. However, the global habitats suitable for these two Icerya species are unclear. The purpose of this study is to determine the potentially suitable habitats of these two species, then to provide scientiﬁc management strategies. Using MaxEnt software, the potential risk maps of I. aegyptiaca and I. purchasi were created based on their occurrence data under di ﬀ erent climatic conditions and topology factors. The results suggested that under current climate conditions, the potentially habitable area of I. aegyptiaca would be much larger than the current distribution and there would be small changes for I. purchasi . In the future climate change scenarios, the suitable habitats of these two insect species will display an increasing trend. Africa, South America and Asia would be more suitable for I. aegyptiaca . South America, Asia and Europe would be more suitable for I. purchasi . Moreover, most of the highly habitat suitability areas of I. aegyptiaca will become concentrated in Southern Asia. The results also suggested that “min temperature of coldest month” was the most important environmental factor a ﬀ ecting the prediction models of these two insects. This research provides a theoretical reference framework for developing policies to manage and control these two invasive pests of the genus Icerya .


Introduction
Due to human activities and globalization, the speed of biologic invasion has been accelerated, leading to the homogenization of biodiversity [1]. Invasive species threaten local ecosystems by competing with native species and endangering local species and human management systems such as those related to agriculture, animal health and forestry [2][3][4][5]. Despite efforts to prevent new intrusions over the past few decades, the record data collected from the 16th century to the present day has not shown any signs of stopping. On the contrary, the invasion speed of most taxa groups increases with time. Insects illustrate this worldwide pattern of biologic invasion [6]. From the 20th century, the first, recorded rate of insect invasions has been steadily rising and remained at a high level, probably due to the global trade and the spread of host plants [7]. Global cost associated with invasive insects was estimated at US$77 billion per year equivalent to the combined cost for all goods and services and health [8].
To date, the increase in yearly invasive rates indicates the need for additional preventive measures, especially for highly invasive groups [9]. Identifying species with high transmission potential, I. purchasi were estimated based on available occurrence data using the MaxEnt software. We identified the main climatic variables that limit the potential distribution of these two scale insects and predicted the trends of their suitable habitat range under different climatic change scenarios. Our results provide a theoretical reference framework for the development of policies to manage and control these two invasive pests.

Species Occurrence Data
The occurrence data of I. aegyptiaca and I. purchasi were compiled from three sources: (1) The occurrence data recorded in the globally published papers (I. aegyptiaca [22,27,[38][39][40], I. purchasi [41][42][43]). Google Earth was used to obtain the geographic coordinates of the records when the locality names were provided in the literature; (2) The distribution data recorded by Global Biodiversity Information Facility (GBIF, https://www.gbif.org/), the Center for Agriculture and Bioscience International (CABI, https://www.cabi.org/), Taiwan Biodiversity Information Facility (http://taibif.tw/en), China Forestry Network (http://www.forestry.gov.cn/); (3) Occurrence data provided by the Forestry Department of China and on-the-spot investigation by our research group. We preliminarily screened the collected occurrence data to remove the occurrence data from greenhouse environment and enclosed buildings, only outdoor occurrence data were retained. Occurrence records tend to be in the areas that are easy to aggregate, such as near cities or other areas with high population density [44], so occurrence data can significantly affect the modeling results [45]. In order to reduce sampling deviation and eliminate spatial autocorrelation, a coarse resolution (9 km × 9 km) was created using the workflow of Li, Du & Guo [46] and a single point was randomly selected from each grid containing one or more sampling points. This part of work was completed by ArcGIS 10.2 software [47]. Overall, 46 and 167 unique presence records of I. aegyptiaca and I. purchasi were collected, respectively.

Environmental Variables
In this study, 19 bioclimatic environmental variables with effects on the growth and distribution of plants and animals were downloaded from the WorldClim Global Climate Database (http://www. worldclim.org) ( Table 1). Too many environmental variables tend to increase the dimension of the ecological space, which is not conducive to the prediction of a given model. Thus, we first used the multivariate tool in spatial analyst tools from the ArcToolbox of ArcGIS software and executed the band collection statistics command to perform the multi-collinearity analysis on each bioclimatic environmental variable layer. Then, the correlation matrix was obtained. The correlation matrix shows the values of the correlation coefficients that depict the relationship between two dataset layers and range from +1 to −1. A positive correlation indicates a direct relationship between two layers, such as when the cell values of one-layer increase, the cell values of another layer are also likely to increase. A negative correlation means that one variable inversely influences the other. A correlation of zero means that two layers are independent of each other. In this study, if the absolute value of the correlation between two variables is >0.8, then only one of the variables was selected. Second, the relative importance of each variable's contribution was assessed by sequential variable removal using Jack's knife of the MaxEnt software. In order to avoid the inaccuracy caused by the spatial autocorrelation between the training data and the testing data, we selected 25% of the occurrence data as independent data for model accuracy test. Finally, five bioclimatic environmental variables for I. aegyptiaca were selected, namely bio2, bio6, bio8, bio13 and bio18; and five bioclimatic environmental variables for I. purchasi were selected as bio5, bio6, bio8, bio13 and bio19 (Table 1). Given the uncertainties of future climate scenarios, impact assessments were incorporated with the data from a range of climate models that can effectively simulate the historical climate of various regions of the world [48]. Therefore, to predict the future potential geographical distribution of I. aegyptiaca and I. purchasi on a global scale, the community climate system model 4 (CCSM4) [49] and the model for interdisciplinary research on climate 5 (MIROC5) [50,51] for 2050 (average for 2041-2060) and 2070 (average for 2061-2080) were obtained from the global climate models (GCM) to climate projections of the International Panel on Climate Change (IPPC) (http://www.worldclim.org). For each of these two GCMs, we considered two different representative concentration pathways (RCPs) [52][53][54][55][56], which are two cumulative measures of human emissions of greenhouse gases from all sources in Watts per square meter. These pathways were developed for the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [53] and corresponded to a total anthropogenic radiative forcing of RCP = 4.5 W/m 2 (low) and RCP = 8.5 W/m 2 (high). Finally, there are eight future climate scenarios (2050RCP4.5-CCSM, 2050RCP8.5-CCSM, 2070RCP4.5-CCSM, 2070RCP8.5-MIROC, 2050RCP4.5-MIROC, 2050RCP8.5-MIROC, 2070RCP4.5-MIROC and 2070RCP8.5-MIROC) for each species' habitat prediction model. All future environmental variables were downloaded from the WorldClim Database (http://www.worldclim.org). All data of current and future bioclimate variables had 5-min (of a longitude/latitude degree) spatial resolution (this is about 9 km at the equator) [55,56].

Setting of MaxEnt Software Parameters
The MaxEnt parameter settings in this study were: the "random test percentage" of the Basic part was 25 (representing randomly selecting 75% of the distribution points as the training dataset to build the model, leaving 25% of the distribution points as the test dataset); the "replicates" was 10 for the model to run repeatedly for ten times under the same setting, "subsample" was selected as the "replicated run type", and the final result taken 10 repetitions (avg). In 2017, Phillips et al. compared the "Cloglog" output with the previous "raw", "logistic" and "cumulative" output methods, and finally proved that the "Cloglog" output was the optimal output mode for predicting the suitable area [41]. Therefore, we used "Cloglog" as the output format in this study. Because I. aegyptiaca is naturally distributed in Asia and Africa, we set MaxEnt to generate background points for I. aegyptiaca in Asia and Africa, and set the projection layers to global (excluding Antarctica); (I) purchasi is distributed in six continents, so the background points generation and the projection layers were set to global (excluding Antarctica).

GIS Analysis
We used ArcGIS software to divide the suitable area and visualize the layers. The range of probability of species' existence was set between 0 and 1, the lowest presence threshold (LPT) was used to define the suitable and unsuitable habitats [63]. The range of suitable habitats was divided into 4 categories as unsuitable habitat area (0.00-LPT), low suitability habitat area (LPT-4), moderate suitability habitat area (0.4-0.6), high suitability habitat area (0.6-1.00). A world vector map from Natural Earth Dataset (http://www.naturalearthdata.com/) was superimposed on the layer processed above. We used DIVA-GIS to combine CCSM and MIROC layers for each of the climate change scenario/year combinations (2050RCP4.5, 2050RCP8.5, 2070RCP4.5 and 2070RCP8.5) when plotting the map of the suitability habitat areas for the future climate. We also generated a potential distribution map based on "maximum training sensitivity plus specificity Cloglog threshold" using the threshold of 0.2888 and 0.2763, respectively, for I. aegyptiaca and I. purchasi. The potential distribution of I. aegyptiaca is shown in Figure S1 under the current climate, and in Figure S2

GIS Analysis
We used ArcGIS software to divide the suitable area and visualize the layers. The range of probability of species' existence was set between 0 and 1, the lowest presence threshold (LPT) was used to define the suitable and unsuitable habitats [63]. The range of suitable habitats was divided into 4 categories as unsuitable habitat area (0.00-LPT), low suitability habitat area (LPT-4), moderate suitability habitat area (0.4-0.6), high suitability habitat area (0.6-1.00). A world vector map from Natural Earth Dataset (http://www.naturalearthdata.com/) was superimposed on the layer processed above. We used DIVA-GIS to combine CCSM and MIROC layers for each of the climate change scenario/year combinations (2050RCP4.5, 2050RCP8.5, 2070RCP4.5 and 2070RCP8.5) when plotting the map of the suitability habitat areas for the future climate. We also generated a potential distribution map based on "maximum training sensitivity plus specificity Cloglog threshold" using the threshold of 0.2888 and 0.2763, respectively, for I. aegyptiaca and I. purchasi. The potential distribution of I. aegyptiaca is shown in Figure S1 under the current climate, and in Figure S2 in the future climate. Figures  S3 and S4 show the potential distribution of I. purchasi under the current climate and in the future climate, respectively.

Model Result Evaluation
The accuracy test of the MaxEnt prediction was judged by the area under the curve (AUC) of the receiver operating characteristic (ROC) curve and the partial receiver operating characteristic (pROC) metric approach [64,65]. Considering that indoor or greenhouse data may be stored at the selected occurrence point, we used 5% error rate (E = 0.05) to judge the model results, and calculated the AUC ratios (AUC ratios = AUC E /AUC 0.5 ) with E = 0.05 through niche analyst platform [65]. The evaluation criteria are: a AUC of 0.5-0.6 indicates failure model performance; a AUC of 0.6-0.7 indicates poor model performance; a AUC of 0.7-0.8 indicates general model performance; a AUC of 0.8-0.9 indicates good model performance and a AUC of 0.9-1.0 indicates very good model performance [66][67][68][69]. The AUC ratios >1 indicates that the model has a good random prediction [66].

Model Assessment
The mean AUC and AUC ratios of I. aegyptiaca are 0.883 and 1.4828. For I. purchasi, the mean AUC and AUC ratios are 0.849 and 1.3813, respectively. The LPTs are 0.0757 and 0.0507 for I. aegyptiaca and I. purchasi, respectively ( Table 2). The AUCs for all modeling results are greater than 0.8 and the AUC ratios are greater than 1, indicating that the prediction result of this experiment is particularly good and has a high credibility. Figure 2 shows the ROC curve when E = 0.05, the ROC curve under the model extends to the upper left side, indicating that the omission rate of the test sample is low. receiver operating characteristic (ROC) curve and the partial receiver operating characteristic (pROC) metric approach [64,65]. Considering that indoor or greenhouse data may be stored at the selected occurrence point, we used 5% error rate (E = 0.05) to judge the model results, and calculated the AUC ratios (AUC ratios = AUCE/AUC0.5) with E = 0.05 through niche analyst platform [65]. The evaluation criteria are: a AUC of 0.5-0.6 indicates failure model performance; a AUC of 0.6-0.7 indicates poor model performance; a AUC of 0.7-0.8 indicates general model performance; a AUC of 0.8-0.9 indicates good model performance and a AUC of 0.9-1.0 indicates very good model performance [66][67][68][69]. The AUC ratios >1 indicates that the model has a good random prediction [66].

Model Assessment
The mean AUC and AUC ratios of I. aegyptiaca are 0.883 and 1.4828. For I. purchasi, the mean AUC and AUC ratios are 0.849 and 1.3813, respectively. The LPTs are 0.0757 and 0.0507 for I. aegyptiaca and I. purchasi, respectively ( Table 2). The AUCs for all modeling results are greater than 0.8 and the AUC ratios are greater than 1, indicating that the prediction result of this experiment is particularly good and has a high credibility. Figure 2 shows the ROC curve when E = 0.05, the ROC curve under the model extends to the upper left side, indicating that the omission rate of the test sample is low.

Important Environmental Variables
Five climatic variables were selected from the previous experiments which affected the distribution of two Icerya species, and their importance was verified by Jack's knife-cutting method. The result is shown in Figure 3. The longer the blue band of a variable is, the closer between the predicted model result with the variable and the result of model prediction with all variables is, indicating that the effect of this climatic variable on the species distribution is more important. Three environmental variables of bio6, bio13 and bio18 had the greatest influence on I. aegyptiaca distribution and those of bio6, bio5 and bio19 had the greatest influence on I. purchasi distribution. The shorter the green band of a variable is, the more information and the greater the impact of this variable is on the species distribution than the other variables. The environmental variable bio6 had more specific information on the predictions of the I. aegyptiaca and I. purchasi distributions and was the most important variable to the prediction models. It can be seen that the min temperature of coldest month had the greatest influence on the distribution of these two Icerya species. more specific information on the predictions of the I. aegyptiaca and I. purchasi distributions and was the most important variable to the prediction models. It can be seen that the min temperature of coldest month had the greatest influence on the distribution of these two Icerya species.
Combining the bioclimate variables with narrow adaptive ranges to the distribution of I. aegyptiaca and I. purchasi, the partial variables were selected to plot the response curve with the existence probability: bio6 and bio13 were selected for I. aegyptiaca; bio5 and bio8 were selected for I. purchasi. Based on the response curves (Figure 4), the climatic conditions associated with high habitat suitability were >5 ℃ for bio6, 200 mm-1300 mm for bio13 for I. aegyptiaca and 25 ℃-35 ℃ for bio5, 7 ℃-28 ℃ for bio8 for I. purchasi.  Combining the bioclimate variables with narrow adaptive ranges to the distribution of I. aegyptiaca and I. purchasi, the partial variables were selected to plot the response curve with the existence probability: bio6 and bio13 were selected for I. aegyptiaca; bio5 and bio8 were selected for I. purchasi. Based on the response curves (Figure 4), the climatic conditions associated with high habitat suitability were >5 • C for bio6, 200-1300 mm for bio13 for I. aegyptiaca and 25-35 • C for bio5, 7-28 • C for bio8 for I. purchasi.

Predicted Distribution of Icerya Species
The predicted distribution of I. aegyptiaca based on the current bioclimate variables and occurrence records is shown in Figure 5. The potential distribution areas of I. aegyptiaca included a wide range of tropical and subtropical regions, mainly concentrated between 35 • S-35 • N. The main areas included southern Asia, Africa, south of North America, South America, Oceania and tropical Pacific. Based on current climatic variables, the total potential habitat suitability area for I. aegyptiaca was approximately 6533.31 × 10 4 km 2 , of which 809.14 × 10 4 km 2 had a high habitat suitability (high risk), 1102.09 × 10 4 km 2 had a moderate habitat suitability, and 4622.07 × 10 4 km 2 had a low habitat suitability ( Figure 6). The three continents with the largest habitat area were Africa, South America and Asia. The high and moderate habitat suitable areas concentrated in the south of Asia, the north of South America and the central of Africa.

Predicted Distribution of Icerya Species
The predicted distribution of I. aegyptiaca based on the current bioclimate variables and occurrence records is shown in Figure 5. The potential distribution areas of I. aegyptiaca included a wide range of tropical and subtropical regions, mainly concentrated between 35° S-35° N. The main areas included southern Asia, Africa, south of North America, South America, Oceania and tropical Pacific. Based on current climatic variables, the total potential habitat suitability area for I. aegyptiaca was approximately 6533.31 × 10 4 km 2 , of which 809.14 × 10 4 km 2 had a high habitat suitability (high risk), 1102.09 × 10 4 km 2 had a moderate habitat suitability, and 4622.07 × 10 4 km 2 had a low habitat suitability ( Figure 6). The three continents with the largest habitat area were Africa, South America and Asia. The high and moderate habitat suitable areas concentrated in the south of Asia, the north of South America and the central of Africa.

Future Invasion Risk
The MaxEnt models based on two different representative concentration pathways (RCP4.5 and RCP8.5) for the potential distribution of I. aegyptiaca in 2050 and 2070 are presented in Figure 9. The position of the potentially suitable habitat on the global scale of I. aegyptiaca was not changed much from that of the current climate, but the total suitable area was increased. The area change of the suitable habitat of I. aegyptiaca under different climate scenarios is shown in Figure 6

Future Invasion Risk
The MaxEnt models based on two different representative concentration pathways (RCP4.5 and RCP8.5) for the potential distribution of I. aegyptiaca in 2050 and 2070 are presented in Figure 9. The position of the potentially suitable habitat on the global scale of I. aegyptiaca was not changed much from that of the current climate, but the total suitable area was increased. The area change of the suitable habitat of I. aegyptiaca under different climate scenarios is shown in Figure 6 The future potential distribution of I. purchasi is presented in Figure 10. Compared with the current climate model, the total area of suitable habitats would increase in the future, but the area of suitable habitats in 2070 would become larger than that in 2050. The area change of the suitable habitat of I. purchasi under different climate scenarios is shown in Figure 8  Asia, mainly including southern China, Vietnam, Thailand, Laos, Myanmar, northeastern India, Malaysia, Indonesia, the Philippines and Ryukyu in Japan. The future potential distribution of I. purchasi is presented in Figure 10. Compared with the current climate model, the total area of suitable habitats would increase in the future, but the area of suitable habitats in 2070 would become larger than that in 2050. The area change of the suitable habitat of I. purchasi under different climate scenarios is shown in Figure 8   × 104 km2) > Current (7238.83 × 104 km2). Furthermore, the total areas of suitable habitat would increase in the future, but the areas of high habitat would decrease considerably in the future. The area of the high habitat area would reach a minimum of 869.32 × 104 km2 under the 2050 RCP85 model, a reduction of 40.30% compared to the current climate.

Discussion
In this study, we used MaxEnt niche models to estimate the global invasion risk of two important forestry Iceryas species. Based on recent recommendations published in the literature [60,64,70,71], we selected the best models for each species using different MaxEnt settings and statistical criteria. Our results support the findings by other studies and demonstrate that the model performance is significantly affected by changing feature classes and regularization multiplier values [70,71].
We used Jack's knife-cut method to measure the importance of environmental factors in the prediction models. The results showed that "min temperature of coldest month" (bio6) had the greatest impact on the distribution of these two species. The model predicted a wider tolerance range of I. purchasi than that of I. aegyptiaca under the optimum range of each bioclimate factor for the survival of the two pests. It is also predicted that I. purchasi would occupy a larger and more invasive area in the world than I. aegyptiaca. Combining with a practical example indicated that the natural distribution area of I. purchasi in China would be in the south, and the distribution area would be bounded by the line of Qinling Mountains and Huaihe River, where is the boundary between the South and the North of China. However, the distribution area of I. aegyptiaca in China would concentrate in the tropics and adjacent areas, mainly in Yunnan, Guangdong and Taiwan [27]. Therefore, in this respect, I. purchasi would occupy a larger and more invasive area in the world than I. aegyptiaca.
From the geographical distribution map of I. aegyptiaca (Figure 5), the potential distribution area concentrated between 35 • S-35 • N, which is the distribution characteristics of most insects according to latitude zones, while the distribution of I. purchasi did not have this characteristic. We can also see that Asia was the origin of I. aegyptiaca [22] and the region with the largest number of highly habitat suitability areas, suggesting Asia would be more suitable for its survival. South America, North America and Oceania had no records of the existence of the scale insects but would have a large portion of suitable areas. Once the insect invades, it may cause serious harm, so these areas should have strict inspection and quarantine to prevent the invasion of the insect. The results of future suitable habitats of I. aegyptiaca showed that environmental changes would have a greater impact on the pest distribution, and the suitable area in 2070 would change more than in 2050, and the area of suitable habitats on the global scale would increase with the drastic environmental changes ( Figure 6). Previous studies had also suggested that some invasive species were likely to experience increase due to climate change [72,73].
The global distribution point of I. purchasi indicated that the insect had successfully invaded Africa, Asia, South America, North America and Europe from its origin in Australia [20], which may be related to its polyphagism and wide environmental adaptability. The highly suitable area of I. purchasi was large, accounting for 20.11% of the total suitable area, and all six continents had highly suitable areas, while the area of the moderate suitable area and the highly suitable area accounted for 35.37% of the total suitable area, which indicates that the region suitable for the pest's survival in the world was large. The map of future suitable habitats of I. purchasi, predicts that the total suitable habitats would increase in the future, and the change in 2070 would be greater than in 2050. This shows that the changing trend of the suitable area is consistent with the change of the climate variable. Although the area of suitable habitat would increase in the future, the high suitable area could have a decreasing trend in the future. This phenomenon may indicate that the damage level by the pest may decrease with environmental changes. Previous studies have also suggested that some invasive species are likely to experience reductions due to climate change [37,72]. The retreat by invasive pests could lead to opportunities for the restoration of currently invaded landscapes, but further research is needed to identify these opportunities and to provide sound guidance for ecological restoration.
An interesting feature was found from the distribution characteristics of the suitable areas in different grades and the invasion records of these two Icerya species. The highly suitable areas had the first successful invasion record, the moderately suitable areas had the second record, and the lowly suitable areas had the third record. This phenomenon could be explained from the following viewpoints. Each species needs a certain population density to colonize in a new environment, i.e., minimum viable population [74], which is related to the fitness of the insects. The higher the fitness the insect has, the easier it is for the insect to colonize, thus fewer minimum viable population is needed for the colonization. Therefore, it is possible that these species have been carried to many new areas along with seedling transports [26], the new areas with low adaptability cannot be colonized under low population density, but with the accumulation of time and the increase of seedling transports, once the minimum viable population of these species required for colonization is reached, they will colonize the area. Thus, there is a situation in which the diffusion trend and the suitable range are consistent. This phenomenon can also be explained by propagule pressure, which is the concept of invasive biology and used to measure the population base of an alien species and its relationship with colonization success [75]. The propagule pressure affects not only the colonization of alien species, but also the stage of population growth and diffusion after colonization. When the propagule pressure is high, the population latency is relatively short, the number and spatial expansion are relatively fast, which promotes invasion. Warren et al. quantitatively described the relationship between propagule pressure and colonization. The colonization possibility of alien species was regarded as a probability event. In addition to the initial number of introduced reproductive bodies, individual reproductive viability was also considered in the analysis. The formula [76] was as follows: where P(n) is the colonization probability; B is the birth rate, D is the mortality rate, assuming b > d; R is the initial growth rate of alien species, r = b − d; n is the initial propagule number of alien species.
The above formula suggests that: (1) when the birth rate and mortality remain unchanged, the colonization probability is higher if the initial propagator number (n) of alien species is larger. Otherwise, the colonization probability is lower; (2) when the mortality rate of alien species is lower (i.e., high R value), the dependence of colonization on the number of initial propagules is smaller; (3) when the mortality rate is higher (i.e., small R value), a higher number of propagules is required to colonize successfully. In other words, it is easier for invasive species to colonize successfully in the areas with high suitability because these areas have low mortality rates and high R values, where the colonization is less dependent on the number of initial propagators. In areas with low suitability, the mortality rate of invasive species is higher, the R value is smaller, the number of reproductive bodies needed for colonization is higher, and multiple invasions are needed to colonize successfully. Thus, the tend to spread for some invasive species and the range of suitability is consistent.
The ecological niche models (ENMs) have been used in several studies to identify suitable areas for the occurrence of different species. These studies relied primarily on abiotic factors to generate models, but historical and biologic interactions as well as dispersal constraints also play an important role in species distribution [67,77]. Although these factors may limit the geographical range of a given species in its native region, in new invasion areas they may not be a constraint [78]. We consider the important influence of host plants on these two Icerya species. These two scale insects are polyphagous insects. The host plants of I. aegyptiaca are from 113 genera of 59 families [16]. The host plants of I. purchasi are from 167 genera of 68 families [16]. The distribution range of host plants is much larger than that of the suitable habitats of the two Icerya species, so we did not overlap the layers of host plant distribution with that of the suitable habitats of these two species. The existence of natural enemies will also have a greater impact on the establishment of alien insects. Previous studies suggested that the Rodolia pumila was the main natural enemy of I. aegyptiaca, but the natural distribution of R. pumila only existed in southern China, Philippines, Guam, Japan and some tropical areas of the Pacific Ocean [79]. Therefore, most of the suitable areas are at risk of being invaded by the I. aegyptiaca because they are not protected by the main natural enemies. As the main natural enemy of I. purchasi, Rodolia cardinalis has been introduced into many regions around the world to control I. purchase [80], so the natural enemy has little influence on the next invasion of I. purchasi. Due to many adaptation mechanisms of insects, different environmental factors have different strength limits for different insects. Different environmental factors act on the distribution of species at different spatial scales. On large scales, interactions between species are often being weakened, in which climate variables play a major role [81]. Therefore, our research on the prediction of the suitable area of species based on bioclimate variables can provide accurate prediction results.
The MaxEnt software has a high simulation accuracy, but it has some shortcomings; (1) There is no unified and complete evaluation system for parameter settings. Further research is needed to determine the effectiveness of parameter settings to avoid overfitting; (2) The prediction results use a probability index model with no fixed upper limit, which can increase the predicted value beyond the climatic conditions of the study area. Therefore, special attention is required when extrapolating to another research area or future (past) climatic conditions. This experiment predicts the potential distribution area on a large scale in the world, some special factors in local areas may make the two Icerya species unable to survive even when the result is expressed the local areas as the suitable area, under these special circumstances, a comprehensive analysis on a small scale is needed in subsequent experiments. Despite these limitations, MaxEnt is still an essential tool to identify suitable areas for invasive species, which ultimately represents regions more vulnerable to invasion than those with unsuitable conditions [78].

Conclusions
With the increase of global trade, logistics and transportation provide a good way of transmission for biologic invasion. Invasive species are considered as a major threat to ecosystem functions and local biodiversity. Climate change can affect the invasion process from the initial introduction to the establishment and dissemination, and thus has a profound impact on the ecosystem. I. aegyptiaca and I. purchasi are two species of polyphagous and sucking pests with strong invasiveness in the world. This study provides the first potential distribution map based on current and future climate scenarios. The results show that under current climate conditions, the potentially suitable habitat area of I. aegyptiaca will be much larger than the current suitable habitat area, and there will be little change in the habitat area of I. purchasi. In the future climate change scenarios tested, the suitable habitats of the two scale insect species will increase. This study provides a reference and a starting point for forecasting, managing and controlling two invasive species of the genus Icerya, which are seriously harmful.
Supplementary Materials: The following materials are available online at http://www.mdpi.com/1999-4907/ 11/6/684/s1. Figure S1: Potential geographical distribution of Icerya purchasi on global scale in current climate modeling base on "Maximum training sensitivity plus specificity Cloglog threshold", Figure S2: Future species distribution models of Icerya aegyptiaca on global scale under different climate scenarios base on "Maximum training sensitivity plus specificity Cloglog threshold", Figure S3: Potential geographical distribution of I. purchasi on global scale in current climate modeling base on "Maximum training sensitivity plus specificity Cloglog threshold", Figure S4: Future species distribution models of I. purchasi on global scale under different climate scenarios base on "Maximum training sensitivity plus specificity Cloglog threshold", Table S1: Correlation matrix of 19 bioclimatic environmental variables.