Suitable Habitats for Cicadella viridis and Evacanthus interruptus (Hemiptera: Cicadellidae) with Global Climate Change

Climate change refers to changes in global and regional climate over time and is an important concept in the study of organism distribution, especially among insects (Deutsch et al. 2008, Liu and Shi 2020). Climate change affects growth, development, survival, migration, and spread of pests (Wang et al. 2020a). Thus, it causes significant changes in the distribution patterns of pest habitat suitability (Zhu et al. 2017, Fan et al. 2020, Xu et al. 2020) and promotes instability in agricultural and forest production. Identifying changes in potential habitat suitability areas of pests resulting from climate change is crucial to improving the effectiveness of pest management (Wei et al. 2019, Wang et al. 2020b, Liu and Shi 2020).

Cicadella viridis and Evacanthus interruptus L. (Hemiptera: Cicadellidae) are two important leafhopper pests in agricultural and/or forestry production. Both are mainly distributed in Europe and East Asia. Cicadella viridis harms many host plants, including willow and poplar trees, corn, rice, soybean, and potatoes, whereas E. interruptus damages a relatively small number of host plants, primarily corn and wheat (Li 1991, Jin et al. 2007, Yang et al. 2017). Both insects damage host plants by feeding on plant tissue and by lacerating plant structures during oviposition and in severe cases cause the death of the host plant. Additionally, C. viridis also indirectly harms host plants by transmitting plant diseases (Yang et al. 2017). Economic losses caused by C. viridis reach up to 200,000 yuan/yr in Altay (about 11,500 km²), China (Jin et al. 2007).

Identification of areas in which these two pests are distributed is needed to assist in formulating effective management measures. However, to date, the current and future habitat suitability areas of these pests are unclear, resulting in reduced efficiency in management or eradication efforts.

Species distribution models (SDMs) relate environmental variables to species occurrence records to forecast the current and future habitat suitability areas (Elith and Leathwick 2009). These models provide a tool to facilitate pest management. For example, Wei et al. (2018) identified potential habitat suitability areas of the scale insect Aulacaspis yasumatsui Takagi (Homoptera: Diaspididae) under different climate change scenarios using SDMs. These results were used to provide a theoretical reference framework for developing policies to manage this pest. During the past few decades, many SDMs have been developed, including maximum entropy (MaxEnt) (Phillips et al. 2006), domain environmental envelope (DOMAIN) (Carpenter et al. 1993), and bioclimatic modeling (BIOCLIM) (Beaumont et al. 2005). The emergence of many SDMs provides researchers with more options; however, different models yield different results. Thus, the choice of model is crucial. The MaxEnt model is the most commonly used SDM because it has the advantage of high prediction accuracy, uses presence data, and performs better with low numbers of occurrence records compared with the other models (Hernandez et al. 2006, Pearson et al. 2006, Phillips et al. 2006). MaxEnt forecasts the distribution range of a species by finding the distribution that has maximum entropy subject to constraints derived from environmental variables found at the occurrence sites (Phillips et al. 2017).

In this study, the MaxEnt model was used to forecast the current and future habitat suitability areas for C. viridis and E. interruptus based on their occurrence records and related environmental variables. This information was then used to design efficient strategies for controlling them in the future. We focused on identifying (a) the primary environmental factors affecting the distribution of the two species and (b) the current and future habitat suitability areas of the two species.

Species occurrence data. Global occurrence records of C. viridis and E. interruptus were compiled from the large reference datasets of the Global Biodiversity Information Facility (GBIF; http://www.gbif.org/). The occurrence records without coordinates were georeferenced using Google Earth. Notably, sampling bias may give rise to strong inaccuracy in the resulting model and lead to incorrect predictions (Fourcade et al. 2014); thus, we used the spatial filtering method to reduce possible sampling bias. Specifically, a grid cell with 0.5° × 0.5°spatial resolution was created in ArcGIS 10.7, and a random occurrence was chosen from each grid cell containing one or more occurrence records. After filtering the data, we obtained 1,130 occurrence records for C. viridis and 567 occurrence records for E. interruptus. The global occurrence records of the two species are visualized in Fig.1.

Environmental variables. The 19 bioclimatic variables with a 2.5-arcmin (approximately 5-km resolution at the equator) spatial resolution were obtained from WorldClim version 2.1 (http://www.worldclim.org). To reduce multicollinearity among environmental variables, Spearman's correlation was performed using SPSS 22 software (Zhang et al. 2019, Beeman et al. 2021), and the highly correlated environmental variables (r > ∣0.85∣) were removed. Finally, six environmental variables remained for C. viridis, which were mean diurnal range, isothermality, mean temperature of the wettest quarter of the year, precipitation seasonality, precipitation of the warmest quarter of the year, and precipitation of the coldest quarter of the year. Five environmental variables remained for E. interruptus, as follows: mean diurnal range, isothermality, mean temperature of the wettest quarter of the year, mean temperature of the driest quarter of the year, and precipitation seasonality. Furthermore, we used computed variance inflation factors (VIFs) for selected environmental variables for two species, which were determined to be <10. Thus, there was no significant multicollinearity in the selected variables.

We used the future bioclimatic variables in 2041–2060 and 2061–2080, with two shared socioeconomic pathways (ssp126 and ssp585). Given the uncertainty of future climate models, the data of future bioclimatic variables were represented by the mean value from eight climate models (BCC–CSM2–MR, CNRM–CM6–1, CNRM–ESM2–1, GFDL–ESM4, IPSL–CM6A–LR, MIROC–ES2L, MIROC6, and MRI–ESM2–08).

Modeling procedure and habitat suitability areas. MaxEnt modeling was conducted using MaxEnt 3.4.4 software. Occurrence records were partitioned into two parts, with 75% used to setup the model and the remaining 25% for model testing. The hinge features were chosen to improve model performance in the analysis process (Phillips and Dudík 2008), and regularization multiplier parameter was set to 1, which is a common value used previously in other studies (e.g., Wang et al. 2020a, Liang et al. 2021). The background points with 10,000 were selected. Additionally, subsampling was used to generate a stable model because it is superior to cross-validation (Anderson and Raza 2010) and bootstrapping (Rospleszcz et al. 2014). Next, five replications were selected to run the model. These previously mentioned parameters have been widely used in modelling efforts. The area under curve (AUC) of the receiver operating characteristics was used to evaluate the accuracy of the model. Then, habitat suitability areas were classified into 4 levels based on the occupancy probability (Xu et al. 2020, Li et al. 2020), with 0–0.2 indicating unsuitable habitat area, 0.2–0.4 indicating low habitat suitability area, 0.4–0.6 indicating moderate habitat suitability area, and 0.6–1.0 indicating high habitat suitability area.

Model performance and important environmental variables. The AUC mean value obtained from the MaxEnt model was 0.938 and 0.962 for C. viridis and E. interruptus (Fig. 2), respectively. The model results showed that isothermality was the most significant environmental variable for the distribution of C. viridis, followed by precipitation levels during the coldest quarter of the year (Table 1). Additionally, the most important and second most important environmental factors affecting the distribution of E. interruptus were mean temperature during the driest quarter of the year and mean diurnal range (Table 1). The responses of environment variables to C. viridis and E. interruptus are illustrated in Figs. 3 and 4, respectively.

Current habitat suitability areas. The habitat suitability areas of C. viridis were concentrated mainly in European countries, such as Poland, Switzerland, and Russia, and in many provinces of southern China, such as Yunnan, Guangdong, Guizhou, and Jiangxi (Fig. 5). Additionally, the northeastern fringes of the United States, the southwestern and southeastern fringes of Canada, and Japan are additional locations of habitat suitability.

The distribution range of habitat suitability areas for E. interruptus is narrower than that for C. viridis, with the primary areas in Europe and southern China (Fig. 5). Western Europe was the main area of distribution of the two species high habitat suitability area, but C. viridis had a wider distribution range than E. interruptus (Fig. 5). Additionally, the highly suitable areas of C. viridis were also concentrated in southern China. The total habitat suitability area for C. viridis was 1.62 × 10⁷ km², of which 5.46 × 10⁶ km² showed high habitat suitability, 4.2 × 10⁶ km² showed moderate habitat suitability, and 6.49 × 10⁶ km² showed low habitat suitability (Table 2). The total habitat suitability area (9.46 × 10⁶ km²) of E. interruptus was smaller than that of C. viridis. Among the three levels of habitat suitability areas, the high habitat suitability area was 2.92 × 10⁶ km² (mainly in western Europe), the moderate habitat suitability area was 2.44 × 10⁶ km², and the low habitat suitability area was 4.28 × 10⁶ km² (Table 2).

Future habitat suitability areas. Future habitat suitability areas in 2041–2060 and 2061–2080 were predicted using the MaxEnt model (Figs. 6 and 7). The results showed that the locations of future habitat suitability for the two species had not changed much compared to that under current climatic conditions, but areas of habitat suitability had increased. Furthermore, the area of high habitat suitability had decreased, increased or was unchanged, whereas the areas of moderate and low habitat suitability had increased (Table 2).

The total habitat suitability area of C. viridis increased the most under climate scenario ssp585 during 2041–2060 and 2061–2080, which is an increase of 2.46% from the current extent, relative to the total habitat suitability areas in the current climate (Table 2). The high habitat suitability areas of this species decreased under the 4 climate scenarios, with 2.01%, 14.83%, 3.29%, and 19.78% decreases in 2041–2060 ssp126, 2041–2060 ssp 585, 2041–2060 ssp126, and 2061–2080 ssp585, respectively (Table 2). This pattern is concentrated in western Europe and southern China (Fig. 6).

Under future climate scenarios, total habitat suitability expansion areas for E. interruptus were primarily located on the east and west coasts of the United States and Canada (Fig. 7). Moreover, the area of high habitat suitability expanded in the 2041–2060 ssp585 scenario but decreased in 2041–2060 ssp126 and 2061–2080 ssp126 scenarios (Table 2). The area of both low and moderate habitat suitability under the four future climate scenarios is larger than that under current climatic conditions (Table 2), and reaches its maximum in 2061–2080 ssp126 (25.7% area increase) and 2061–2080 ssp 585 (12.29% area increase), respectively (Table 2).

In this study, MaxEnt modeling was used to predict the current and future global habitat suitability areas for C. viridis and E. interruptus based on occurrence records and environmental variables. The future climate data used in the predictions were obtained using different models, which adds uncertainty. Thus, environmental data from a single model may bias predictions. To avoid this problem, the mean of the eight climate models was calculated as the “true” value. Recently, many studies have used future climate data from the Coupled Model Intercomparison Project Phase 5 (CMIP5), to conduct species distribution modeling to predict future habitat suitability areas for species. However, recently, future climate data provided by the Coupled Model Intercomparison Project Phase 6 (CMIP6) have been released. The primary difference between CMIP5 and CMIP6 involves future scenarios (Kamruzzaman et al. 2021). The ssp of CMIP6 consider future changes in socioeconomic conditions (e.g., populations and ecosystems) to assess emission scenarios and include climate change mitigation and adaptation efforts (O'Neill et al. 2016), which are considered more realistic future scenarios than the concentration pathways of CMIP5 (Song et al. 2021). Therefore, future habitat suitability areas based on future climate data from CMIP6 are relatively more accurate. Furthermore, in this study, ssp126 and ssp585 were adopted to obtain the maximum difference of habitat suitability areas for the two pests under future climate scenarios.

According to the habitat suitability area maps of the two species, the distribution range of suitable habitat predicted by the MaxEnt model is larger than the distribution range of currently known occurrence records and is mainly concentrated in Europe and southern China (Fig. 6). Clearly, some regions with no occurrences are also suitable survival locations for the two species, such as in the United States and Canada, indicating that these areas should be monitored for accidental introductions. Currently, many parts of the world do not have suitable habitats for the two species (e.g., Africa), which implies that these areas do not have suitable environmental conditions for the survival of these species. As shown in previous studies (Biber-Freudenberger et al. 2016, Huang et al. 2019, Wei et al. 2020, Wang et al. 2020b, Lee et al. 2021), climate change will lead to an increase or decrease in the range of suitable habitat areas for pests, which is related to the ability of species to adapt to climate change. This study found that the areas of suitable habitats for both pests under all future climate scenarios were larger than those under current climatic conditions. However, under climate change, the increase in the area of the suitable habitats is not large (Table 2), and the distribution range of the suitable habitats is basically the same as that of the current range (Figs. 5, 6, and 7). The results reflect that the two species will hardly spread to other areas in the future.

We also found that the high habitat suitability areas of the two species showed different changes under future climate changes. Compared with the current climate condition, the high habitat suitability area of C. viridis will decrease under future climate changes, whereas the high habitat suitability area of E. interruptus will increase in future climate scenario ssp585 in 2041–2060 and in 2061–2080 and decrease in future climate scenario ssp126 in 2041–2060 and in 2061–2080. More attention is needed to monitor and control the two species in high habitat suitability areas.

The model results showed that the most important environmental factors affecting the distribution of the two species were isothermality and mean temperature during the driest quarter of the year, of which both were temperature-related variables. These results are consistent with previous studies involving other insects, such as scale insects (Wei et al. 2017, Zhao et al. 2019), Diaphorina citri (Kuwayama) (Wang et al. 2020a), beetles (Marchioroa and Krechemer 2018, Tang et al. 2019), leafhoppers (Santana et al. 2019), Opisina arenosella Walker (Zhao et al. 2015), and Spodoptera frugiperda (J.E. Smith) (Lin et al. 2019), which may be due to the fact that insects are poikilotherms and, thus, temperature has a greater effect on their growth and development. Generally, within the proper temperature range, insect growth and fecundity will accelerate with the increase of temperature and reach maximums at the optimum temperatures (Zhang and Lin 2015, Li et al. 2019).

Although our study predicted the habitat suitability areas of two species, there are still some limitations. We randomly sampled background points all around the world. In this case, MaxEnt model will underestimate the potential probability of distribution of the target species at locations that have a suitable environment but no occurrence records were reported. Moreover, it is well known that the distribution of species is determined by multiple abiotic and biotic factors. In this study, only the climatic factor was considered, and other abiotic factors (e.g., topography) and biotic factors (e.g., host plants, competition, and natural enemies) were not considered. This may have resulted in an overestimation of the distribution range of habitat suitability areas for the two species. Thus, to improve the accuracy of model predictions, the various factors mentioned above should be considered.

This research was funded by the Science and Technology Foundation of Guizhou Province (Qiankehe foundation [2017] 1001) and Guizhou Province Educational Science and Technology Research (Guizhou Province KY Character [2016] 279, [2016] 278).