Forecasting the spread associated with climate change in Eastern Europe of the invasive Asiatic flea beetle, Luperomorpha xanthodera (Coleoptera: Chrysomelidae)

The current climate has facilitated biological invasions worldwide, with the spread of invasive species accelerating over the last few decades. Introductions of species are often accidental and facilitated by many factors, including intercontinental trade. In this paper, we explore the secondary range of an adventive fl ea beetle native to Asia, Luperomorpha xanthodera, which was recently found in many European countries on several occurrences. This species has detrimental effects mainly on ornamental plants, being destructive to leaves and fl owers, which can have serious negative economic consequences. We inferred current and future potential distribution of this alien species using Ecological Niche Modelling (ENM) and analysed the future increase in suitable areas in European countries. Suitable climatic conditions for L. xanthodera are assumed to be similar to that in areas where this species currently occurs, which includes Ireland and some Balkan countries, where the species is not recorded yet. In the future, a North-eastern expansion is predicted, with many countries currently lacking suitable climatic conditions for L. xanthodera becoming suitable and potentially colonizable by this species. The geographic trend revealed and quantifi ed in our analysis follows the increase in the mean temperature in the coldest quarter of the year, which is predicted to rise in the next 30 years. This will result in this species spreading from south-western Europe to North-Eastern countries, such as Poland, Ukraine, Romania, Belarus and Latvia, which should adopt preventive measures to avoid the accidental introduction of L. xanthodera.


INTRODUCTION
The ongoing climate change is one of the most important phenomena shaping the current distribution of species, of both conservation and invasive concern (Bellard et al., 2012;Cerasoli et al., 2019;Wieczynski et al., 2019). The spread of invasive species of all taxa has increased during the last few decades (Seebens et al., 2017). Although preventive measures were adopted in both high-and low-income countries, such controls on trade in pets, stowaways in international trade and marine ports in most of them are not able to limit the increase in the number of possible unwanted introductions (Early et al., 2016). Alien species can negatively affect both biodiversity and agriculture, which generate expensive management issues (Allison et al., 2009;Paini et al., 2016;Gilioli et al., 2017).
Luperomorpha xanthodera is of Asian origin and an adventive species in Europe since 2003(Johnson & Booth, 2004Bieńkowski & Orlova-Bienkowskaja, 2018b). Its area of origin (primary range) includes China and the Korean Peninsula [even though this last is doubtful according to Löbl & Smetana (2010)], while the secondary range currently includes Great Britain, Spain, France, Italy, Germany, The Netherlands, Hungary, Austria, European Russia (Bieńkowski & Orlova-Bienkowskaja, 2018a, b) and Montenegro, even though for this country the establishment of this species is doubtful (Radonjić & Hrnčić, 2017) (Fig. 1A). This fl ea beetle has the potential for becoming locally abundant and destroying ornamental plants, especially the fl owers (Kozłowski & Legutowska, 2014). In Central Italy (Tuscany) it produces two, partially overlapping generations and overwinters in all post-embryonic stages (Del Bene & Conti, 2009). Adults appear in spring and are present on plants throughout the whole vegetative season. Larval development takes place in the soil and the damage they do the roots does not seem to affect infested plants (Kozłowski & Legutowska, 2014). More distinct damage is caused by adults that feed on plants, especially the fl owers and leaves of seedlings (Del Bene & Conti, 2009).
In this paper, we infer potential future distribution of L. xanthodera in terms of gains and losses associated with the predicted changes in the geographical trends in the emissions of different greenhouse gases (GHGs). Further, the results are mainly for the countries within the species' secondary range and will help in the development of national management strategies where this species currently occurs and is predicted to occur in the future.

Species studied
Chrysomelidae is one of the largest Coleopteran families, including over 37,000 species worldwide (Jolivet & Verma, 2002).

Fig. 1.
A -localities where Luperomorpha xanthodera occurred in its primary (blue) and secondary (red) ranges; the question mark indicates the doubtful occurrence in Montenegro. B -predicted habitat suitability for L. xanthodera under current climatic conditions ranging from low (0 -blue) to high (1 -red).

Dataset and study area
The localities where L. xanthodera occurs were obtained from the two most recent papers on this species in both its primary (Wang et al., 2010) and secondary (Bieńkowski & Orlova-Bienkowskaja, 2018a) ranges, along with some new data from Italy (Lombardia and Veneto, M. Biondi, unpubl.). These papers were chosen because they are the most recent and comprehensive studies on this species. The dataset includes 39 localities for the primary and 46 for the secondary ranges, respectively, of this species (Fig. 1A).
The whole dataset is available in Table S1, including coordinates (in WGS84) and indication of their occurrence in either the primary or secondary ranges.

Ecological niche modelling
Recently Ecological Niche Models have become more popular among ecologists for predicting the potential distribution in space and time of a species, and identifying the environmental variables associated with the distribution of taxa, even using the presence-only data reported in bibliographic sources, online databases and museum records (e.g. Ferrer-Sánchez & Rodríguez-Estrella, 2016;D'Alessandro et al., 2018;Brunetti et al., 2019). Included in the Ecological Niche Models is a set of nineteen bioclimatic variables, which were downloaded from the online repository Worldclim.org, at a resolution of 30 arc-seconds, for the 'current' climatic conditions (ver. 1.4) (Hijmans et al., 2005); names and codes for each variable are in Table S2. For predicting this species' future distributions, three different sets of bioclimatic variables for future climatic conditions (2050 and 2070) were also downloaded, in order to address differences among Global Climate Models (GCMs, models represent the interac-  (Brunetti et al., 2019;Iannella et al., 2019a). For both 2050 and 2070, the Representative Concentration Pathways (RCPs, represent the trajectories climate will take due to radiative forcing) 6.0 and 8.5 were used to model this species' response to different scenarios of greenhouse gas (GHGs) emissions. The RCP_6.0 is the result of the increased radiative forcing due to GHGs emissions decreasing only after 2080, while RCP_8.5 considers the emissions continuing to at least 2100 (Meinshausen et al., 2011;Riahi et al., 2011), thus, considering both a "plausible" and an "extreme" climatic scenario.
To avoid multicollinearity among predictors, a correlation matrix was built from which pairs of variables exceeding Pearson's |r| > 0.85 were removed (Elith et al., 2006). In addition, the occurrences database (both for primary and secondary ranges) was also checked for correlation among records using the Moran test in ArcMap 10.0 (ESRI, 2010).
Models were built using the Maxent algorithm implemented in the SDMtoolbox ver. 2.4 (Brown et al., 2017) in ArcMap 10.0; this tool permits the use of the high performance Maxent algorithm (Elith et al., 2006;Wisz et al., 2008) coupled with spatial jack-knifi ng, a procedure which geographically performs the cross-validation (Brown, 2014) used to model several conservation and biogeographical issues (Kumar et al., 2015;Iannella et al., 2018b;Hill & Winder, 2019). A sampling bias mask was generated using the 'Gaussian kernel density of sampling localities' (sampling bias distance = 50 km) in order to avoid the more sam- pled areas having a negative effect on the models (Phillips et al., 2009). Further, all predictors were projected into an Equal Area Projection (WGS84 -Cylindrical Equal Area), to correct for any latitudinal bias (Brown, 2014); spatial jack-knife was set at 5 replicates per parameter and 3 spatial groups were created. Models were calibrated for both primary and secondary range localities of this species, a technique used to improve performance of predictions and reduce uncertainties (Broennimann & Guisan, 2008) when modelling the distribution of invasive species (Cerasoli et al., 2019;Iannella et al., 2019a;Perret et al., 2019;Rodríguez-Rey et al., 2019).
The occurrences of spatial autocorrelations were tested for using the Moran test in ArcMap 10.0 ['Spatial Autocorrelation (Global Moran's I)' toolbox]; this test identifi es patterns of regular occurrence, clustering or spatial independence of the localities where this species is present (ESRI, 2010).
This workfl ow was used to both lower geographic and possible self-correlation biases in the occurrence data and facilitate tests of the different types of features and regularization of multipliers, resulting in higher model performances (Brown et al., 2017).
For each model based on the different GCMs used, the Area Under the Curve (AUC) of the Receiver Characteristics Operator (ROC) was used as a weight in the Multivariate Environmental Dissimilarity Index (MEDI) algorithm , which proportionally combines different GCMs projections, starting with the MESS (Elith et al., 2010) maps produced by Maxent, which down-weights the extrapolations of the models. Both current and future predicted scenarios were binarized using the 10 th percentile of training presences, which is a reliable and conservative threshold when data are collected at different times by several observers (Freeman & Moisen, 2008;Rebelo & Jones, 2010;Urbani et al., 2015Urbani et al., , 2017.

Post-modelling analyses
To determine possible changes in range, current and future binarized maps were compared using the 'BIOMOD_RangeSize' function in the 'biomod2' package (Thuiller et al., 2016) in R (R Core Team, 2016); this function calculates loss, similarity and gain of area when comparing two binarized maps, for instance a current and a future map, as in our case. Gains were analysed using ArcMap 10.0 on a latitude-longitude grid (1 square degree) to measure possible shifts and areas in each cell. All geographic data and the corresponding maps were managed in ArcMap 10.0.

RESULTS
The Moran's test of the localities of the occurrences resulted in I = -0.0262, with z-score = 0.000880 and p = 0.9992 for the primary range and I = -0.0182, with z-score = 0.1112 and p = 0.9114 for the secondary range, which indicates random distributions in both ranges. The correlation matrix of all the candidate predictors resulted in eight variables being selected for modelling, namely BIO1  (Table  S2). Models had high AUC values, with an average value = 0.915 and standard deviation = 0.014, which indicates that habitats mainly in Central Europe, British Isles, parts of Italy and the Balkans, and a small are in Southern Russia (Fig. 1B) are currently suitable. The most important variables are: Mean Temperature in Coldest Quarter (BIO11 = 63.4%) and Precipitation in Driest Month (BIO14 = 22.2%), with the fi rst showing a preference for mean temperatures between 5.5° and 11.1°C during the coldest months and the second a preference for precipitation of > 50 mm (Fig. S3).
The 10th percentile training presence calculated for each replicate resulted in an average value of 0.384, which was used as the threshold for the binarization of current and future predictions. The future predictions (Fig. 2) show an increase in the area suitable for L. xanthodera in Central Europe and a noticeable trend in a greater suitability in Eastern European Countries and European Russia, apparently resulting from an increase in the most important variable BIO11 (Fig. 3). The trend is obvious in the range-shift maps (Fig. 4A-D), where areas gained are mainly in the aforementioned areas, especially in RCP_8.5 (Figs 4B, D). The latitude and longitude analysis using these maps further confi rmed this trend (Fig. 4E, F) and shows that a parallel increase in areas gained is predicted in terms of both latitude and longitude. Considering this species' secondary range, the results also indicate a north-eastern shift in terms of countries hosting current and future suitable conditions (Fig. 4). The countries currently with highly suitable areas are mainly France, Germany, the British Isles and Italy (Fig. 5), which host 69.3% of the total binarized suitable area. The future predictions indicate that the country's most likely to be colonized by this species are the East-European ones, mainly Poland, Ukraine and Romania (Figs 5A, B), with noticeable peaks in the RCP_8.5 of Latvia and Belarus (Fig. 5B). In terms of the net gain for each country (net gain = gain -loss) and corresponding total area of the secondary range colonized, Poland, Ukraine and Romania reach 32.5% in 2050 and the 41.5% in 2070 based on RCP_6.0. By adding the RCP_8.5 peaks of Latvia and Belarus, the sum is 53.9% in 2050 and 51.5% in 2070. Negative net gains are reported only for Spain and Turkey for both RCP_6.0 and RCP_8.5 scenarios, with the exception of a slight gain in 2070 based on RCP_8.5 (Fig. 5B).

DISCUSSION
The ENMs for current climatic conditions in the secondary range of this species reveal that currently suitable territories do not extend as far as the occurrences currently recorded, except for Ireland and some Balkan countries (Fig. 1B). On the contrary, Poland does not seem to have any suitable areas for L. xanthodera (Fig. 1B), even though Kozłowski & Legutowska (2014) state that this species could become locally abundant. Interestingly, occurrences in the British Isles do not fully match the climatic potential of L. xanthodera, a condition recorded also for the Balkans; therefore, an accidental introduction of this species into these areas could be particularly dangerous as the climate there is suitable for the establishment of this species. Indeed, more attention should be paid to the northern Atlantic coasts of Spain and Portugal, and Italy and Balkan countries, to identify possible introductions of this species as our results indicate that it is highly likely that L. xanthodera is already in these countries. Such unrecorded occurrences limit the present study, as these could confi rm (or change) the ENMs' results. This issue is common to invasive species of leaf-beetles, whose distributions are scattered at particular localities and are not recorded every year (Bieńkowski & Orlova-Bienkowskaja, 2018a).
The mean temperature in the coldest quarter (BIO11) seems to be important in determining a species' potential distribution, based on both its high contribution in the ENMs and its similar importance in the future (Fig.  3). This result is in line with other published trends, in which low temperatures are important in determining the distribution of insects with some invasive alien insects favoured by an increase in the low temperatures that are likely to occur with global warming in the future (Orlova-Bienkowskaja & Bienkowski, 2016;Iannella et al., 2019b). In marked contrast, high altitude and endemic species are not favoured by the predicted temperature increase (Urbani et al., 2017;Brunetti et al., 2019;Cerasoli et al., 2020). Depending on the climatic scenario and RCP used, an increase in BIO11 is observable from south-western to north-eastern Europe, which is very similar and concurrent with the future predictions (Fig. 4). L. xanthodera, which currently produces two partially overlapping generations in Italy (Del Bene & Conti, 2009), may thus become more invasive in the future at high latitudes in Italy, where it currently overwinters. Considering the trophic preferences of this species and the plants it is found on (Del Bene & Conti, 2009;Vincent & Doguet, 2011), the hypothesis of an invasion of Eastern Europe's market gardens (Beenen & Roques, 2010) is corroborated by our future projections, in which the climatic suitability clearly increases throughout this species' secondary range. The geographic trend in area that will become suitable in the future shows an impressive increase in Eastern Europe, with countries currently not suitable becoming highly suitable like France, Germany and U.K. Poland, Ukraine and Romania, as well as Belarus and Latvia should adopt special measures to prevent the accidental introduction of L. xanthodera. In fact, on average, all the future predictions (2050 and 2070) and different climatic scenario (6.0 and 8.5 RCPs) forecast a high and diffuse increase in the area suitable for this species (Figs 2, 4).
Finally, this trend is similar to other invasions (Grapputo et al., 2005;Orlova-Bienkowskaja, 2013Ukrainsky & Orlova-Bienkowskaja, 2014), which leads us to conclude that special attention should be given to controlling the pests on ornamental plants, especially those being traded with eastern Europe.