The New Dominator of the World: Modeling the Global Distribution of the Japanese Beetle under Land Use and Climate Change Scenarios

: The spread of invasive species is a threat to global biodiversity. The Japanese beetle is native to Japan, but alien populations of this insect occur in North America, and recently, also in southern Europe. This beetle was recently included on the list of priority species of European concern, as it is a highly invasive agricultural pest. Thus, in this study, we aimed at (i) assessing its current distribution range, and identifying areas of potential invasion, and (ii) predicting its distribution using future climatic and land-use change scenarios for 2050. We collected species occurrences available on the citizen science platform iNaturalist, and we combined species data with climatic and land-use predictors using a Bayesian framework, speciﬁcally the integrated nested Laplace approximation, with a stochastic partial differential equation. We found that the current distribution of the Japanese beetle was mainly, and positively, driven by the percentage of croplands, the annual range of temperature, habitat diversity, percentage of human settlements, and human population density; it was negatively related to the distance to airports, elevation, mean temperature diurnal range, wetlands, and waters. As a result, based on current conditions, the Japanese beetle is likely to occur in 47,970,200 km 2 , while its distribution will range from between 53,418,200 and 59,126,825 km 2 , according to the 2050 climatic and land-use change scenarios. We concluded that the Japanese beetle is a high-risk invasive species, able to ﬁnd suitable conditions for its colonization in several regions around the globe, especially in light of ongoing climatic change. Thus, we strongly recommend strict biosecurity checks and quarantines, as well as regular pest management surveys, in order to reduce its spread.


Introduction
The spread of invasive species is a threat to global biodiversity. Recently, globalization has promoted the ongoing increase in alien species [1], and in addition to being one of the biggest threats to ecosystems [2], biological invasions are also very costly to the global economy [3]. Moreover, the impact of biological invasions on biodiversity and ecosystems has resulted in numerous management and control issues [4][5][6].
Generally, the distribution and spread of exotic species depend on several factors, such as landscape characteristics, climatic conditions, biotic resistance by native taxa, and human aided dispersal, and thus, the identification of areas at high risk of invasion, and the subsequent monitoring of these sites to prevent further incursion, are fundamental [7]. Preventing the establishment of an invasive species, and its consequent spread, is considered an efficient management strategy compared to eradication, containment, and control [8].
For these reasons, species distribution models (SDMs) have been widely and successfully used to predict current and future distribution of invasive species [9][10][11][12]. SDMs need accurate occurrences throughout the whole species range (including both native and

Presences and Observer-Oriented Pseudo-Absences
We considered all occurrences of Pj collected by citizen scientists during the years 2010-2020 around the globe, extracted from the iNaturalist platform [19]. We downloaded this data using the 'get_inat_obs' function in the R package 'rinat' [20]. We only considered species locations (with geographic coordinates) collected between June and September, as this corresponds to the active biological period of our target species [21]. iNaturalist is a citizen science-based, open-access platform aimed at recording biodiversity worldwide, and it allows for the easy downloading of species occurrences using specific queries (i.e., place, date, taxon, observer, etc.).
Similar to [14], to select pseudo-absences we listed all the observers of Pj, and then downloaded, from iNaturalist, all the occurrences of all the species (i.e., including both plants and animals, but excluding those of Pj) collected by those observers. We downloaded this data using the 'get_inat_obs_user' function in the R package 'rinat' [20]. Similar to species' occurrences, we only considered observer-oriented (oo-) pseudo-absences collected during the years 2010-2020, also between June and September, with geographic coordinates. Thus, in order to avoid introducing false pseudo-absences in which the target species had not yet colonized the area, we applied two spatial filters: (1) we limited oo-pseudo-absences to countries in which the Pj was observed (Table S1); and (2) to those occurring inside the minimum convex polygons (shown in Figure S1, derived using the function 'mcp' in the R package 'adehabitatHR') [22] estimated around our target species locations, after spatial clustering was carried out using the 'Mclust' function in the R package 'mclust' [23].

Predictor Variables
For the current period (2010-2020), we considered an initial number of 34 predictors (2 topographic, 10 land cover, 19 bioclimatic, and 3 anthropogenic variables) to represent the habitat characteristics of Pj [7] (Table 1). We derived topographic variables from ASTER GDEM [24] available at a spatial resolution of 15 m), while land cover features were derived from the recent European Space Agency Climate Change Initiative Land Cover layers (ESA-CCI 2019, with a spatial resolution of 300 m) [25]. This database provides yearly spatial data for the period 1992-2019, at a spatial resolution of 300 m. We used data from the period 2010-2019, as it best matched the data of our target species. Relating temporally mismatching species locations and predictors could lead to a biased estimation, due to the annual variation in the climatic and land-use conditions [26]. We grouped the initial 21 classes available in the ESA CCI dataset into 10 land cover types (including percentage of human settlements, see below), and from these derived the Shannon habitat diversity index (Table 1). Moreover, we considered 19 bioclimatic predictors derived from the WorldClim2 dataset [27] for the current period (2010-2018 average), at a 2.5 min resolution (≈5 km). Finally, we included the anthropogenic variables encompassed percentages (see above) of human settlements, human population density, available at the spatial resolution of 1 km [28] (averaging the values of the years 2010-2020), and the distance to airports derived from [29,30]. All these predictors were resampled at a spatial resolution of 5 km.
Thus, we estimated the variance inflation factor (VIF) [31], considering all values of the predictors around the globe (not only at the occurrence of the Japanese beetle) to avoid the negative effect of multicollinearity among predictors on SDMs. In detail, we carried out a stepwise procedure that removed predictors until the highest VIF value was <3 [31].
For the year 2050, we considered as constant the two topographic predictors used for the current period, but estimated future land cover based on a recent dataset developed by Clark labs for ESRI (land cover projection 2050, at a spatial resolution of 300 m [32]. This dataset, derived from the ESA CCI dataset mentioned above, provided spatial information of land cover types (the same considered for the current period) from which we also derived the Shannon habitat diversity index for the year 2050 (Table 1). To the best of our knowledge, this is the first time such a high-resolution land cover dataset, available for a global scale for the period 2050, was used to forecast species distribution. Moreover, we considered the same 19 bioclimatic predictors collected for the current period, but estimated for the year 2050 at 2.5 min resolution (≈5 km). These were derived from future climate change scenarios available for the period 2040-2060 [33]. Similar to [34,35], to reduce single GCM uncertainty [36,37], we also considered four climate change scenarios derived from eight general circulation models (GCMs: BCC-CSM2-MR, CNRM-CM6-1, CNRM-ESM2-1, CanESM5, IPSL-CM6A-LR, MIROC-ES2L, MIROC6, and MRI-ESM2-0), representing four widely used representative concentration pathways (RCP 2.6, RCP 4.5, RCP 7, and RCP 8.5) for 2050. These scenarios derived from the sixth assessment of the Intergovernmental Panel for Climate Change [38]. The selected RCPs represent four possible greenhouse gas emission trajectories, ranging from low (RCP 2.6) to high (RCP 8.5), corresponding to increases in global radiative forcing [34]. Specifically, the increase of the global mean surface temperature in 2050 is likely to be ≈1.67 • C under RCP 2.6, ≈1.99 • C under RCP 4.5, ≈2.18 • C under RCP 7, and ≈2.42 • C under RCP8.5 [39]. The anthropogenic variables encompassed percentages of human settlements (derived from land cover projection 2050) and human population density for the year 2050 were extracted from the SEDAC 2000-2100 1-km grid dataset (see above). Similar to the two topographic predictors, we considered distance to airports as constant, as there are currently no projections for the future distribution of airports for the year 2050. All these predictors were resampled at a spatial resolution of 5 km.
Similar to the current period, for the year 2050 we calculated VIF ( Table 1) for each of the RCP scenarios considered.

Species Distribution Models in INLA
To estimate Pj distribution on a global scale, we combined presence and oo-pseudoabsence locations with predictor variables for the current period, using the recently developed method of integrated nested Laplace approximation (INLA) [40]. Compared to other existing Bayesian methods (i.e., the Monte Carlo Markov chain), INLA has several benefits. For example, (i) INLA provides robust and accurate results [41], (ii) is extremely fast compared to other existing frequentist and Bayesian modelling approaches (especially in light of the growing availability of big datasets) [42], and (iii) is one of the very few modelling frameworks that explicitly accounts for spatial autocorrelation among species locations (i.e., dependency among species locations) [43], by incorporating spatial random effects into binomial models (effective in producing SDM-type spatial predictions) [44].
We specifically developed a binomial model in INLA, considering Pj presence/oopseudo-absence as the response variable, with uncorrelated predictor variables as the fixed effect, and we also took into account spatial dependency among species locations using the stochastic partial differential equation (SPDE) approach [45], based on computations using a Gaussian Markov random field representation of the Gaussian field [42].
To avoid a huge number of oo-pseudo-absences bias in the results, we considered the same number of species occurrences to serve as oo-pseudo-absences (see Results section below). We repeated this procedure 10 times and found consistent results for the further analyses. We specified that we did not select a total of 10,000 pseudo-absences (often considered in SDMs) [46] as the number of Pj occurrences exceeded this number (see Results section below).
We tested the predictive accuracy of our model using 10-fold cross-validations, splitting our initial dataset in two random groups, one to train the model, with ≈90% of the locations, and the second to validate them, with ≈10% of the locations [47]. We specifically considered two validation statistics: (i) the area under the receiver operating characteristic curve (AUC), and (ii) the true skills statistic (TSS). The former ranged between 0 and 1 (worse than a random model and best discriminating model, respectively) while the latter ranged between −1 and 1 (higher values indicate a good predictive accuracy, while 0 indicates random prediction).
Assuming that the current occurrence of Pj adequately reflects the habitat requirements of this pest species, we projected its occurrence around the globe, according to the future environmental and anthropogenic conditions of the year 2050. Thus, we converted the resulting continuous maps of current and future distribution into binary models, considering threshold values estimated by maximizing TSS [48,49]. Values higher or lower than this threshold represented sites where Pj is likely to occur or not occur, respectively.

Definition of the Areas Suitable for the Colonization of the Japanese Beetle
SDMs carried out for very large geographic areas do not always completely capture all the details in species distributions, and could potentially project some false presences beyond the actual species ranges [50]. Thus, we combined Pj occurrences with the resulting non-multi-correlated predictors for the current situation, to identify novel environments in which predictions based on SDMs were uncertain and were not considered for extrapolations. Similar to [12], instead of using the multivariate environmental similarity surface (MESS) [51], we used a modified version of MESS (mMESS) [52]. We estimated using the mMESS model, instead of MESS, as the former does not rely on the use of the most dissimilar variable as an indicator of overall similarity, but rather considers all predictors [52]. In mMESS, values equal to or higher than one indicate that a pixel has at least one predictor variable with values outside the range of the species locations, and thus, it should not be considered for the extrapolations of SDMs [52].
Moreover, since the sum of degree days the beetle needs to complete its development from larva into an adult is a strong limiting factor for the distribution of our target species [21,53], we derived a global map of degree days (DD) from the WorldClim2 monthly dataset. As Pj need a minimum of 711 DD, with a soil temperature greater than 10 • C, to complete its life cycle [21,54] we estimated DD as the sum of all monthly temperatures greater than 10 • C, multiplied by the total number of days using the function 'growingDeg-Days' in the R package 'envirem' [55]. Then, based on [54], we considered only areas with a minimum number of degree days above 711.
Similar to the current period, for the year 2050 we calculated mMESS and DD for each of the RCP scenarios considered. Moreover, similar to [35,56], we identified suitable areas accessible to Pj in the future by multiplying the yearly dispersal distance of 59.463 km/year, derived from 1321.4 m/day [57] × 45 days of life/year [58], for the total number of years (n = 30, from 2020 to 2050). The resulting distances were divided by the number of years larvae would take to develop in the adult stage (711 < DD < 1422 two years, DD > 1422 one year) [21], and then used to define accessible areas around current occurrences.

Results
Using the data available in iNaturalist between June and September in the years 2010-2020, we collected a total of 24,721 occurrences of Pj, uploaded by a total of 14,911 observers who also collected a total of 892,121 non-target species occurrences (Figure 1), which we initially considered as oo-pseudo-absences. Thus, our dataset consisted of 11,204 cells (at 5 × 5 km 2 resolution) in which Pj occurred, and a total of 33,376 cells in which the same observers of our target species collected occurrences of species other than our target species (after applying spatial filters related to countries and minimum convex polygon around Pj occurrences).
Land 2022, 11, x FOR PEER REVIEW Figure 1. Study area, the world. Target species locations in red, total observer-oriented pseudo-absences (i.e., observances from locations other than species, collected by the observers of the target species) in green. Figure 1. Study area, the world. Target species locations in red, total observer-oriented pseudoabsences (i.e., observances from locations other than target species, collected by the observers of the target species) in green.
Among the 16 predictors used in the analysis, we found that in 2050, those linked to natural and semi-natural habitats increase, while those linked to water and to human activity decrease (Table S2). This trend is similar in all continents with two exceptions, Africa and Australia, with increasing urban areas and sparse vegetation, respectively (Table S2).
We found that Pj is significantly and positively related to the percentage of croplands, annual temperature range, habitat diversity, percentage of human settlements, and human population density; it is negatively related to the distance to airports, elevation, mean temperature diurnal range, wetlands, and waters (Table 2).  Figure 2). This value included the 6,097,741 km 2 already occupied by the species, and corresponded to the 12.71% of the total suitable areas available. Among the non-native territories, the United States and Canada (North America) are those most affected by Pj, with around 50% of the suitable areas already occupied.   Considering future conditions, we found that Pj occurrence increases in 2050. Indeed, according to land-use and climatic change scenarios, suitable areas for Pj would account for between 53,418,200 km 2 and 59,126,825 km 2 of the globe (Table 4; Figure 3). Most of the suitable territories are in North America and increase according to the increase in temperature, while in Central and South America there are fewer suitable areas available in the future, but their expansion remains constant while temperature increases (Table 4). Table 4. Future distribution (km 2 ) of the Japanese beetle. Suitable and reachable (from the territory currently occupied) areas were estimated alternatively, excluding and accounting for the dispersal ability of the species. Percentage of reachable suitable areas is also shown.  Figure 4).

Discussion
Our results highlight the current distribution and the future worldwide expansion of the invasive destructive pest Pj as a consequence of both climate and land-use change. Due to its strong adaptability, and its highly polyphagous behavior once outside its range, this species has expanded very rapidly, firstly in North America and secondly in Europe, becoming a more serious pest than in its area of origin [54,59,60]. From our models, it emerges that the species will rapidly expand its distribution range in the future, especially because of the general increase in annual temperatures, which both favor the species reproduction, and increase the surface of suitable areas [61]. The future changes in land use seemed not to have a crucial effect in the spread of Pj.

Current Potential Distribution and Effect of Environmental Variables
Based on our results, most of the areas suitable for Pj are not still occupied by the species, especially those located in the Southern Hemisphere. In these areas, Pj has not yet spread [62], although many suitable areas are available. This is most likely due to the asynchrony between the Northern Hemisphere adult activity season (June and July), and the optimal weather conditions for the survival of the adults in the Southern Hemisphere, which is December and January. Therefore, adults that were accidentally transported from the Northern Hemisphere during summer would have found winter conditions in the Southern Hemisphere, with cold temperatures and very few available and edible plants.
North America is the area of the Northern Hemisphere most affected by the species [59]. In particular, in the United States this invasive beetle was accidentally introduced in 1916, and despite management actions to limit its colonization, within the last 100 years it had successfully spread across much of the east of the country [17]. Moreover, according to our results, Pj is extending its range northward into Canada [63], occupying about half of the whole suitable territories of North America. In Europe, Pj was introduced in 2014 and, according to our results, it currently occupies less than 1% of the suitable areas in this continent. However, the aggressiveness of this species, and its rapid spread [57], alarmed European governments, which immediately activated actions to contain the pest and prevent its spread. Currently, Pj is a quarantine pest, designated as a high priority candidate in the new phytosanitary legislation of the European Union [18,64], and is listed in Annex I Part A/1 of Council Directive 2000/29/EC4.
Our results suggest that the areas occupied by the species are those most impacted by anthropogenic activities. In general, human disturbance influences many biological invasions, especially during the earlier stages after introduction [65]. Territories with a high anthropogenic footprint allow the survival of invasive pests in areas with suboptimal (even unsuitable) climate, e.g., via propagule pressure or by the creation of microclimates [66,67]. Indeed, as expected, we found a strong relationship between suitable areas for Pj and their closeness to airports, confirming the crucial role of human-assisted movement in the rapid expansion of this invasive species [68][69][70][71]. In this regard, in Europe the first observations of Pj referred to a small area located in the Ticino Valley Natural Park in Italy, extremely close to Malpensa airport [72].
Human actions, through the modification of landscapes and the conversion of natural areas into croplands, favor the spread of pests [73]. Currently, agriculture is a dominant form of land management globally [74], and therefore, it is not surprising that pest species are rapidly expanding. Pj is not an exception. Based on our results, Pj tends to occupy those territories with a large extension of croplands, but also with a high habitat diversity. This beetle is a generalist pest, thriving in large areas of turf and pasture grass for grubs developing [61], feeding with over 300 host plant species in 79 families . Therefore, Pj most likely finds those areas with a greater variety of habitat and plant species, including those of non-economic importance, more suitable to optimize the survival of both grubs and adults. The ability to feed on several hosts also ensures the survival of the species, and explains its ability to adapt to new environments [76,77].
According to already published literature [17,78], among the climatic factors, temperature is the parameter that most influences the expansion of this invasive beetle. Moreover, specific temperature is required for oviposition and larval development [15,53,79]. Thus, our models confirm the role of temperature in affecting the spread of Pj, and show that areas with high annual temperatures, and low daily temperature fluctuations, are highly suitable for the beetles. Consequently, we can also include areas with temperate (maritime or sub-continental), Mediterranean, and tropical climates in the areas that are suitable for Pj.
Although the role of soil moisture has been recognized as a key parameter to limit the potential spread of Pj [15,80], in our models precipitations have little effect on its distribution. Therefore, based on our climatic projections, we identified suitable areas for the establishment of Pj that also included areas previously considered unsuitable because of the lack of summer rainfall, such as the Mediterranean regions [80]. Soil moisture depends on many factors besides rainfall, such as soil properties [81,82], proximity to water table [83], and the intensity of water consumption by plants [84]. All these aspects are strictly linked to crops productivity [85], and strongly influence the development of Pj grubs [18,86,87]. Therefore, fertile areas, such as the Po Valley (Italy), are certainly more suitable than others for the proliferation of the beetle, regardless of the intensity of atmospheric precipitations.

Future Potential Distribution
One of the peculiarities of Pj is the ability to undergo a progressive acclimatization to a broad range of environmental conditions and human pressure [7]. In 2050, the new environmental conditions, combined with a rise in temperature between 1.6 • and 2.4 • C, will lead to an increase of territories suitable for Pj of up to more than 20% of the current territories. North America is the area with most of the new suitable regions for Pj (up to 59% of the current areas), followed by Asia and Europe. Based on our results, predictor variables such as croplands and human density, positively related to Pj occurrence, decrease in the future, without producing the expected decrease in areas suitable for this species. This is most likely due to the contrasting effects of climate over anthropogenic factors, as already reported in the literature [7]. Indeed, the rise in winter temperature, reducing cold stress especially in areas of the Northern Hemisphere such as Canada and Russia, would change previously unsuitable areas into new climatically suitable ones [17]. Thus, future conditions enable the beetle to complete its life cycle in a single year instead of two [88,89]. This last aspect leads to a considerable increase in the beetle dispersal capacity, with an increase in its ability to colonize new suitable areas in a short time.
Conversely, increasingly hot and dry future conditions reduce the previously suitable range of areas (Petty et al., 2015 [17,90]. However, in contrast to previous studies [17], the narrowing of the suitable range for Pj does not involve all continents in the Southern Hemisphere, but only Central and South America experience a decrease of about 10% of their current suitable territories. In the rest of the Southern Hemisphere, suitable areas increase but only slightly in Australia. In 2050, even considering the high dispersal capacities of Pj, less than half of its future available areas are reachable. However, our predictions do not take into account future accidental introductions of Pj into new territories (this is not statistically predictable) that could lead to a significant increase in distribution.
Considering that the surface currently occupied by Pj is only 12.71% of the total available, any expansion in the next 30 years is alarming. In North America, Pj would colonize almost all suitable areas, while in Europe it would reach about half of its potential range. In Asia, in the next 30 years, Pj would reach only 10% of the suitable areas. Indeed, in this continent, despite having a suitable surface greater than that of Europe, and almost equivalent to that of North America, large future suitable areas are not accessible to Pj, currently located in the islands (Japanese and Russian), or in sub-optimal areas (Southern India). However, Pj would be able to successfully establish itself in the mainland of Asia, by way of accidental introduction, due to the new climatically favorable habitats outside the distribution of their natural enemies [91].

Conclusions
In this study, we identified with great precision the areas that the pest species could reach in the next 30 years. Therefore, we provided useful information to direct all the resources and effort necessary to prevent the introduction and spread of Pj.
However, the distribution and spread of exotic species also depends on a number of factors not considered in our models, such as biotic resistance by native taxa, human aided dispersal, and to a lesser extent, the availability of host plants that limit Pj expansion [7,49], and which have not been well studied or understood. Therefore, we encourage further research, especially aimed at deepening the ecology of the species, and we advocate the further development of new statistical models in which additional biotic factors can be incorporated. Finally, we strongly suggest adopting approaches similar to those used in this study, possibly combining species occurrences collected during standardized sampling on a local scale when available (e.g., in INLA SPDE, including the source of data as a random effect), in order to model invasive species distribution through space and time, and provide ecologically and statistically robust estimates of species occurrence.