Predicting the potential distributions of the invasive cycad scale Aulacaspis yasumatsui (Hemiptera: Diaspididae) under different climate change scenarios and the implications for management

Species data

The occurrence data for CAS were compiled from three main sources: the Centre for Agriculture and Bioscience International (CABI, https://www.cabi.org/), the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/) and the existing literature (Table S1). Geo-coordinates for each chosen data-point were either referenced from information in the literature or by using Google Earth coordinates. Duplicate or unclear locations were excluded, and data were manually checked for accuracy.

Occurrence records are often biased towards areas where assemblages are easily accessible, such as near cities or other areas with high population density (Kadmon, Farber & Danin, 2004); therefore, occurrence data can significantly impacts the modelling results of the model (Elith et al., 2011). To minimize sampling bias and remove spatial autocorrelation, a coarse resolution (5 km × 5 km) was created using the workflow of Li, Du & Guo (2015) and a single point was randomly selected from each cell containing one or more sampling points (Rodriguez-Castaneda et al., 2012). The workflow was conducted in ArcGIS 10.1 (ESRI, Redlands, CA, USA. http://www.esri.com/), and after filtering, 94 locations remained, including the native range (6 points) and the invaded regions in North America (59 points), Europe (7 points), Asia (16 points) and Africa (6 points). The list of distribution sites is shown in Fig. 1 and Table S2.

Environmental variables

The environmental variables used to characterize ecological niches were selected considering the climate and topography. Climatic variables are probably the main determinant of species niches at large scale (Luoto, Virkkala & Heikkinen, 2007), and they have previously been used for insect niche modelling (Medley, 2010). All environmental variables used in this study were download from the WorldClim Global Climate Database (http://www.worldclim.org) (Hijmans et al., 2005), and current climate conditions were represented by minima, maxima and average values of monthly, quarterly and annual ambient temperatures as well as precipitation values recorded from 1950 to 2000. These parameters provided a combination of means, extremes and seasonal differences in variables known to influence species distribution (Root et al., 2003). Topographic variables represented by altitude (ALT) were derived from HYDRO1K (http://eros.usgs.gov). Multicollinearity among predictor variables may hamper the analysis of species-environment relationships (Heikkinen et al., 2006). Thus, principal component analysis (PCA) in SPSS Statistics was performed to remove highly correlated variables (r ≥ |0.85|) (Table S3). Finally, six climatic environmental and topographic variables were selected: Mean Diurnal Range (Mean of monthly (max temp–min temp)) (Bio2), Mean Temperature of Wettest Quarter (Bio8), Mean Temperature of Driest Quarter (Bio9), Precipitation Seasonality (Coefficient of Variation) (Bio15), Precipitation of Coldest Quarter (Bio19) and Altitude (ALT).

Given the uncertainty in scenarios of future climate scenarios, impact assessments should incorporate data from a range of climate models that effectively simulate historical climate over the area of interest (Guisan et al., 2013). Therefore, to predict the future potential distribution of CAS on a global scale, the Hadley Global Environment Model 2-Atmosphere Ocean (HADGEM2-AO), Geophysical Fluid Dynamics Laboratory Coupled Model v3 (GFDL-CM3) and Model for Interdisciplinary Research on Climate (MIROC5) for 2050 (average for 2041–2060) and 2070 (average for 2061–2080) were obtained from the Coupled Model Intercomparison Project Phase 5 (CMIP5) of the fifth assessment of the International Panel on Climate Change (IPCC) in 2010 (Moss et al., 2010). Four possible greenhouse gas concentration trajectories were selected that ranged from low (RCP 2.6) to high (RCP 8.5) and corresponded to increases in global radiative forcing values in the year 2100 relative to preindustrial values (2.6, 4.5, 6.0 and 8.5 w/m², respectively). Additionally, it was assumed that the global mean temperature will increase by 0.3 °C–4.8 °C by 2100 across all RCPs (Van Vuuren et al., 2011). Four climate change scenario/year combinations were selected to simulate the future influence of potentially extreme climatic variations on the geographical distribution patterns: RCP 2.6-2050, RCP 8.5-2050, RCP 2.6-2070 and RCP 8.5-2070. All future climate projection data were download from the WorldClim Database (http://www.worldclim.org). All current and future data had a spatial resolution of 2.5 arc min (approx. ∼5 km resolution at the equator).

MaxEnt modelling

In recent years, MaxEnt software (MaxEnt version 3.3.3k, http://www.cs.princeton.edu/s˜chapire/maxent/) has frequently been used to simulate shifts in species ranges under current and future climate scenarios (Ponce-Reyes et al., 2012; Wong et al., 2013; Cooper et al., 2016). This software is popular because it requires only presence-only data. Comparative studies have consistently shown that MaxEnt has excellent performance and outperforms than many other methods (such as GARP) in estimating potential species distributions, particularly when sample sizes are small (Elith et al., 2006; Bosso et al., 2016; Sultana et al., 2017). Moreover, some studies found that this software performed well regardless of the number or geographical extent of species records, compared to the performance of mahalannobis typicalities, random forests methods and other methods (Bosso et al., 2016; Hernandez et al., 2006; Elith et al., 2011).

Building models from MaxEnt with an appropriate amount of complexity for the study objectives is critical for robust inference (Merow et al., 2014). Moreover, for several reasons, the default settings of MaxEnt are expected to produce overfit models (Radosavijevic & Anderson, 2014). Recent studies have demonstrated the importance of both balancing model complexity and predictive power and evaluating model’s performance with spatially independent data. Hence, in order to produce the best possible model for CAS, avoiding overfitting while maximizing predictive power, we employed the R package ENMeval to select the optimal combination of two important MaxEnt’s parameters, the value of the regularization multiplier and the combination of feature classes (Muscarella et al., 2014; Warren & Seifert, 2011). In this study, the “checkerboard2” approach was used by calculating the standardized Akaike information criterion coefficient (AICc), and the parameterizations that resulted in the model with the lowest delta AICc score were selected to run the final MaxEnt models. The regularization multiplier was varied from 0.5 to 4 in increments of 0.5, and the following six feature combinations were tested: (1) Linear (L); (2) Linear (L) and Quadratic (Q); (3) Hinge (H); (4) Linear (L), Quadratic (Q) and Hinge (H); (5) Linear (L), Quadratic (Q), Hinge (H) and Product (P); (6) Linear (L), Quadratic(Q), Hinge (H), Product (P) and Threshold (T). The ENMeval package were implemented in R 3.1.3 (R Core Team, 2015). The result were shown in Table S4 and Fig. S1. Thus, a regularization multiplier = 2.5; feature combinations = Linear, Quadratic and Hinge (LQH) were choose in final MaxEnt configuration. Background points were randomly chosen within the area enclosed by a minimum convex polygon comprising all records of the species of interest as suggested by Phillips (2008). The logistic output of MaxEnt was used for all the analyses. The 10th percentile training presence logistic threshold was used to define the suitable and unsuitable habitats for CAS. This threshold is a conservative value that is widely adopted in species distribution modelling, especially when data have been collected by different observers and have been collected using different methods over a long time (Bosso et al., 2016). By applying this threshold, binary presence/absence maps were created using the reclassify module from ArcGIS 10.1 (ESRI, Redlands, CA, USA). To maximize the predictive information and simplification of future analysis, the suitable habitat areas for CAS were reclassified into four levels: unsuitable habitat (no risk), low habitat suitability (low risk), moderate habitat suitability (medium risk) and high habitat suitability (high risk). The area under the curve (AUC) of the receiver operating characteristic (ROC) (Fielding & Bell, 1997) and the true skill statistic (TSS; sensitivity + specificity − 1) (Allouche, Tsoar & Kadmon, 2006) were used to estimate the performance of the model according to Swets (1988). The AUC value ranges from 0 to 1, where a value below 0.5 can be interpreted as a random prediction; 0.5–0.7 indicates poor model performance; 0.7–0.9 indicates moderate performance; and a value above 0.9 is considered to have “good” discrimination abilities (Peterson et al., 2011). The TSS accounts for both omission and commission errors, and the value range from −1 to 1, where 1 indicates perfect agreement and 0 represents a random fit. Excellent model performances are expressed by TSS values close to 1 (Allouche, Tsoar & Kadmon, 2006).

A 10-fold cross-validation was used to run MaxEnt to prevent random errors from affecting the selection of the validation and prediction sample. To assess the influence of environmental variables on species, the jackknife test was used to measure the importance and the percent contributions of each variable. To account for variation among GCMs, a final projection environmental suitability map was created for CAS by averaging the projections from the all future climate scenarios.

Model performance for potential distribution

The model performance for A. yasumatsui was better than random based on 10-fold cross validation, with a mean AUC value of 0.939 ± 0.02 and a TSS value of 0.701 ± 0.11. A “maximum training sensitivity plus specificity” threshold value of 0.112 was obtained from the 10th percentile training presence of the species occurrences. The suitable habitat areas for the pest were reclassified into four levels: <0.112 which indicated unsuitable habitat (no risk); 0.112–0.4, which indicated low habitat suitability (low risk); 0.4–0.6, which indicated moderate habitat suitability (medium risk); and >0.6, which indicated high habitat suitability (high risk).

Important environmental variables

Among the six environmental variables, Mean Temperature of Driest Quarter (Bio9) and Altitude (Alt) had the largest contributions to the distribution model for A. yasumatsui (Table 1). These two factors could explain 74.8% of the modeled distribution. The contributions of the other factors, i.e., Precipitation of Coldest Quarter (Bio19), Mean Temperature of Wettest Quarter (Bio8), Mean Diurnal Range (Bio2) and Precipitation Seasonality (Coefficient of Variation) (Bio15) had contributions of the model were 15.2%, 6.8%, 2.7% and 0.5%, respectively, to the model. It appeared thermal conditions and topography were more important than other variables in creating the distribution map for this pest.

Current invasion pattern

The current distribution pattern of A. yasumatsui based on the current environmental variables and occurrence records (including invasive records and native records) is shown in Fig. 1. The current climatic variables in the invasive range of the pest, as predict from MaxEnt, were used to identify potential distribution areas around the world, which were especially concentrated in North America, South America, Europe and Asia. In the USA, Florida, Georgia, Alabama, Louisiana and Texas had high invasion risks under the current climate scenario and areas with moderate invasion probabilities were also found in eastern Mexico. In addition, the Hawaiian Islands were also identified as having high potential distribution areas on the current map. In Europe, moderate-risk areas included France and Italy. In Asia, India, South China, Malaysia, Vietnam and Japan as well as portions of eastern Laos and Indonesia were the areas with the highest risk on the current potential distribution map. In Africa, the potential distribution area will extend to the east of Africa. Based on current climatic variables, the total area of potential suitable habitat for this pest is approximately 46,432,375 km², of which 2,390,775 km² (i.e., 5.1% of the total potentially invadable area) has high habitat suitability (high risk), and 5,557,900 km² (approximately 11.9% of the total potentially invadable area) has moderate habitat suitability (Table 2).

Based on the response curves (Figs. 2A–2F), the climatic conditions associated with high habitat suitability were 25–28 °C for Bio8, 15–20 °C for Bio9, 200–400 mm for Bio19. However, both Bio2 and Bio15 shown no relation to high habitat suitability.

Future invasion risk

The MaxEnt models based on the RCP 2.6 emission scenarios for the potential distribution of A. yasumatsui in 2050 are presented in Fig. 3A and Table 2. The potential area of invasion with suitable habitat was 40,077,550 km², which represents an shrink of 13.6% over the current area of suitable habitat. The area of highly suitable habitat expand to 2,734,675 km² (an increase of approximately 14.3% over the current highly suitable area).

The potential area of distribution with all suitable habitat will be decrease, but area of highly suitable habitat expand in some area during this period. In Asia (Fig. S2 and Table 2), the model predicts that the amount of suitable habitat in this area will be 10,441,600 km² based on current climatic variables. Suitable conditions would expand under the RCP 2.6-2050 future climate scenario (Fig. S2 and Table 2), but the areas with high potential invasive ranges were mostly concentrated at low latitudes, where the species will continue to expand to an area of 1,350,900 km² (an increase of approximately 88.9% over the current highly suitable area) under RCP 8.5-2050. Expansions in the total suitable area and highly suitable area also occurred in Europe (Fig. S3 and Table 2). However, the suitable and highly suitable areas decreased in Africa (Fig. S4 and Table 2) and in North America (Fig. S5 and Table 2) under these climatic scenarios.

Under the RCP 2.6-2070 future climate scenario, MaxEnt predicted a gain in the area of suitable habitat of 39,006,250 km², approximately 15.9% decrease over the current status (Fig. 3B and Table 2). The highly suitable area might notable increase by 17.1%% under this scenario, expanding to 2,801,300 km². The restricted range expansion of the suitable habitat area under this climatic scenario also occurred in Italy, the Kingdom of Bhutan, north-west Africa, northern Bangladesh, Southeast India, South China and some regions in Russia (Fig. 3 and Table 2). Moreover, the expansions in the total and highly suitable areas also occurred in Europe (Fig. S3 and Table 2) and Asia (Fig. S2 and Table 2). However, the suitable and highly suitable areas decreased in Africa (Fig. S4 and Table 2) and in North America (Fig. S5 and Table 2) under these climatic scenarios. Under RCP 8.5-2050, the model-predicted area of suitable habitat was 38,217,200 km², which represents an decrease of approximately 15.9% over the current suitable habitat area (Fig. 3C and Table 2), and the highly suitable area expanded to 2,801,300 km² (an increase of 17.1% over the current extents). In this scenario, the potential highly potential suitable area in Asia will be increase to 1,369,925 km² (approximately 91.5% increase over the current area) (Fig. S2 and Table 2) and 129,925 km² in Europe (approximately 41.4% increase over the current area) (Fig. S3 and Table 2). However, in Africa (Fig. S4 and Table 2) and North America (Fig. S5 and Table 2), the highly suitable area will shrink to 54,625 km² (i.e., a decrease of about 60% compared to the current extent) and 605,675 km² (i.e., a decrease of about 1.5% compared to the current extent), respectively.

MaxEnt predicted that the area of suitable habitat under the RCP 8.5-2070 future climate scenario would be 38,189,200 km², which is actually an decrease approximately 17.7% over the current area (Fig. 3D and Table 2), but the highly suitable area expanded to 2,859,175 km² (a 19.5% increase). In this scenario (Fig. S2 and Table 2), the potential highly suitable area in Asia will continue to expand to 9,091,925 km² (an increase of approximately 11.9% over the current extent); in Africa (Fig. S4 and Table 2), the area will shrink to 37,800 km² (a decrease of approximately 72.3%); and in Europe (Fig. S3 and Table 2), the area will expand to 141,775 km² (an increase of approximately 54% over the current extent). However, the potential highly suitable area in North America (Fig. S5 and Table 2) will further decrease to 501,250 km².

In short, the potential distribution of this pest will increase with climate change. Specifically, the area of highly suitable habitat will increase in most areas, but decreases in both the areas of suitable and highly suitable habitats will occur in Africa (Fig. S4 and Table 2) and also in North America (Fig. S5 and Table 2) under the future climatic scenario.

Change in habitat suitability

In this paper, the potential distribution of the invasive CAS was first simulated using MaxEnt and presence data. The CAS is native to Thailand and has already invaded other continents. The results of the simulation of the potential distribution of this invasive scale insect showed that climate change would affect its worldwide distribution. Our model also suggests that the CAS has a high invasive potential and is capable of invading large portions of the global landmass on several continents, particularly in South China, North India, north-western Africa and some areas of Russia. The occurrence of unsuitable areas on the current distribution map of A. yasumatsui suggests that these areas are constrained by inappropriate ecological conditions due to natural enemies or local adaptation, but these areas are expected to change over time with the increase in climate warming. This pest is already causing serious damage to the mature leaves and trunks of many cycad species (Normark et al., 2017) in Southern China, Guam, and Florida, where the pest is invasive. Interestingly, our results suggests that western Australia, northern New Zealand, western Madagascar, southern Brazil, Uruguay and northeastern Argentina have suitable habitat for this pest; however, fortunately, there have not been any reports of its presence at these sites. Thus, strict quarantine measures should be enforced across these regions to prohibit the spread of this pest from its invasive and native ranges. The potential distribution map for this pest provides a reference for developing monitoring strategies to detect future infestations in currently uninfected regions.

However, not all regions showed an expansion of the potential distribution range under climate change. Based on two climate scenarios, the current research showed that climate change will reduce the risk of invasion on some continents, such as Africa and previous studies have also suggested that some invasive species are likely to experience range reductions due to climate change (Sarma, Munsi & Ananthram, 2015; Bradley et al., 2010; Roura-Pascual et al., 2004). 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 provide sound guidance for ecological restoration.

Factors influencing invasion

The results of the potential CAS distribution modelling by MaxEnt showed that the Mean Temperature of Driest Quarter is the most important climatic variable defining the current global distribution of this pest. Previous results have demonstrated that CAS eggs hatch to nymphs in 8–12 days when the temperature is 24.5 °C (Howard et al., 1999) and adult female would be not laid eggs when temperature at 18 °C and 35 °C (Cave, Sciacchetano & Diaz, 2009). Moreover, 25 °C–28 °C has been confirmed to be the most suitable temperature range for the development and reproduction of CAS (Ravuiwasa et al., 2012). However, temperatures below 15 °C have been shown to be unsuitable for this insect, as indicated by the high mortality observed in laboratory research (Teaaro, 2009). Moreover, some research has suggested there were significant correlations between population growth and temperature variations. Our results confirmed that the mean temperature of driest quarter is the most important factor affecting the model, accounting for 48.1% of the contribution of all climatic variables. Our model suggested that the population of this pest would increase under an mean temperature of 15–20 °C in the driest quarter and a mean temperature of 25–28 °C in the wettest quarter. Therefore, these temperature values represent the most suitable climatic conditions for the population growth of this pest and some currently unsuitable habitats, such as southern China, may be come suitable in the future as currently suitable areas become too hot for CAS development. In contrast, some currently highly suitable habitats, such as North America, may become unsuitable due to global warming. This decrease in habitat suitability in these regions may constraint the spread of CAS outside its current suitable range or limit its development.

Although the potential distribution of CAS under climate change was predicted by MaxEnt, other factors may have limited the accuracy of our results. Only constraints caused by climate conditions were considered in the current research; however, other factors, such as host-plant availability (Ning, Wei & Feng, 2017), interspecific interactions (Gao & Reitz, 2017), and dispersal ability (Guisan & Thuiller, 2005), also affected the precision of the model; these factors should considered in future research.

Implications for management strategies

Reasonable management strategies and effective control methods are the key to controlling the spread of CAS and they should be based on a broad consideration of the relationship between invasive species and other factors such as climate change, land use, and human activities (Pyke et al., 2008). Based on the results of the current and past studies, one possible method and management strategy for controlling this species is proposed here.

First, identifying the areas experiencing climate change is a critical step (Pyke et al., 2008). Species distribution modeling (SDM) is a cost-effective, simply early warning system that enables the identification of potential invasion areas (Warren, 2012). Forecasting future areas of invasion by SDM provides an opportunity for governments to prioritize those regions for management and then develop strategies to control pests in these areas. For example, our results suggest that CAS has a highly suitable environment in southern China, which has severely harmed Cycas revoluta in Shenzhen City of Guangdong Province (Yang et al., 2009) and Cycas guizhouensis in Guizhou Province (Wu, Li & Luo, 2008); however, unfortunately, no management strategy for controlling this pest is exists in these region. Moreover, some regions, such as western Australia, northern New Zealand, western Madagascar, southern Brazil, Uruguay and north-eastern Argentina, were identified as having suitable habitat for the CAS in our work. Therefore, these areas should formulate strict quarantine measures to prevent invasion by this pest.

Second, many studies have shown that the use of oil (organocide or paraffin-based ultra-fine horticultural oils) on foliage can significantly reduce CAS population size, but several applications may be required (Hodges, Howard & Buss, 2003). Conventional insecticides, such as cygon and pyriproxifen (Weissling, Howard & Hamon, 1999; Emshousen, Mannion & Glenn, 2004; Howard et al., 1999), can also control CAS in some instance. However, the chemical control efficiency is low, because the scale can protect itself by hiding between the leaves or stems of cycads (Hodges, Howard & Buss, 2003). Therefore, a combination of biological control agents and chemical control can improve management efforts.

Finally, natural enemies are likely another effective way to control CAS. The first biological control agents, a parasitic wasp, Coccobius fulvus and a predatory beetle, Cybocephalus nipponicus, were released in 1998 (Cave, 2006), but a satisfactory level of control has not been achieved. Since then, many natural enemies have been used to control CAS at certain times and over certain ranges; for example, both Cybocephalus nipponicus and C. flavocapitis are effective biological control agents of CAS in Taiwan Province, China (Bailey, Chang & Lai, 2011). Approximately 15,000 Coccobius fulvus individuals were released in 14 Florida counties during February–April 2002, and good CAS control results were observed (Cave, 2006). Moreover, some new natural enemies have been found to feed on CAS, but their control efficacy needs to be further verified. These biological control agents included Rhyzobius lophanthae, Aphytis lingnanensis and Aphytis lepidosaphes.

In brief, areas that have been invaded by CAS must combine chemical control with natural predators to manage this pest. Additionally, strict quarantine measures should be implemented in areas where this pest has not yet been recorded but where climatic conditions appear to be suitable to CAS, as suggested by our results.