Predicting the global distribution and invasion scenarios of the Spotted Lanternfly, Lycorma delicatula (White, 1845) (Hemiptera, Fulgoridae)

267 Predicting the global distribution and invasion scenarios of the Spotted Lanternfly, Lycorma delicatula (White, 1845) (Hemiptera, Fulgoridae) Enrico Ruzzier1,2 , Davide Scaccini3, Pietro Tirozzi4, Valerio Orioli4, Olivia Dondina2,4 , Andrea Di Giulio1,2 , Alberto Pozzebon3, Luciano Bani2,3 1 Department of Science, Roma Tre University, Rome, Italy 2 National Biodiversity Future Center (NBFC), Palermo, Italy 3 Department of Agronomy, Food, Natural resources, Animals, Environment (DAFNAE), University of Padua, Padova, Italy 4 Department of Earth and Environmental Sciences, University of Milano-Bicocca, Milan, Italy Corresponding authors: Enrico Ruzzier ([email protected]); Olivia Dondina ([email protected]) Copyright: © Enrico Ruzzier et al. This is an open access article distributed under terms of the Creative Commons Attribution License (Attribution 4.0 International – CC BY 4.0). Research Article Abstract Lycorma delicatula is a species native to the PR of China that has become invasive in South Korea, Japan, and the United States. It is considered a significant threat to agriculture, particularly the viticulture industry, and its further spread into new areas could exacerbate its economic impact. Species Distribution Models (SDMs) are commonly used to analyze the distribution of invasive species, and while L. delicatula has already been studied using this spatial modeling approach, previous research has mostly focused on the effects of bioclimatic variables, often overlooking the role of habitat characteristics, including the presence of host plants. In our study, to assess the presence of other suitable habitats for the species on a global scale, we developed a two-step SDM calibrated within its native range, within biogeographic barriers. In the first step, we built two separate models: one for habitat suitability (HSM), incorporating land cover, elevation, and host plants, and another for bioclimatic suitability (BSM). The HSM calibration included the allocation of a portion of pseudo-absences (PAs) along background points (BPs) to mitigate the sampling bias of occurrences concentrated in highly urbanized areas (one bias-controlled PA for every four random BPs). In contrast, the BSM calibration was performed using a fully random allocation of BPs. In the second step, the two models were combined to produce an overall suitability map for the native range, which demonstrated excellent validation performance. When projected globally, the model confirmed that the species has already colonized suitable areas. Beyond the currently invaded regions (i.e., North America and East Asia), the model identifies additional colonizable areas exclusively in Europe. To evaluate the invasion dynamics in both currently invaded regions and potential future invasion areas in Europe, we developed a resistance-distance constrained dispersal model based on a) environmental resistance, derived from the overall suitability map produced by the SDM, and b) intrinsic dispersal distance, estimated from the temporal progression of occurrences in the primary invasion area in the eastern United States. The maximum annual dispersal distance was estimated to be approximately 25 kilometers. Given its relatively low annual dispersal capacity, the species appears capable of naturally spreading within areas of medium to high suitability, while its potential for expansion into low-suitability areas remains quite limited. Mediumand high-suitability regions where the species is already present are likely to undergo significant colonization over the next 10 years. In contrast, adventive populations in low-suitability areas seem unable to expand successfully without human-mediated translocations. In the absence of an adaptive capacity to thrive in ecological contexts different from its native range (niche change or niche shift), the species appears to be a relatively modest invader beyond its currently invaded Academic editor: Samuel Ward Received: 29 March 2025 Accepted: 11 September 2025 Published: 4 November 2025 Citation: Ruzzier E, Scaccini D, Tirozzi P, Orioli V, Dondina O, Di Giulio A, Pozzebon A, Bani L (2025) Predicting the global distribution and invasion scenarios of the Spotted Lanternfly, Lycorma delicatula (White, 1845) (Hemiptera, Fulgoridae). NeoBiota 103: 267–298. https://doi. org/10.3897/neobiota.103.154246 NeoBiota 103: 267–298 (2025) DOI: 10.3897/neobiota.103.154246 Advancing research on alien species and biological invasions A peer-reviewed open-access journal NeoBiota 268 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly areas or Europe. Successfully managing the invasion of L. delicatula therefore depends on the implementation of early detection and eradication measures, which should be deployed as promptly as possible, especially in highly suitable introduction areas. Key words: Auchenorrhyncha, bioclimatic suitability model, habitat suitability model, host plants, insect pest, invasive species, Resistance-Distance Constrained Dispersal Model, Species Distribution Model (SDM) Introduction The Spotted Lanternfly, Lycorma delicatula (White, 1845) (Hemiptera: Fulgoridae) is a planthopper endemic to the PR of China that has become widely invasive in Japan, South Korea, and the United States (Lee et al. 2019; Urban and Leach 2023). The species is of primary phytosanitary concern given its impact on crop production caused by feeding nymphs and adults (Lee et al. 2019; Urban and Leach 2023; Zhang et al. 2023). Indeed, both nymphs and adults of L. delicatula feed on plant sap, affecting the health of the plant (Dara et al. 2015). Furthermore, the gregarious behavior of this species, with aggregations of hundreds of individuals, further amplifies the damage to the host plant (Lee et al. 2019). Lycorma delicatula is widely polyphagous, feeding on multiple plant species, both in its native and invaded range (CABI 2020; Jung et al. 2022; EPPO 2023), although Ailanthus altissima (Mill.) Swingle (commonly known as Tree of Heaven) is often referred to as the primary or preferred host (Nixon et al. 2022; Elsensohn et al. 2023). As a consequence, this insect currently occurs in various anthropogenic environments, such as agricultural and urban areas, as well as in natural and semi-natural habitats (Urban and Leach 2023). Although the Spotted Lanternfly is potentially harmful to many plant species, multiple studies carried out in the invaded areas have shown that L. delicatula is only moderately harmful to hardwood and fruit trees (Hoover et al. 2023; Nixon et al. 2023a) while it is a primary pest on grapevine plants (Vitis spp.) (Harner et al. 2022). In the eastern United States, L. delicatula prefers vineyard-dominated landscapes (Harner et al. 2022) where it significantly contributes to yield losses and grapevine decline (Urban 2020). These conditions, therefore, make this pest a real threat to the wine and grape industry on a global scale, both in countries where it has already been introduced and in those where it could be in the near future (e.g., Reynolds et al. 2021; Jones et al. 2022; EPPO 2023). Management of this pest predominantly relies on increased use of broad-spectrum insecticides with high residual toxicity, posing significant threats to beneficial non-target organisms and the environment (Lee et al. 2019; Urban and Leach 2023). However, this type of management may be ineffective over the long term, as treated areas can be recolonized by resident populations in adjacent non-crop habitats (Park et al. 2009). It is crucial to determine the potential distribution of L. delicatula on a global scale to identify the areas with the highest establishment risk and to guide more effective early detection strategies and pest management actions (Urban and Leach 2023). In this regard, species distribution models (SDMs) have become effective tools for forecasting the potential distribution of invasive insects (e.g., Elith 2017). To date, several studies predicted the distribution of L. delicatula through SDMs, both locally and globally, albeit through different approaches and datasets (Jung et 269 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly al. 2017; Namgung et al. 2020; Wakie et al. 2020; Huron et al. 2022; Chartois et al. 2024; Zhao et al. 2024). In addition, Jones et al. (2022) and Montgomery et al. (2023), beyond assessing species suitability, also estimated the risk of invasion by accounting for connectivity under current and future scenarios. These studies, in some cases, estimated suitability based on bioclimatic variables, while in others, they employed a modelling framework that tended to overemphasize the influence of these variables over habitat features (e.g., elevation, host plants). As a result, the interpretation of L. delicatula potential distribution is often confined to a broad scenopoetic context rather than a local scale. Recent studies highlight that incorporating habitat variables into a two-step framework, separating a habitat suitability model (HSM) from a bioclimatic suitability model (BSM), enhances SDM forecasting (Adde et al. 2023; Ruzzier et al. 2024). For phytophagous species, integrating host plants into habitat modelling improves the reliability of predictions for potentially suitable areas, both within and outside the native range of the target species (Dang et al. 2021; Ruzzier et al. 2024, 2025). Studies analyzing the movements of adult females in search of oviposition sites (Wolfin et al. 2019) or food resources (Wolfin et al. 2020) highlight the low propensity of L. delicatula to cover long distances. However, the dispersal potential of propagules remains unknown. Gaining insight into the spatial expansion dynamics of propagules could support the development of effective containment strategies, particularly in relation to already infested areas and the suitability of adjacent habitats. The aims of this research were: 1) to separately calibrate a BSM and a HSM within the native range of L. delicatula as a basis to produce an overall environmental suitability map; 2) to forecast suitable areas for the species at global scale based on the models (BSM and HSM) calibrated to the native range, with a particular focus on invaded and newly colonized areas; 3) to evaluate the model effectiveness in predicting the current distribution of this species in invaded areas (USA and South Korea); 4) to estimate the annual dispersal distance of L. delicatula based on repository occurrences; 5) to model the invasion range expansion and forecast the future USA, South Korea and Japan invasion scenarios, as well as the potential invasion scenario in Europe. Materials and methods Species occurrences Occurrences of L. delicatula in its native and invaded ranges were retrieved from the Global Biodiversity Information Facility (GBIF) (https://doi.org/10.15468/ dl.azzdjh; 32,752 occurrences) and iNaturalist (41,867 occurrences), accounting for 74,619 occurrences worldwide. The occurrence data used in the SDMs were updated to October 11, 2024. The dataset was checked for all potential issues (e.g., duplicate records, misidentifications, etc.). For the BSM calibration on native range (mainland PR of China), we used all available occurrences, while for the HSM, we used a subset of the previous data, selecting only those with coordinates in degrees with at least two decimal places of resolution or accuracy lower than 270 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly 2.5 km (half of the raster cell size of the environmental layers used in the HSM model) (Ruzzier et al. 2025). The rational of the procedure relies on the capability of associating each occurrence with habitat data derived from land cover cartography with a grid resolution of 0.04167° (corresponding to approximately 5 km at the equator). For BSM, occurrences with higher uncertainty were also included, as bioclimatic variables tend to exhibit greater spatial inertia in their variability. The overall dataset was thinned through the function “ensemble.spatialThin” of the package “BiodiversityR” (Kindt and Coe 2005), with a thinning parameter of 2.5 km in R, version 4.4.2 (R Core Team 2023) to exclude duplicated data and remove occurrences within the same raster cell. For HSM and BSM 324 and 402 occurrences were used, respectively. Occurrences from Japan, the Korean peninsula, and Taiwan were excluded from the analyses since these areas are outside the putative native range of this species; in these countries, L. delicatula is naturally absent, because of the existence of major biogeographic barriers (ocean and mountain ranges) that the species is unable to overcome, whereas present occurrences are the result of human introduction. Bioclimatic, elevation and habitat variables Bioclimatic (19 strata) and elevation (i.e., digital elevation model, DEM) layers were downloaded from WorldClim 2.1 (available at: http://www.worldclim.com/ version2) and averaged for the years 1970–2000, with a grid resolution of 0.00833°. Since most Bioclimatic strata are strongly collinear, a preselection of bioclimatic variables may help in reducing multicollinearity, improving model interpretability, and enhancing predictive accuracy (Dormann et al. 2013; Pliscoff et al. 2014). We preselected a subset of 10 among them: mean diurnal range (mean of monthly (max temperature – min temperature)) [BIO2]; max temperature of warmest month [BIO5]; min temperature of coldest month [BIO6]; mean temperature of warmest quarter [BIO10]; mean temperature of coldest quarter [BIO11]; annual precipitation [BIO12]; precipitation of wettest quarter [BIO16]; precipitation of driest quarter [BIO17]; precipitation of warmest quarter [BIO18]; precipitation of coldest quarter [BIO19]. Habitat variables were retrieved from The Global 1-km Consensus Land Cover available at http://www.earthenv.org. Data were represented by 12 layers (evergreen/deciduous needleleaf trees [EDNT], evergreen broadleaf trees [EBT], deciduous broadleaf trees [DBT], mixed/other trees [MOT], shrubs [SHR], herbaceous vegetation [HERB], cultivated and managed vegetation [CMV], regularly flooded vegetation [RFV], urban/built-up [BU], barren lands [BR], open waters [OW] and Snow/Ice [SI]), each of which provides the prevalence of one land-cover class, expressed as a percentage in a pixel having a resolution of 30 arc-seconds (i.e., 0.00833° resolution). To reduce computational effort, bioclimatic, DEM, and land cover layers were upscaled to a 0.04167° resolution via the “aggregate” function included in the terra package (Hijmans 2024). Among the 12 land-use layers, we excluded a priori SI and OW from the HSM modelling, as they represent unsuitable environments for the target species and for insects in general. To habitat variables, we added a layer of host plants to test whether this covariate could improve the predictive effectiveness of HSM (see Ruzzier et al. 2024, 2025). To this aim, we compiled the host plants list of L. delicatula (see Suppl. 271 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly material 1: table S1) using information obtained from the European and Mediterranean Plant Protection Organization (EPPO 2015), whose dataset is based on scientific literature, expert assessments, and official reports provided by National Plant Protection Organizations of its member countries. If more than one species belonging to the same genus was known as a host plant, or if congeneric species could be used as vicariant host plant species in areas at risk of invasion or already invaded, the genus was used as the host taxon (e.g., Acer). Otherwise, a single species was used as the host plant (e.g., A. altissima). For host plants identified at the genus level, all records of each species within the genus were downloaded from GBIF using the “occ_data” function in the rgbif R package (Chamberlain et al. 2024). This same approach was applied for host plants at the species level. Data for each host plant taxon were combined into a single layer of georeferenced points using the “as.ppp” function to create a density map via kernel interpolation with the “density.ppp” function; the kernel interpolation radius was estimated with the “bw.ppl”; all functions pertaining to spatstat R package (Baddeley and Turner 2005; Baddeley et al. 2005). To reduce potential oversampling, we used the logarithmic transformation of the kernel density maps (representing the quantitative range map of each host plant taxon). To prevent overestimation of each host taxon distribution, we multiplied the logarithmic density layer by the land cover type associated with each taxon habit (e.g., deciduous broadleaf trees) or by the land cover type relevant for economically valuable taxa (e.g., cultivated and managed vegetation); the association between host plant taxa and land cover is indicated in Suppl. material 1: table S1. Finally, we produced an overall host plants layer by extracting the maximum density value found in each grid cell across all host plant layers. Since host plants can locally replace each other, we used the maximum density value rather than summing the density values of all layers to avoid overestimating host plant density. Data on host plant taxa were downloaded from GBIF on September 30, 2024. It should be acknowledged that host plant occurrences are shaped by bioclimatic conditions, and that the kernel density map implicitly incorporates these effects. As a result, bioclimatic influences on the insect species may be indirectly accounted for in the HSM. To avoid collinearity issues, all environmental covariates used to build the HSM and BSM were subjected to variance inflation factor (VIF) analysis using the “vifcor” function from the usdm R package (Naimi et al. 2014). Highly correlated variables, those with a Pearson’s correlation coefficient |r| > 0.7, were excluded (Dormann et al. 2013). This analysis was conducted separately within the calibration extents of BSM (broad-range extent: 85°E to 125°E; 0°N to 70°N) and HSM (narrow-range or native extent: 90°E to 125°E; 25°N to 45°N). None of the land use covariates used in the SDMs were pairwise correlated. Following the VIF analysis, six out of ten bioclimatic variables were retained for the BSM procedure: three related to temperature (BIO2, BIO5, and BIO11) and three related to precipitation (BIO17, BIO18, and BIO19). Species Distribution Modelling framework The modelling approach was based on an ensemble procedure (Araújo and New 2007; Guisan et al. 2017). To better assess the influence of environmental variables on L. delicatula spatial distribution, bioclimatic and habitat variables (including digital elevation model) were modeled separately in two distinct ensemble models 272 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly (HSM and BSM), and their resulting maps combined as a spatial product, following the framework outlined by Ruzzier et al. (2024). We included DEM among the habitat variables to account for local microclimatic variations. This approach is based on the idea that integrating wide-scale (bioclimatic) and narrow-scale (habitat) variables into a single model may result in underestimating the contribution of variables acting at different geographic scales (scenopoetic vs local; sensu Hortal et al. 2010, see Franklin 2010). From a probabilistic perspective, this strategy relies on the assumption that bioclimatic suitability, i.e., the probability of finding an occurrence at a given site according to the BSM [P(BSM)], represents a necessary but not sufficient condition (a climatic constraint), while habitat suitability, as estimated by the HSM and intrinsically conditioned by bioclimatic factors [P(HSM|BSM)], refines the local probability. Accordingly, the conditional probability is obtained by multiplying the HSM and BSM probability (suitability) maps [P(BSM|HSM)=P(BSM)•P(HSM|BSM)] (see Grimmett and Stirzaker 2001). We acknowledge that the overall suitability map, obtained as the product of the BSM and HSM probability maps, does not represent a formal conditional probability but rather a composite suitability index, since the two models were developed independently. Nevertheless, because HSM is intrinsically conditioned by bioclimatic factors, the composite suitability index can be considered conceptually close to a conditional probability. When modelling presence-only data, the selection of pseudo-absences (PAs) or background points (BPs) is one of the most critical aspects to address (Barbet-Massin et al. 2012); pseudo-absences should be intended as locations where the species is assumed to be absent, though it was never confirmed through observation, while background points represent randomly sampled locations across the study area where the presence or absence of the species is unknown (Guisan et al. 2017). Sampling PAs or BPs exclusively within a species’ native range prevents a proper evaluation of the influence of bioclimatic variables, which primarily operate at broad spatial scales. As a result, the inability to isolate the effects of these bioclimatic factors would hinder the identification of key variables that limit the species’ native range and reduce the accuracy of predicting potentially suitable areas beyond that range. Therefore, when working at large spatial scales (e.g., continental or global), it is essential to sample PAs or BPs within the occurrences range and slightly beyond its limits in order to determine whether, and to what extent, bioclimatic variables constrain the species’ distribution (see Guisan et al. 2017; Adde et al. 2023). In contrast, variables such as land cover, DEM, and host plants exhibit fine-scale variation and act primarily at the local level, influencing habitat suitability. To assess the effects of these local-scale covariates, PAs or BPs should be selected as close as possible to known presence points (Fournier et al. 2017; Mateo et al. 2019; Adde et al. 2023). This approach minimizes the potential confounding effect of bioclimatic variables on the estimation of habitat-level influences. Since we empirically observed a relatively high density of occurrences of L. delicatula in urban areas, we suspected a possible sampling bias likely due to an easier detection of this species in urban context than in more remote areas (e.g., Beck et al. 2014). To evaluate this observation, we conducted an initial modelling attempt using BPs. As expected, this first output revealed a substantial overestimation of variables associated with anthropogenic contexts, particularly BU (results not shown). To overcome this issue, we decided to use a mixed set of BPs and PAs for HSM calibration, i.e. a random BPs identification combined with PAs identified on the base of a sampling bias map. We first identified an HSM calibration 273 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly area (within the narrow-range extent; Suppl. material 2: fig. S1) defined by surfaces where the density of occurrences exceeded 0.001 points/km2. Subsequently, we determined a set of PAs, in which the number of PAs equaled the number of occurrences. The allocation of PAs was proportional to the fractional cover of BU within each cell of the calibration area. Additionally, we identified a set of BPs, where the number of BPs equaled four times the number of occurrences within the calibration area. Overall, the total number of PAs + BPs for HSM was five times the number of occurrences (324 PAs + 1296 BPs). For the BSM, we identified BPs across the wide-range extent, setting the total number of BPs at five times the number of occurrences (2010 BPs; Suppl. material 2: fig. S2). Both HSM and BSM models return a map of suitability in terms of species presence probability. From the spatial product of the two maps of species presence probabilities, we obtained an overall suitability map, which returns the probability of finding the species in each location based on both habitat and bioclimatic covariates (see Ruzzier et al. 2024 for details). The modelling procedure was performed via the biomod2 package version 4.25-2 (Thuiller et al. 2024). As modelling algorithms, the Artificial Neural Network (ANN), the Flexible Discriminant Analysis (FDA), the Generalized Additive Model (GAM), the Generalized Boosting Model (GBM), Random Forest (RF), the Multiple Adaptive Regression Splines (MARS), and the Maximum Entropy (MAXNET) were used, all under “bigboss” settings (Thuiller et al. 2024). A total of 20 sets of PAs + BPs were used for HSM, and 5 sets of BPs were used for BSM. Given the lack of external data to evaluate the models’ accuracy, a 5-fold cross-validation data-splitting procedure was performed using 80% of the data as the training set with the remaining 20% as the validation set (Ruzzier et al. 2024). Each of the models obtained from both HSM and BSM was validated via two validation metrics: the area under the receiver operating characteristic (ROC) curve (AUC; Hanley and McNeil 1982) and the true skill statistic (TSS, Allouche et al. 2006). From the two metrics, we selected the best 10% models in validation to produce the final ensemble models, one for each metric (Allouche et al. 2006; Mandrekar 2010). Habitat and bioclimatic ensemble models were calculated via the “EMwmean” method in biomod2 package. To account for the performance of each single model used in building the ensemble model, we weighted its contribution by applying a ‘proportional’ option decay. A three-run permutation procedure was used to assess the importance of each variable included in the ensemble model. The same ensemble model framework was adopted for both the HSM and BSM. A total of 1,400 models were developed (7 modelling algorithms × 20 sets of (PAs + BPs) × 5 cross-validation runs × 2 validation metrics) for the HSM, and 350 models (7 modelling algorithms × 5 sets of BPs × 5 cross-validation runs × 2 validation metrics) for the BSM. The HSM and BSM were projected onto the following five geographic extents: 1) native area; 2) global, 3) Eastern and Southeast Asia (ESEA), 4) the United States of America (USA), and 5) Europe. The ESEA projection provides a detailed map of the native range and the closely invaded area (Japan, South Korea and Taiwan). The USA projection provides a map of areas already invaded by L. delicatula and those suitable for species establishment; we opted to project the calibrated model over the European territory since this species seems to be present in Spain (iNaturalist record) and because its potential impact on major European agricultural crops, particularly grapevine (Vitis spp.) (e.g., Huron et al. 2022; Urban and Leach 2023). 274 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly Based on all of the geographic extents, we produced the overall suitability map. Since the output of both HSM and BSM produced suitability maps expressed in terms of probability of occurrence of L. delicatula starting from independent datasets, the overall suitability maps were obtained as the product of the HSM and BSM maps. Using the Boyce index (Boyce et al. 2002), we assessed the effectiveness of the overall suitability map in predicting the distribution of L. delicatula within the calibration area. We also applied the index to evaluate the predictive performance of the overall suitability maps obtained by projecting the calibrated models onto invaded regions, including the eastern United States, where a major self-sustaining population is established, and South Korea. The Boyce index was used in an exploratory manner for the invasive populations, as these are still expanding and many areas identified as suitable by the model have not yet been colonized. This index measures a model’s ability to discriminate among areas of varying suitability by comparing predicted values at presence locations with those expected under a random distribution. Analyses were performed using the modEvA package (Barbosa et al. 2013). To detect possible extrapolation effects of the HSM and BSM at the global scale, we used the multivariate environmental similarity surface index (MESS; Elith et al. 2010), implemented in the dismo package (Hijmans et al. 2024). Resistance-Distance Constrained dispersal models To evaluate the invasion model of L. delicatula in the USA and South Korea, where multiple established populations are present, we developed a framework to assess the impact of environmental resistance on the dispersal capacity during the colonization process. This approach was not applied in Japan due to the recent establishment of this pest in that region. The framework relies on three steps: 1) the estimation of the annual dispersal distance; 2) the identification of the best function defining the resistance surface with respect to species colonization (i.e., species range expansion); 3) the reachable threshold, defined as cumulative cost faced by propagules during the colonization dispersal (this aspect accounts for both the ecological permeability of the environmental matrix and the intrinsic dispersal distance of L. delicatula). The invasion model was then used in forecasting potential invasion scenarios in invaded or potentially invadable areas across the world. Dispersal distance From a management perspective, besides estimating the annual range expansion, it could also be relevant to assess the magnitude of intrinsic linear dispersal of L. delicatula. To this end, from the overall global dataset we used occurrence data from the eastern USA (11,326 records), the only region worldwide where the invasion process is well documented, to estimate the annual range expansion and corresponding expansion rate. These data are consistent with those reported by De Bona et al. (2023). The annual dispersal distance was estimated by analyzing data from two consecutive years, under the assumption that occurrences were detected immediately following local population establishment. For each point in the dataset corresponding to year i, the closest point in the dataset from year i+1 was identified, and the Eu- 275 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly clidean distance between them was calculated. Starting from 2015 (when more than one occurrence was available), we selected and bonded among consecutive years all the closest pairwise Euclidean distances. Then, we bootstrapped the overall median of these Euclidean distances, along with its confidence intervals at 95% (Efron and Tibshirani 1993). In addition, we bootstrapped the 99th percentile of the closest pairwise Euclidean distances, and its confidence intervals at 95%. This procedure allowed us to estimate the maximum dispersal distance for L. delicatula from one year to another. All points from the year i+1, where the value of the closest pairwise Euclidean distances exceeded the 99th percentile, were classified as an event of new introduction due to human-mediated translocation (see “Resistance cost threshold”). Resistance surfaces Resistance surfaces are widely used to assess how much environmental features hinder animal dispersal, to both assess the population ecological connectivity (Bani et al. 2015; Lecis et al. 2022, 2024) and to evaluate the colonization processes (Sherpa et al. 2020). To achieve this, we calculated a set of resistance maps as function of the overall suitability map obtained from the SDM framework. We tested two methods to evaluate the resistance encountered by propagules during dispersal: a) resistance defined as the complementary value to suitability (i.e., additive inverse relationship; Zeller et al. 2012; Bani et al. 2018), and b) resistance that varies non-linearly with suitability, following a logistic function (McRae and Beier 2007; Zeller et al. 2012): a) R = 1−S b) R = 1 / (1 + exp(−k * (1−S−x0))) where R is the resistance, S is the suitability, k is the steepness and x0 is the inflection point that determines a change of concavity. In this case we used three sets of values with different values of x0 (0.25, 0.5, 0.75), while we set k value equal to 15 (Fig. 1). The rationale behind this concept is that colonization begins by prioritizing areas with the lowest resistance, those with the highest suitability (i.e., highest probability of colonization), even if they are geographically farther away but still within the dispersal distance achievable by propagules. Colonization then extends to regions with lower or suboptimal suitability only after these high-suitability areas have been saturated (Lovell et al. 2021). In addition, we assumed that newly colonized cells immediately act as new sources for propagule spread. We then determined the value of resistance between two neighboring cells of the resistance map, producing a transition matrix, through the function “transition” of the gdistance R package (van Etten 2017); in the transition matrix, for a pair of neighboring cells A and B the transition value between them (i.e., the ecological resistance experienced by an organism passing from the cell A to cell B) depends on the resistance values of the two cells and the transition function defining the mathematical relationship between the two cell values; as the transition function we used the reciprocal (i.e., the inverse relationship; Adriaensen et al. 2003; McRae and Beier 2007) of the average value of the two adjacent cells: f(x)=1/mean x, where x is a vector containing the values of two neighboring cells. The spatial distance between two cells is also considered: for a 3×3 group of cells, the distance between 282 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly cost threshold. From just under 2,000 km2 colonized in 2015, the invaded area in 2024 had expanded approximately 70-fold, reaching nearly 130,000 km2. These areas include only occurrences associated with propagules presumed to originate from the initial introduction (i.e., the area enclosed by the isochores, as illustrated in Fig. 8). Table 1. Values of attributes obtained from the functions used in assessing the expansion model of Lycorma delicatula in the eastern USA. In italics are shown the best thresholds selected by the score for each function and related parameters used in defining the resistance surface. In bold italics: the best threshold associated with the best function. Functions Resistance thresholds and spatial attributes Values of the resistance thresholds and derived spatial attributes Linear Resistance cost threshold 10,500 11,000 11,500 12,000 12,500 13,000 Invaded range (km2) 108,418 116,888 125,280 134,093 141,776 149,750 Point included into the range (%) [A] 78.2 82.3 86.4 89.3 91.2 92.9 Points/range (points/100km2) [B] 6.47 6.32 6.19 5.98 5.77 5.57 Score [A×B] 506.5 520.3 534.1 533.7 526.4 517.2 Logistic: Resistance cost threshold 17,000 17,250 17,500 17,750 18,000 18,250 k = 15; x0 = 0.25 Invaded range (km2) 133,350 136,314 139,363 142,570 145,342 148,463 Point included into the range (%) [A] 85.5 86.5 87.9 88.7 89.3 89.9 Points/range (points/100km2) [B] 5.77 5.70 5.66 5.58 5.51 5.44 Score [A×B] 492.2 492.9 497.4 494.8 492.4 488.8 Logistic: Resistance cost threshold 1,100 1,150 1,200 1,250 1,300 1,350 k = 15; x0 = 0.50 Invaded range (km2) 95,782 98,538 101,283 101,283 101,283 108,114 Point included into the range (%) [A] 89.6 91.6 93.1 93.9 94.4 94.8 Points/range (points/100km2) [B] 5.74 5.87 5.97 6.02 6.05 6.08 Score [A×B] 537.1 546.0 548.4 545.4 540.0 532.8 Logistic: Resistance cost threshold 25 30 35 40 45 50 k = 15; x0 = 0.75 Invaded range (km2) 98,368 117,700 140,666 163,061 181,944 201,574 Point included into the range (%) [A] 68.8 83.8 91.5 93.1 93.4 93.6 Points/range (points/100km2) [B] 6.27 6.39 5.83 5.12 4.61 4.17 Score [A×B] 431.5 534.8 533.6 477.0 430.6 389.8 Table 2. Area invaded by L. delicatula in eastern USA from 2015 to 2024 starting from its initial introduction point. The range area was estimated by the range expansion model and refers solely to the year-by-year progression of the populations derived from propagules whose origin is putatively attributed to the initial introduction (it does not include areas occupied by populations believed to have originated after internal translocations mediated by humans). Year Invaded area (km2) 2015 1,797 2016 8,611 2017 21,041 2018 35,724 2019 52,363 2020 72,511 2021 89,480 2022 102,769 2023 115,525 2024 127,721 283 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly Figure 8. Range expansion model: year-by-year 2017 to 2024 invasion dynamic of L. delicatula in eastern USA starting from first introduction. The isochore lines derived from the range expansion model refer solely to the population derived from propagules whose origin is putatively attributed to the initial introduction (it does not include areas occupied by populations believed to have originated from internal translocations mediated by humans). Red points represent the occurrences that are putatively assigned to the population originating from the first introduction event, based on the maximum annual dispersal distance estimated for L. delicatula (~25 km). Future invasion scenarios in the USA, South Korea and Japan, and potential invasion scenario in Europe The expansion scenarios for L. delicatula over the next 10 years in the USA, South Korea and Japan, assessed by the range expansion model and starting from its current distribution, are shown in Figs 9–11, respectively. The European scenario, based on three hypothetical entry points and one iNaturalist recording point for L. delicatula (Madrid, Spain), and developed using the eastern USA range expansion model, is illustrated in Fig. 12. This invasion scenario is based on the assumption that population expansion requires prior local establishment. Based on the eastern USA range expansion model, the potential spread of L. delicatula in the Iberian Peninsula is unlikely due to the low suitability of the introduction area. 284 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly Figure 9. Future invasion scenario of Lycorma delicatula in the USA and Southern Canada over the next 10 years, based on the species’ current distribution. Figure 10. Future invasion scenario of Lycorma delicatula in South Korea over the next 10 years, based on the species’ current distribution. All these scenarios should be considered without accounting for potential human-mediated translocations, which could further facilitate the expansion of L. delicatula range. 285 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly Figure 11. Future invasion scenario of Lycorma delicatula in Japan over the next 10 years, based on the species’ current distribution. Figure 12. Possible European invasion scenario of Lycorma delicatula, assuming three hypothetical entry points 10 years after local population establishment from Milano Malpensa Airport (Italy), Paris Charles de Gaulle Airport (France), and Hamburg Port (Germany). Red points indicate putative occurrences of L. delicatula. 286 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly Discussion Species distribution model: comparative evaluation Although the output of HSM built solely on habitat variables may be partially influenced by bioclimate (host plant taxa and land covers can be affected by climatic conditions to some extent), bioclimatic variables used to build the independent BSM can nonetheless provide a separate contribution to the distribution of L. delicatula, distinct from that offered by habitat variables. Lycorma delicatula seems to be primarily influenced by bioclimatic conditions in its native range. However, within areas that are bioclimatically suitable, habitat availability appears to be a key limiting factor in the species distribution, as indicated by the habitat suitability map. When considering the overall suitability map, the combined impact of both bioclimatic and habitat factors on the species distribution pattern becomes even more evident. Consequently, it seems unlikely that the species could naturally expand where all its ecological needs are not simultaneously met, as can be observed by the overall suitability map within native range. In addition, within its native range, L. delicatula seems unable to expand into nearby areas with suitable habitats and bioclimates when ecological barriers, such as oceans or major mountain ranges, are present. Built-up areas (BU) emerged as an important variable in the calibration area (i.e., the native range), showing a positive effect on the probability of species presence. This effect is also confirmed by other authors, who have observed a similar influence in areas subject to invasion (Owen et al. 2024; Levine et al. 2025). According to these studies, L. delicatula may benefit from the urban environment, particularly due to the presence of host plants used for ornamental purposes, as well as the microclimatic “heat island” effect, which could increase the chances of survival for the species. In addition, the mean temperature of the coldest quarter (BIO11) significantly influences the winter survival of L. delicatula, as its eggs require a period of overwintering to break diapause. However, prolonged and extreme cold events can lead to embryo mortality, as previously observed by several authors (Lee et al. 2011, 2014; Keena and Nielsen 2021; Kreitman et al. 2021). The numerous non-native areas with similar bioclimatic and habitat conditions identified as suitable by our modelling framework can thus be colonized by the species only through human-mediated translocation, due to geographic separation (Ladin et al. 2023). The adopted approach appears to effectively identify the most suitable areas for L. delicatula on a global scale, as all currently invaded areas are classified as suitable or, at least, moderately suitable both in terms of bioclimatic factors and habitat, including host plants (overall probability of species presence ≥ 0.5). For this reason, it can be hypothesized that L. delicatula may not be sufficiently plastic to change or expand its ecological niche. Our results, like several recent studies (Namgung et al. 2020; Wakie et al. 2020; Huron et al. 2022; Chartois et al. 2024; Zhao et al. 2024; Lin and Liao 2025), differ from Jung et al. (2017), who linked Lycorma delicatula mainly to tropical and subtropical humid climates. In contrast, most studies, including ours, associate the species with a broader range from humid continental to humid subtropical climates. These differences likely stem from Jung et al. use of the CLIMEX model, while others employed algorithms like Maxent, Random Forest and others. 287 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly Our study offers important advancements compared to previous analyses in modeling the distribution of L. delicatula. Unlike Namgung et al. (2020), who integrated both bioclimatic and habitat variables (limited to Vitis vinifera L. as the host plant) within a single model focused solely on South Korea, our approach models these factors separately, allowing for a more detailed and differentiated assessment. Furthermore, we utilize a composite layer comprising several dozen host plant taxa, in contrast to the single host considered by Namgung et al., and apply an ensemble modeling approach, enhancing prediction robustness. Finally, the allocation of PAs or BPs likely influenced the discrepancies observed between the two studies, as it is well known that modelling techniques based solely on presence data are highly affected by this aspect (as noted by Steen et al. 2024). Wakie et al. (2020) developed a bioclimatic model incorporating presence data from both the native range and invaded regions, with background points randomly distributed within a large polygon surrounding the occurrence data. Their results are generally consistent with our BSM, except for the native range, where they predict notably lower suitability. Wakie et al. acknowledge that heavy reliance on non-native data in their model likely contributes to the underestimation of suitability in the native range. Our BSM map largely agrees with that of Huron et al. (2022), which was produced by averaging two different maps: the first based on bioclimatic variables, digital elevation model, forest canopy height, and access to cities, and the second related to the suitability of Ailanthus altissima as a host species. The alignment between our BSM map and the overall suitability map by Huron et al. (2022) is largely due to the dominant influence of the bioclimatic covariates in the model that produced the latter map provided by Huron and co-authors. The suitability map produced by Zhao et al. (2024), based on a Maxent model integrating 19 bioclimatic variables and a digital elevation model, and using occurrence data from both native and invaded regions, largely aligns with our BSM map in terms of global suitability patterns. However, their map generally shows lower suitability values compared to ours. A key limitation of their approach is the random global allocation of background points, which risks placing these points in climatically suitable areas that remain uncolonized due to biogeographic barriers or the recent absence of introductions. This leads to an underestimation of the importance of variables influencing true habitat suitability (see Ruzzier et al. 2025; Steen et al. 2024). Nevertheless, since most suitable regions are already partially or fully invaded by L. delicatula, this issue likely had minimal impact on their overall results, as most background points fell within unsuitable areas. Our approach, however, effectively minimizes this caveat. Chartois et al. (2024) used an ensemble model based on bioclimatic variables to produce a global suitability map largely consistent with our BSM and Zhao et al. (2024). However, while Zhao et al.’s map shows generally lower suitability than ours, Chartois et al.’s map indicates higher suitability in both intensity and geographic extent. Their study properly allocated PAs for bioclimatic modeling but did not exclude regions separated by biogeographic barriers, which may be suitable but uncolonized. The key difference with our BSM likely stems from their use of PAs restricted to the north and south of the presence polygon along a latitudinal gradient, whereas we allocated background points (BPs) across the entire gradient. This comprehensive BP allocation in our model enables a more balanced assessment of bioclimatic effects, while their method may overemphasize conditions 288 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly within the presence polygon by assuming it is entirely suitable and areas outside are unsuitable. The model we developed did not identify the U.S. state of California as a critical area for potential invasion risk, despite some detections of L. delicatula reported in literature and databases. We believe these occurrences likely represent adventitious populations, as they have not shown any signs of local range expansion. Such occurrences may have been sustained by commercial transportation (see Jones et al. 2022; Ladin et al. 2023; Montgomery et al. 2023) within a region of low suitability. Conversely, other studies published around the same time (Wakie et al. 2020; Jones et al. 2022; Montgomery et al. 2023) identified California as a suitable area for this pest, based on modelling approaches that included occurrences from populations not yet at equilibrium. However, we argue that calibrating models with occurrences from non-equilibrium populations, such as adventive populations, may distort niche estimation and the influence of covariates, potentially leading to misleading predictions. Similarly, propagules detected in central U.S. states and Florida are in areas of low suitability and thus represent a low invasion risk. The same pattern applies to Canada, where confirmed occurrences occur mostly in areas of low suitability, except for Nova Scotia, the western shores of Lake Ontario, and the northwestern shores of Lake Erie. Ultimately, Namgung et al. (2020), Huron et al. (2022), Jones et al. (2022), Montgomery et al. (2023), and Chartois et al. (2024) should be credited for acknowledging the importance of incorporating or assessing host plants, albeit through different approaches, when evaluating the suitability for potential colonization by L. delicatula. The analytical framework we have developed is, in our view, capable of providing a more robust predictive capacity, thanks to a series of methodological aspects, many of which represent novel implementations. In particular, the calibration of the model using data from populations in equilibrium within the native range of the species allows for a more realistic assessment of the ecological niche occupied by the species. Furthermore, the criteria adopted to identify BPs and PAs, along with their spatial extent and relationship to occurrence records, play a crucial role in refining the estimation of the influence of environmental covariates. The integration of host plant distributions, combined with high-resolution land cover data, further enhances the estimation of habitat-level suitability, providing a more accurate representation of the ecological requirements of the species. Additionally, the separate implementation of a HSM and a BSM, within a two-step modelling approach, enables a more effective distinction between the effects of local-scale covariates (addressed in the HSM) and those associated with broader-scale ecological processes (addressed in the BSM). Altogether, this modelling strategy leads to an improved explanatory understanding of the distribution of L. delicatula within the native range and delivers more effective predictive performance, in line with a niche conservation perspective, in regions outside the native area. Evaluation of SDM projections The Boyce index, employed to assess the effectiveness of the SDM calibrated on the native range in predicting the distribution of L. delicatula in invaded regions (USA and South Korea), reveals two distinct scenarios. In the USA, although this insect has yet to saturate suitable habitats (with populations still expanding), the 289 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly evaluation is highly positive, as a significant proportion of occurrences fall within areas identified as suitable. In contrast, the Boyce index for South Korea indicates that the distribution of occurrences is far from saturating suitable areas and does not exhibit a pattern consistent with natural population expansion. Instead, the occurrences appear highly scattered, likely resulting from multiple introduction events and/or human-mediated translocations, which have initiated or are sustaining various invasion hotspots or may represent adventive populations. Given the recent invasion of Japan by L. delicatula, evaluating the accuracy of the model prediction for this region no longer seems appropriate. Despite the presence of several suitable areas on a global scale, many of these areas, which are currently not invaded, appear to be potentially susceptible to colonization, in agreement with Huron et al. (2022). Several of these areas are characterized by the presence of valuable crops on which this pest could cause a significant economic damage (e.g., grapevine). These areas have a mild climate, with temperatures higher than −11 °C during winter (the minimum survival threshold for the insect, as indicated by Lee et al. 2014), and are also human-dominated, with an extensive presence of non-native plants, both invasive and ornamental, used as host plants by L. delicatula (Nixon et al. 2023a). These conditions thus greatly support L. delicatula populations, creating suitable habitats and facilitating the expansion of this pest by forming invasion corridors, in agreement with Zhang et al. (2023). Considering the invasion process of L. delicatula in ESEA, South Korea is already a largely invaded territory where it can be considered nearly ubiquitous, even at present, with populations ranging from adventive to established. In Taiwan, the risk of invasion is rather limited due to the low suitability across most of the territory, except for the anthropized context of the northern part of the island (e.g., Taipei and the surrounding areas), where adventive populations already exist. In contrast, for Japan, considering the high suitability of much of the territory and the presence of established populations, the risk of invasion appears significant. All these potential risks of invasion should be interpreted in the light of the current climatic and habitat conditions. Range expansion model and invasion scenarios Considering both the suitability and the intrinsic dispersal potential of L. delicatula, we simulated the insect’s range expansion in invaded areas, starting from the current distribution, as well as in potentially invadable areas, over a 10-year period. Our projections for the United States, where invasive populations are more widespread, identify two major regions at risk of colonization. The first is located east of the Appalachian Mountains, extending toward the Atlantic coast, and appears to be already extensively colonized, with L. delicatula having nearly occupied all suitable habitats. The second lies further inland, west of the Appalachians and south of the Great Lakes, where the invasion is ongoing, showing clear directional movement toward the continental interior. Smaller, likely internally relocated populations are already present there, and the model predicts nearly complete colonization of the most suitable areas within the next 10 years. The model also indicates a low likelihood of populations in California, the Central-Western states, and Florida developing into expanding invasive populations. In South Korea, the 10-year expansion scenario projected from currently colonized areas suggests that potential invasion will likely concentrate along the south- 290 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly western coast. Other populations, possibly arising from internal translocations, are unlikely to contribute to further range expansion, including in the Seoul area, despite high occurrence densities recorded there. Given these uncertainties, this invasion pattern should be interpreted with caution, but we deemed it valuable to present this potential scenario in relation to habitat suitability across South Korea. For Japan, where suitability is high, the 10-year scenario indicates that current populations could rapidly spread through the central and southern regions of Honshu Island. However, the model does not predict full saturation of suitable areas further north within this period, so additional range expansion is still possible. In Europe, the range expansion simulation emphasizes that successful invasion depends on established populations within medium to highly suitable areas. This is illustrated by the single recorded occurrence in Madrid, Spain, which our model suggests is unlikely to contribute to further spread. Conversely, introductions and establishments in highly suitable areas could lead to rapid expansion, as demonstrated by the three hypothetical entry points included in our simulations. Pest status perspectives The findings of this study suggest that, given the limited ecological plasticity and relatively slow expansion rate estimated for L. delicatula, the prompt implementation of eradication measures could be particularly effective in eliminating adventive populations. These efforts must be reinforced by stringent measures to prevent the establishment of new invasion hotspots caused by accidental human-mediated dispersal. Additionally, in moderately to highly suitable areas, further efforts are required, including the implementation of effective monitoring protocols and eradication methods as soon as the species is detected. The damage caused by L. delicatula invasion should be assessed considering its potential impact on crops and ornamental plants, derived from its feeding activity and honeydew deposition, especially when it aggregates in huge numbers (Urban and Leach 2023). As highlighted by the suitability model and the invasion scenarios, a substantial part of the USA remains exposed to the risk of invasion with potentially similar impacts to those observed in other invaded areas, as indicated by Harper et al. (2019), Leach and Leach (2020), Urban (2020), Urban and Leach (2023), and Molfini et al. (2024). Global trade increases the likelihood of L. delicatula introductions and associated economic losses in countries with extensive suitable areas for this species, particularly key producers of agricultural and viticultural products, as emphasized by Huron et al. (2022). This scenario appears especially realistic for Europe, given the recent record of this pest in Spain. Our findings also align with broader invasion ecology principles and with evidence from other invasive herbivores. In many species, bioclimatic conditions determine the broad limits of potential spread, while habitat and host availability act as finer-scale filters. For example, the spread of the Spongy Moth, Lymantria dispar (Linnaeus, 1758), across North America depends not only on climatic suitability but also on the distribution of host tree species (Régnière et al. 2012), and the Emerald Ash Borer, Agrilus planipennis Fairmaire, 1888, remains tightly constrained by the availability of Fraxinus hosts despite suitable climates (Herms and McCullough 2014). Similarly, bark beetle, Dendroctonus ponderosae (Hopkins, 1902), dynamics 291 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly (illustrate how climate defines physiological boundaries, while local forest composition governs outbreak severity (Bentz et al. 2010). These parallels underline that invasion risk of L. delicatula should not be evaluated on climatic suitability alone, but on the combined filtering action of climate, habitat, and dispersal opportunities. Climatic conditions set the physiological limits for survival and development, while the presence and density of suitable host plants affect successful colonization and reproductive output. Local habitat features and composition further regulate population growth and expansion, ultimately shaping the speed and pattern of invasion. Recognizing the interplay of these factors is therefore essential to fully understand the mechanisms underlying L. delicatula spread and to accurately assess invasion risk. Our integrated HSM–BSM framework reflects this multi-layered filtering process, emphasizing that climatically suitable areas may not be colonized unless appropriate hosts and habitats are present, and that natural spread remains limited by major geographic barriers unless overcome by human-mediated translocation. As predicted by our scenario, the introduction of L. delicatula into multiple suitable areas across the European subcontinent could trigger a rapid pest expansion, possibly resulting in significant impacts on crop production. Countries like France, Italy, and Germany, among the world’s largest wine producers, are particularly vulnerable to the invasion risk of L. delicatula, as a significant portion of their territory is highly suitable for this pest. The accidental introduction of L. delicatula in these regions could indeed cause severe damage to the wine production sector, leading to substantial losses as documented in Pennsylvania, USA (Leach and Leach 2020), where, in 2018, up to 90% of grapes were lost despite the use of insecticides (Urban 2020). The potential impact on the fruit production sector should also not be underestimated. However, these impacts may depend on the features of the agroecosystem in question, including management strategies, as discussed in Nixon et al. (2023b) and Urban and Leach (2023). Conclusions The SDM and the spread model implemented in this study represent valuable tools for managing L. delicatula and invasive species more broadly, as they enable the precise identification of areas at highest risk of invasion, including both colonization and expansion phases. This would support the development and implementation of early surveillance and targeted monitoring strategies, optimizing resource allocation and enhancing the effectiveness of interventions, while reducing unnecessary costs in low-risk areas. Moreover, mapping environmental suitability is crucial for effectively planning potential introductions of biological control agents, by identifying areas where such actions are most likely to succeed. The outputs of these models also provide a solid basis to support political and management decisions, offering useful insights for defining quarantine zones, regulating the movement of goods, or launching outreach and awareness campaigns for local communities. Looking ahead, integrating these models with socio-economic data, such as transportation flows or commercial activities, could further improve predictive accuracy by more precisely identifying human-mediated dispersal routes and potential entry points. Finally, extending these approaches to future climate scenarios would allow anticipation of invasion dynamics over the medium to long term, enabling proactive and adaptive planning of management strategies. 298 NeoBiota 103: 267–298 (2025), DOI: 10.3897/neobiota.103.154246 Enrico Ruzzier et al.: Global invasion scenarios of Spotted Lanternfly Supplementary material 1 Supplementary tables Authors: Enrico Ruzzier, Davide Scaccini, Pietro Tirozzi, Valerio Orioli, Olivia Dondina, Andrea Di Giulio, Alberto Pozzebon, Luciano Bani Data type: xlsx Explanation note: table S1. Schematic associations between the host plant taxon and land cover. table S2. Evaluation HSM. Table of the three evaluation metrics used for HSM. table S3. Evaluation BSM. Table of the three evaluation metrics used for BSM. table S4. Importance of the covariates included in the HSM. table S5. Importance of the covariates included in the BSM. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.103.154246.suppl1 Supplementary material 2 Supplementary figures Authors: Enrico Ruzzier, Davide Scaccini, Pietro Tirozzi, Valerio Orioli, Olivia Dondina, Andrea Di Giulio, Alberto Pozzebon, Luciano Bani Data type: pdf Explanation note: fig. S1. Distribution of the 20 sets of PAs + BPs used for the habitat suitability model (HSM) calibration in the narrow native calibration range. fig. S2. Distribution of the 5 sets of BPs used for the bioclimatic suitability model (BSM) calibration in the wide calibration range. fig. S3. Response curves of covariates included in the Lycorma delicatula habitat suitability models (HSM). fig. S4. Response curves of covariates included in the Lycorma delicatula bioclimatic suitability models (BSM). fig. S5. Habitat suitability map for Lycorma delicatula in the narrow native calibration range obtained projecting the habitat suitability model (HSM). fig. S6. Bioclimatic suitability map for Lycorma delicatula in the wide calibration range obtained projecting the bioclimatic suitability model (BSM). fig. S7. Overall suitability map for Lycorma delicatula in the narrow native calibration range. fig. S8. Boyce of the overall suitability map for Lycorma delicatula in the narrow native calibration range. fig. S9. Global (a) HSM and (b) BSM Multivariate Environmental Similarity Surfaces (MESS) index. fig. S10. Overall suitability map for Lycorma delicatula in the Eastern USA. In red species occurrences. fig. S11. Boyce of the overall suitability map for Lycorma delicatula in the Eastern USA. fig. S12. Overall suitability map for Lycorma delicatula in South Korea. In red species occurrences. fig. S13. Boyce of the overall suitability map for Lycorma delicatula in South Korea. fig. S14. World projections of Lycorma delicatula suitability maps; a) habitat suitability map; b) bioclimatic suitability map. fig. S15. Eastern and southeast Asia (ESA) projections of Lycorma delicatula suitability maps. fig. S16. USA projections of Lycorma delicatula suitability maps. fig. S17. Details of USA (East Coast) projections of Lycorma delicatula suitability maps. fig. S18. European projections of Lycorma delicatula suitability maps. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.103.154246.suppl2