Geographic records of distribution of H. halys were retrieved from Global Biodiversity Information Facility (GBIF.org 2022). Geographic records of T. japonicus and T. mitsukurii were retrieved from GPS latitude and longitude coordinates of: i) monitoring sites of DISAFA (University of Turin, Italy) and DAFNAE (University of Padova, Italy), ii) available data on Hymenoptera Online database (HOL 2021), and iii) literature data (Suppl. material 1: table S1). Localities lacking geographic coordinates were georeferenced using Google Maps. Records with doubtful species identification, unspecified or unknown localities were deleted. In localities where T. japonicus was used in augmentative release programs as a biological control agent (BCA) of H. halys, data for a year were considered valid if confirmed by observations during the following spring-summer, indicating wasp overwintering after the release. These various sources provided a total of 14,489, 393, and 356 point locations recorded for H. halys, T. japonicus and T. mitsukurii, respectively (Suppl. material 4: fig. S1).

In order to verify the absence of potential divergences in the climatic niches of populations in primary and secondary ranges, which are sometimes introduced by bottleneck effects in introduced populations, we analyzed the climatic preference of the three species across their distribution ranges (Fig. 1). The results suggested that we should try to model the distribution of the three species using the entire global dataset rather than considering native and invaded ranges separately.

Figure 1.

Climatic preferences of Halyomorpha halys (A), Trissolcus japonicus (B), and T. mitsukurii (C) in the native range (i.e., Asia) and invaded range. Colored dots represent the climatic conditions in the species presence sites and gray shade indicates the climatic conditions available across the species distribution.

We used the distribution data for modeling habitat suitability for H. halys, T. japonicus, and T. mitsukurii at three spatial scales. All three species’ global data were used to fit the models at the three scales. These models were fitted with three different spatial resolutions and projected into three different geographic areas. The maps were obtained from data combined and processed at different resolutions. Therefore, the global, European and Italian scale maps are the result of different processing. We used environmental predictors with a resolution of 10 minutes of geographic degrees (i.e., pixels about 20 km large) for predicting the habitat suitability at the global scale, 5 minutes (i.e., about 10 km) for the European scale, and 2.5 minutes (i.e., about 5 km) for the Italian scale. In order to mitigate the geographical bias associated with the non-random process of data collection, presence data were thinned to remove duplicate points in the pixels (Verbruggen et al. 2013). This multi-scale approach provided practical recommendations for managing the H. halys invasion and simultaneously overcoming the ‘tyranny of the local’ (Groves 2003), which consists of overlooking important areas at the large scale when reserves are planned locally.

We used several parameters representing climate, land morphology, land cover, and water availability as predictor variables. Climate variables were downloaded from the WorldClim 2.1 databank (Fick and Hijmans 2017). We selected a subset of non-correlated climatic variables for fitting the models by calculating the variance inflation factor (VIF) (Belsley 1991) and excluding those with VIF ≥ 5 (Dormann et al. 2013). Thus, we used eight, seven, and eight climatic variables at global, continental, and national scales, respectively. Land morphology variables were calculated from the altitude information associated to the WorldClim climatic data (Fick and Hijmans 2017). In particular, we calculated slope, roughness, eastness (calculated as the sine of terrain aspect), and northness (cosine of terrain aspect) (Kumar et al. 2006) for each pixel by comparing elevation values of the pixel with those of the eight surrounding pixels. The land cover was downloaded from the GLC2000 database (European Commission Joint Research Centre 2003), which is available at the resolution of 1 km at Equator. Data on water availability were extracted from the Digital Global Map of Irrigation Areas (Siebert et al. 2013). In addition, for the two parasitoids, we incorporated the effect of the interaction with H. halys by adding the suitability for the stink bug as a predictor. Each variable was resampled at the three resolutions, tested for collinearity (see Suppl. material 2: table S2), and used for fitting and projecting the models.

The habitat suitability for the three species at the different scales was calculated using the MaxEnt algorithm (Phillips et al. 2004, 2006). MaxEnt minimizes the relative entropy between the probability densities in the covariate space estimated from the presence data and from the background (Elith et al. 2011). In order to correctly define unsuitable conditions, we defined the background conditions by selecting a number of random points across the study area (i.e., the globe) equal to five times the number of presence points. In order to account for the uncertainty in the definition of the background conditions associated with the randomicity of the background points, we repeated the analysis 100 times and generated a final consensus calculating the average value from the 100 models. All the models were 10-fold cross validated by deriving the Area Under the ROC Curve (AUC) value (Fielding and Bell 1997; Faraggi and Reiser 2002) and the Boyce index (Boyce et al. 2002). The AUC value can be interpreted as the probability that a presence site, randomly chosen from the dataset, has a higher predicted value than an absence site (Elith et al. 2006). The Boyce index measures how much model predictions differ from random distribution of the observed presences across the prediction gradients and is an appropriate metric for presence-only models (Hirzel et al. 2006).

To derive discrete categories from the continuous values of habitat suitability, we used a data-driven multi-thresholds approach (Nenzén and Araújo 2011). For each model we applied five different binarization criteria. First, the highest threshold at which there is no omission (Cantor et al. 1999). Second, the threshold at which the sum of the sensitivity (true positive rate) and specificity (true negative rate) is highest (Cantor et al. 1999). Third, the threshold at which sensitivity and specificity are equal (Fielding and Bell 1997). Fourth, the threshold at which the modeled prevalence is closest to the observed prevalence (Dormann et al. 2008). Fifth, the threshold at which the sensitivity of the model is equal to 0.5 (Cantor et al. 1999). These thresholds were used for dividing the suitability values in six categories of increasing suitability, from no suitability (not suitable according to all five criteria) to high suitability (suitable according to all five criteria).

In order to explore the geographic dimension of the interaction between H. halys and its parasitoids, we generated the ‘host-parasitoids co-suitability maps’ by overlapping the habitat suitability models of the three species. In these maps, pixels are categorized based on their contemporary suitability for H. halys and for the parasitoids. We defined: (1) areas with scarce suitability for all the three species, (2) areas with scarce suitability for H. halys but good suitability for one or (3) two parasitoids, (4) areas with good suitability for H. halys but scarce suitability for the two parasitoids, and (5) areas with good suitability for H. halys and for one or (6) two parasitoids. In addition, for clarifying which parasitoid has the highest probability of success in the different areas, we produced conflict maps by overlapping the habitat suitability of the two Trissolcus species. In particular, we defined areas with scarce suitability for both the two species, areas with good suitability for one species and scarce for the other, and areas with good suitability for both the species. In the host-parasitoids and in the parasitoids co-suitability maps, we defined pixels with scarce suitability as those considered suitable according to less than three criteria and defined good suitability pixels as those considered suitable according to three criteria or more. In addition, we summarized the overlaps at the three scales by continent/country/region as the percentage of surface in the different conditions. All the analyses, as well as these maps, were done in R (R Core Team 2021), using the packages dismo (Hijmans et al. 2021), raster (Hijmans 2022), rgdal (Bivand et al. 2022), rJava (Urbanek 2021), and usdm (Naimi et al. 2014).

A fundamental phase in species distribution modelling is the validation of the outputs, which measures how accurately the model predicts the presence or absence of the species. All our models have been 10-fold cross-validated; this means that each model was fitted 10 times with 90% of the available data, randomly chosen at each repetition, and used the remaining 10% as pseudo-independent data for validating the predictions. To do that, the true positive and the false positive rates (i.e., how many presence or absence points are correctly predicted) were measured for different thresholds, and the definite integral of this curve was calculated. This value represents the AUC index, which is one of the most frequently used methods for model validation and varies between 0 (no test point is correctly predicted) and 1 (100% of the test points are correctly predicted). In addition to this, we also calculated the Boyce index, which is an appropriate metric for presence-only models and varies between -1 and +1, with positive values indicating predictions that are consistent with the distribution of presences in the evaluation dataset (Hirzel et al. 2006). Our models obtained, on average, an AUC score of 0.987 (± 0.01) and a Boyce index of 0.994 (± 0.09), these results indicate an extremely high ability to predict the species presence/absence correctly.

While the GBIF dataset of H. halys can be considered robust, thanks to a large number of reliable records (citizens usually identified this stink bug with good accuracy), the distributional data of T. japonicus and T. mitsukurii needed a careful review due to the high risk of misidentification. Avoiding erroneous data is, therefore, crucial. Consequently, we performed a detailed screening of the bibliographic records (Suppl. material 1: table S1), also taking into account the historical taxonomic complications presented by Talamas et al. (2017) and Buffington et al. (2018). The dataset was supplemented by many directly verified field collections.

As expected, the model of potential distribution of H. halys here proposed is similar to models previously estimated using both MaxEnt (Zhu et al. 2012) and CLIMEX (Kriticos et al. 2017; Yonow et al. 2021). At a global scale, our model shows lower suitability for expansion of this pest to tropical regions (Fig. 2A) than reported by using CLIMEX. At the European scale, compared with lower climate suitability proposed by Kriticos et al. (2017) and Yonow et al. (2021), our model shows that the habitat suitability for H. halys largely reaches a medium-high level (Fig. 2D) with patterns of high suitability in eastern France, northern Italy, Slovenia, and Croatia, as confirmed by several field studies (Hemala and Kment 2017; Streito et al. 2021; Delbac et al. 2022; Hess et al. 2022). In Italy, our model shows high potential suitability and distribution of the pest mainly in the north, aligning with observations of economic damage from previous studies (Bosco et al. 2018; Moraglio et al. 2018; Francati et al. 2021). The maps produced at different scales are not simple zooms of the same model. Rather, they are different models, fitted and projected at different resolutions. Thus, they used differently aggregated data and, as a consequence, their outputs are not perfectly identical. Nevertheless, the general pattern is congruent among different models.

Avila and Charles (2018) and Yonow et al. (2021) made assumptions on the potential distribution of T. japonicus – also supported by data on the wasp distribution in North America and Europe (Talamas et al. 2015; Milnes et al. 2016; Sabbatini Peverieri et al. 2018; Abram et al. 2019; Stahl et al. 2019b) – resulting in a similar finding as reported by our study, except for Spain and the entire Mediterranean coast, where our model indicates low or no suitability.

According to our results, T. mitsukurii has a broad, but low, suitability in Australia, New Zealand, North and South America (Fig. 2C). Nevertheless, the potential suitable extent we predicted is more limited with respect to the assumptions provided in a previous study (Yonow et al. 2021). The rapid expansion of adventive populations of T. mitsukurii in Europe (Sabbatini Peverieri et al. 2018; Scaccini et al. 2020; Bout et al. 2021; Rot et al. 2021; Mele et al. 2022) confirmed our hypothesis of a high level of potential habitat suitability: this species has the potential to spread throughout Europe and beyond.

In Italy, T. japonicus shows greater potential for range expansion at lower latitudes in southern Italy and the islands of Sardinia and Sicily than T. mitsukurii, due to denser patterns from medium to high habitat suitability (Fig. 4B, C). The expectation for the future could be the potential spread of both species in these areas, even considering the current presence at lower latitudes of H. halys (Cianferoni et al. 2018; Tassini and Mifsud 2019).

The co-suitability map of parasitoids against H. halys at a global scale (Fig. 5A) shows that areas with good suitability for H. halys and scarce suitability for Trissolcus spp. are absent in the Asian region where, instead, there is a predominance of areas with good suitability for the host and for one or two parasitoids (marked with light blue and dark blue patterns, respectively), and areas with scarce suitability for H. halys but good suitability for one or two parasitoids (marked with light red and dark red patterns, respectively). The lack of yellow zones (i.e., good suitability for H. halys but scarce for Trissolcus spp.) in the Asian region means that the host distribution area completely overlaps the suitable areas of the parasitoids. It is possible that the two parasitoids and their host, all native to eastern Asia, show a more complex interspecific interaction. In native areas, the two parasitoids are able to exploit a larger host range (Zhang et al. 2017). Confirming this, in Australia, where H. halys was only intercepted and eradicated (Horwood et al. 2019; Hess et al. 2022), and shows low suitability (Fig. 2A), the presence of T. mitsukurii has been known since at least 1914 (Caron et al. 2021), suggesting the possible establishment of this parasitoid by exploiting other yet unknown hosts. Concerning T. japonicus, information on possible secondary hosts in the new areas in Europe and the United States is still scarce (Haye et al. 2020; Moraglio et al. 2020; Zapponi et al. 2021). Certainly, further investigations are needed to identify possible secondary hosts of both parasitoids, especially to assess their attractiveness in comparison with the coevolved host and, consequently, to track the new interactions that will be established in favorable areas for parasitoid establishment.

The areas where one or both parasitoids show suitability to coexist with H. halys (marked with light blue pattern and with dark blue pattern, respectively; Fig. 5A) reveal good chances of success in controlling the pest using them as BCAs. This means that the use of either one of the two parasitoids (light blue pattern, in the map) or both (dark blue pattern) can be effective in controlling H. halys depending on the individual habitat suitability (Fig. 2). Trissolcus japonicus shows a higher potential habitat suitability than T. mitsukurii at a global scale (Fig. 2B, C). Moreover, in Europe and Italy, where more structured surveys of the two parasitoids were performed (Zapponi et al. 2021), and according to the co-suitability map (Fig. 6), T. japonicus shows greater potential suitability than T. mitsukurii at lower latitudes. Therefore, these results further support the decision to field release T. japonicus as a part of a three-year national program initiated in northern Italy and extended throughout the peninsula (MATTM 2020). In light of the parasitoid co-suitability map (Fig. 6), the recent proposal to release T. mitsukurii against H. halys in the Northwest of France (Martel G., Hamidi, R., Thomas M., ANPN-RIPPOSTE, Région Nouvelle-Aquitaine, France, N° 15994520), in the area where this parasitoid was found in 2020 (Bout et al. 2021), appears to be the best solution considering its higher suitability, especially if compared with the lower suitability of T. japonicus in most of that country.

The predominance of the yellow and blue pattern (Fig. 5B) in the adventive range in Europe and North America suggests the possibility of successfully using one of the two Trissolcus species as a BCA, also from the perspective of a classical biological control program. To optimize the success of any BCA at the local scale, this model should be supplemented with an assessment of the local climatic parameters best suited to the parasitoid species. In particular, this is even more important in areas characterized by high habitat richness, in terms of geographic and climate variables (Cervellini et al. 2021), such as in the Iberian and Scandinavian peninsulas, eastern Europe, southern Greece (Fig. 6B), and, at Italian scale, in Sardinia and most of southern peninsular Italy (Fig. 7B).

Environmental factors play an important role in the performance of a BCA, because they can determine the success of any biological control program. Therefore, preliminary studies on the species habitat suitability can help to choose better performing BCAs to introduce into an ecosystem in a classical biological control program. Predicting suitable ecological niches for H. halys and its two parasitoids is a critical approach for crop management and biological control of this pest. The suitable core areas for T. japonicus and T. mitsukurii are quite similar to each other. The niches of these species are generally smaller than the area indicated by model-based predictions, and this is due to, among others, climatic variables that are not the only determinants of habitat suitability. Several natural and anthropogenic factors can also influence the potential habitat distribution of these species as well as their interaction. Therefore, the information provided by the model can help in the selection of the parasitoid to be used in relation to the suitability of the area, but it will still need to be substantially validated through field surveys, which can also further document the spatial intensities of the parasitoid species in potentially suitable areas.

Francesco Tortorici: Conceptualization, Data curation, Species identification, Investigation, Methodology, Software, Writing – Original draft. Pierluigi Bombi: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Writing – Original draft. Laura Loru: Funding acquisition, Resources Supervision, Writing – Review and Editing. Alberto Mele: Species identification, Investigation, Writing – Review and Editing. Silvia Teresa Moraglio: Investigation, Writing – Original draft. Davide Scaccini: Data curation, Species identification, Investigation, Writing – Review and Editing. Alberto Pozzebon: Funding acquisition, Resources Supervision, Writing – Review and Editing. Roberto Antonio Pantaleoni: Conceptualization, Funding acquisition, Resources Supervision, Writing – Original draft. Luciana Tavella: Conceptualization, Funding acquisition, Resources Supervision, Writing – Original draft. All authors have read and agreed to the published version of the manuscript.