Tritrophic analysis of the prospective biological control of brown marmorated stink bug, Halyomorpha halys, under extant weather and climate change

The highly destructive Asian brown marmorated stink bug (Halyomorpha halys, BMSB) invaded Europe, Caucasus region, and North and South America. Efforts to control it are ongoing in the Palearctic European-Mediterranean Basin and North America by introducing and redistributing two Asian stenophagous scelionid egg parasitoids (Trissolcus japonicus and T. mitsukurii) that are attacked by an adventive oligophagous pteromalid Asian hyperparasitoid (Acroclisoides sinicus). Large BMSB nymphs and adults may be parasitized by new associations of oligophagous tachinid flies and immature stages by egg parasitoids and predators. The terms stenophagous and oligophagous are commonly used to define narrow and wider ranges, respectively, of food eaten, but here they refer to the range of hosts attacked by adult female parasitoids. A holistic weather-driven physiologically based demographic model (PBDM) of the tritrophic interactions was developed to evaluate prospectively the impact of natural enemies on the biological control of BMSB under current and climate change weather. Our study focuses on the European-Mediterranean region, with the results for the USA, Mexico, and Central America reported as supplementary information. The PBDM analysis suggests that biotypes of the egg parasitoids T. japonicus and T. mitsukurii with high search capacity could suppress BMSB regionally, but the requisite levels of parasitism by these parasitoids for economic control are not observed in their native range nor in invaded areas. The model suggests that the action of T. japonicus is greater than that of T. mitsukurii, but that joint interactions of the two egg parasitoids would provide higher mortality of BMSB. Field data and model results suggest that the egg hyperparasitoid A. sinicus has a modest negative impact on the suppression of BMSB. Moreover, tachinid parasitoids of adults could have an important supplemental role in suppressing BMSB densities. Analysis suggests that new biotypes of egg parasitoids and species of tachinid parasitoids of large nymphs and adults be sought.


Introduction
The invasive polyphagous Asian brown marmorated stink bug Halyomorpha halys (Stål) (Heteroptera: Pentatomidae, BMSB) is native to eastern China, Japan, and Korea and has since been detected in most US states, Canada, Europe and Caucasus region, and recently in South America (see Rider et al. 2002;Leskey and Nielsen 2018). In invaded areas, new associations with native natural enemies have arisen (e.g., Dieckhoff et al. 2017;Joshi et al. 2019) but are insufficient for economic control (Abram et al. 2017). Further, chemicals widely used for control often provide unsatisfactory results (Leskey and Nielsen 2018). Classical biological control is seen as the long-term solution, and major efforts are ongoing in Europe and North America with the major focus being on the introduction and redistribution of two oligophagous Asian egg parasitoids: Trissolcus japonicus (Ashmead) and Trissolcus mitsukurii (Ashmead) (Hymenoptera: Scelionidae) that also attack eggs of other pentatomid species at low rates (Zhang et al. 2017;Hedstrom et al. 2017;Botch and Delfosse 2018;Sabbatini-Peverieri et al. 2018Maistrello 2019;Milnes and Beers 2019;Charles et al. 2019;Haye et al. 2020;Zapponi et al. 2021;Yonow et al. 2021;Kamiyama et al. 2021;Giovannini et al. 2022b). In native areas in Northern China, T. japonicus parasitism of H. halys eggs is 50-80% during July with an average of ~ 60% and comprised 77.2% field egg parasitoid collections (Zhang et al. 2017). Adventive populations of one or both egg parasitoids are present in Europe, USA and Canada (Talamas et al. 2015;Abram et al. 2017;Sabbatini-Peverieri et al. 2018;Stahl et al. 2019). In northern Italy, adventive populations of T. mitsukurii occur (Benvenuto et al. 2020;Zapponi et al. 2021), and T. japonicus is being reared and released in the area (MiPAAF 2020). Trissolcus mitsukurii in northern Italy has a higher impact compared to new associations of native egg parasitoids that attack BMSB (e.g., Anastatus bifasciatus (Geoffroy) (Hymenoptera: Eupelmidae), Trissolcus basalis (Wollaston), and Trissolcus kozlovi Rjachovskij (Hymenoptera: Scelionidae)) Scaccini et al. 2020). However, Abram et al. (2020) in review of the biological control efforts against four stink bugs cautioned about overestimating the impact of high levels of egg parasitism on the biological control of stink bugs.
In invaded areas, species of tachinid flies are known to form new associations with invasive stink bugs (see Abram et al. 2020). In North America, the tachinid Trichopoda pennipes (Fabricius) has high attack rates on squash bug (Anasa tristis (DeGeer)) and has been reported attacking BMSB nymphs and adults in the eastern USA, with overall parasitism rates of 2.38% to 4.5% during spring (Joshi et al. 2019). Trichopoda pennipes is recorded from Europe but has not been recovered from BMSB (Colazza et al. 1996). In our study, we assume the action of a generic tachinid species (i.e., Trichopoda spp.).
The focus of our study is to evaluate prospectively the potential impact of the natural enemies depicted in Fig. 1 on the biological control of BMSB across the Palearctic and Nearctic regions under current and climate change weather using holistic physiologically based demographic models (PBDMs) in a geographic information system (GIS) context given low background predation and parasitism of BMSB eggs and nymphs by new associations of native natural enemies. For heuristic purposes, we assume all the natural enemies in Fig. 1 are sympatric across the geographic range of our analyses.

Methods
Density-dependent weather-driven PBDMs of the species in the BMSB system enable analyses of their interacting dynamics at local and regional scales including the effects of climate change. PBDMs capture the daily weather-driven biology of the species making their predictions independent of time and place, and of detection records. PBDMs fall under the ambit of time varying life tables (Gutierrez 1996). Biological and ecological data available in the literature and developed by coauthors were used to formulate tritrophic age structured PBDMs of the four species each operating on its own physiological time scale (Gutierrez and Baumgärtner 1984; Appendix 1; see Gutierrez 1996). The biology of each species is captured by biodemographic functions described in the text and summarized in Table 1. Abbreviations Tj for T. japonicus, Tm for T. mitsukurii, As for the hyperparasitoid A. sinicus, and Tp for tachinid parasitoids are used in the text. The effects on BMSB of new associations with extant natural enemies (e.g., parasitoids and predators) and tachinid parasitoids are included in the system model as simple density-dependent mortality models.
Annual summaries of the simulation results are mapped using GRASS GIS, and the data are analyzed using multiple linear regressions (MLR). The analyses focus on the European Palearctic-Mediterranean Basin, with applications of the models to the continental USA, Mexico, and Central America summarized as supplementary information.

Diapause and winter survival
The biology and ecology of BMSB were reviewed by Leskey and Nielsen (2018). The hemimetabolous BMSB has egg, nymphal, and preovigenic and ovigenic adult life stages. BMSB overwinters as adults that enter diapause during fall with declining temperatures and photoperiod below 12.7 h daylength (dl) (Nielsen et al. 2016;McDougall et al. 2021). The proportion of adults entering diapause ( In (t) ) starting at time t when dl F = 12.7 dl is computed using a linear model (Eq. 1).
Diapausing adults are chill intolerant and seek refuge in sheltered places including human made structures (Cira et al. 2016(Cira et al. , 2018Ciancio et al. 2021). In West Virginia and Maryland, Lee et al. (2014) found overwintering diapausing BMSB adults in dry crevices in dead, standing trees with thick bark, particularly oak (Quercus spp.) and locust (Robinia spp.), but none from downed trees or leaf litter. Lee et al. (2014) opined that BMSB uses chemical or tactile cues to form overwintering aggregations. In cold north temperate areas, adults unable to find suitable shelter have extremely low survival rates, whereas survival in human structures may be quite high (Cira et al. 2016(Cira et al. , 2018Ciancio et al. 2021), enabling BMSB to persist in areas with extreme freezing minimum temperatures. The proportion of the overwintering population in human structures is unknown. Vermunt et al. (2012) developed a large data set of under-bark temperatures of ash trees (Fraxinus spp.) in woodlot environments at six different Ontario, Canada, locations. Rough estimates of minimum temperatures under bark were developed using linear regression of their data; specifically, T underbark = 0.414T ambient + 9.213; R 2 = 0.60 . The mortality rate per day of unsheltered diapause adults on under bark minimum temperature was estimated as the dashed line (gray symbols •), while the mortality rate of BMSB at various mean temperatures during the warmer season is characterized by the dotted line (black symbol •) (Fig. 2).
As temperatures warm in spring above 12 °C and photoperiod increases above dl S =13.5 dl, BMSB adults begin to emerge from diapause, feed, and enter the preovigenic stage, the duration of which is ~ 1.3 times longer than the preovigenic period of summer adults (Nielsen et al. 2017;McDougall et al. 2021). The proportion of adults leaving diapause at time t is a linear function of increasing dl ( Out (t))

Developmental biology
The developmental rate of BMSB life stages ( BMSB R = 1/ days) as a function of mean temperatures (T) has been well studied (Nielsen et al. 2008;Baek et al. 2017; field studies of Costi et al. 2017). The data for the nymphal stages are illustrated in Fig. 3a. A simple function was fit to the data with an estimated lower thermal threshold of 12.1 °C and an upper threshold of about 33.5 °C where the developmental rate function begins to decline sharply to zero at ~ 35 °C (see also Briére et al. 1999). Table 1 summarizes the developmental times of the different stages in degree days (dd > 12.11 °C) in the linear range of favorable temperatures.

Oviposition
Adult BMSB females are polyphagous and lay eggs in an age-temperature-dependent manner. Because BMSB is polyphagous, oviposition sites are not considered limiting. The profile of average per capita egg masses per female age in weeks at 26 °C after completion of the preoviposition period is depicted in Fig. 3b. The oviposition profile was converted to average eggs per day (i.e., f(i); Bieri et al. 1983) at age i = 1, …, A K with A K equal the mean 80-day life expectancy of adults. The concave effect of temperature on reproduction is introduced as a scalar (0 ≤ ϕ(T) ≤ 1 in the interval 18.25 °C < T < 34 °C; Fig. 3c) as ϕ(T)f(i)) ( Fig. 3d). Total oviposition (E(t)) by all females at time t is the sum of the product of the age-specific per capita fecundity and the number of adults ( A N (i,t)) of age i corrected for the sex ratio (sr assumed 0.5) (Eq. 4). (3) BMSB R(T) = 1∕days(T) = 0.0018(T − 12.11) (1 + 4.8 T−33.5 ) , R 2 = 0.96 Biology of the egg parasitoids assumed to be the same as for field sheltered BMSB adults. Parasitoid adults emerge from overwintering in relative synchrony with the beginning of reproduction by spring BMSB adults emerging from diapause. Further, the predation rate of parasitized and unparasitized eggs by native species during the season is assumed the same.
The per capita age-specific reproductive profiles of adult females of T. japonicus and T. mitsukurii are illustrated in Fig. 4a (Sabbatini-Peverieri et al. 2020); the effect of temperature on oviposition in Fig. 4c,f; the rates of egg to adult development in Fig. 4b,e; and the heavily female biased sex ratios on temperature in Fig. 4d,g. The lower thermal threshold of T. japonicus is 11.16 °C and that of T. mitsukurii is 12.25 °C, the upper thermal threshold of both is ~ 34 °C, and the oviposition profile of T. japonicus > T. mitsukurii, with total fecundity over a 10-day period being 91 and 82 eggs, respectively (Sabbatini-Peverieri et al. 2020). The fitted biodemographic functions are summarized in Table 1. Biodemographic functions for BMSB parasitoids T. japonicus (Tj (•)) and T. mitsukurii (Tm (○)). a Age-specific fecundity profiles at 25 °C for both species. For T. japonicus: b developmental rate of immature stages, c temperature-dependent oviposition, and d temper-ature-dependent sex ratios. Similarly for T. mitsukurii: e developmental rate of immature stages, f temperature-dependent oviposition, and g temperature-dependent sex ratios (data cf. Sabbatini-Peverieri et al. 2020) Primary parasitoid attack rates Salt (1935) outlined the steps in parasitoid attack biology as host finding (search), acceptance (preference), and host suitability (physiological suitability) that are components of parasitoid functional and numerical responses. Information on these factors for the parasitoids is incomplete, and hence we assume equal preference and suitability of BMSB eggs and explore the effects of search rates via simulation scenarios (see below).
The functional responses of T. japonicus and T. mitsukurii females are assumed similar with differences in age-specific demands for hosts and responses to temperature. Trissolcus japonicus is considered the most viable candidate for the biological control of BMSB (Zhang et al. 2017;Conti et al. 2021), and hence we use its biology as a template for T. mitsukurii. Ignoring the time subscript, the egg load (i.e., demand, D) of parasitoid females drives search for hosts, and for T. japonicus (left superscript Tj), the total age structured population demand ( Tj D) by all females is where Tj f(i) is the age (i)-dependent per capita oviposition profile (Fig. 4a), Tj(i) is the number of T. japonicus adults of age i, Tj sr(T) is the temperature (T)dependent sex ratio (Fig. 4d), and 0 < Tj ϕ(T) < 1 (Fig. 4c) is the normalized concave function that corrects oviposition demand for the effects of temperature. The combined demand for hosts by both species is Both parasitoids prefer to attack younger-age, unparasitized host eggs (E, i = 1, …,0.75 E K where E K is the maximum age of BMSB eggs and 0.75 E K approximates the maximum age attacked) (Qiu et al. 2007;Yang et al. 2018). Hence, the number of eggs available for attack by both parasitoids is Egg parasitoids rely heavily on the presence of semiochemicals to search the immense spatial gaps between egg masses in their environment (Conti and Colazza 2012). Search (parameter α with left superscript for species) is embedded in the related demand-driven ratio-dependent Gilbert-Fraser parasitoid form of the functional response model (Eq. 7i) and the related Gutierrez-Baumgärtner predator form of the models (Frazer and Gilbert 1976;see Gutierrez 1996). Specifically, the search parameter Tj α = 1 implies random search, Tj α < 1 is less than random, and Tj α > 1 is greater than random as one might expect from a highly efficient parasitoid using host chemical cues to increase search capacity. In field studies in the USA, BMSB egg masses were paired with those of three native pentatomids and exposed to adventive T. japonicus populations resulting in 67% parasitism of H. halys egg masses and 0.4-8% on the other species (Milnes and Beers 2019). This confirmed the oligophagous biology of T. japonicus, and the weaker suitability/preference of other hosts. Bertoldi et al. (2019) found a strong preference in T. japonicus for olfactory cues of its putative coevolved BMSB host, even when it emerged from other species (Boyle et al. 2020). In northern Italy, adventive populations of T. mitsukurii also display a strong preference for BMSB (Rondoni et al. 2022).
In the model, host preferences of both parasitoids for BMSB eggs are assumed the same (i.e., Tj φ = Tm φ = 1), but differences can be included explicitly in the model if differences are documented. The egg parasitoids compete, and super-and multiple parasitism occur (Giovannini et al. 2022a), and the Gilbert-Fraser parasitoid model (G-F) accommodates this biology. In the G-F model, the number of hosts attacked (Na) by the joint action of the two parasitoids is However, to partition the hosts attacked by each species, we first compute the number attacked by each species acting independently with search rates ( Tj α, Tm α) yielding Tj Na and Tm Na.
The proportion of multiple parasitism is Tj ∩ Tm . If there is no competitive advantage, the cases of multiple parasitism are divided proportionally as Tj Na = Na ⋅ Tj ∕( Tj + Tm ) and, Tm Na = Na ⋅ Tm ∕( Tj + Tm ) . However, if complete competitive advantage exists, the division will favor the dominant parasitoid (say T. japonicus) as Tj Na = Tj ⌢ H and Tm Na = ( Tm − Tj Tm )Ĥ . Partial dominance can also be accommodated as can difference in host preference and search rates.

Biology of the hyperparasitoid A. sinicus
The population dynamics of BMSB, T. japonicus, and T. mitsukurii are impacted by the obligate pteromalid hyperparasitoid A. sinicus that attacks pre-pupal and pupal stages of primary egg parasitoids ( Fig. 1) in the family Scelionidae (Giovannini et al. 2021 The figures for the demographic functions are like those of the egg parasitoids and hence are not illustrated (see Table 1). The lower thermal threshold for development is 12.65 °C, and the upper threshold is ~ 35 °C. The nonlinear developmental rate As R(T) at mean temperature T is At time t, the population demand for primary parasitoid hosts (Mele et al. 2021) is where i is age in days after a 2-day preoviposition period, As sr is the sex ratio, and the symmetrical concave scalar As (T) corrects for the effects of temperature on oviposition in the range 15-33.5 °C. In laboratory studies, the sex ratio of A. sinicus offspring from T. mitsukurii-parasitized eggs ranged from 62 to 91% females and 50-80% females from T. japonicus (Giovannini et al. 2021). We use average values in the model. The number of T. japonicus and T. mitsukurii parasitized BMSB eggs available as hosts for A. sinicus with preference values As Tj and As Tm is We use As Tj = As Tm , but the model anticipates differences in preference. The combined number of both species hyperparasitized ( As Na ) is computed using the G-F (8) As R = 0.0095(T − 12.65) 1 + 4.5 T−35 .
The number of T. japonicus parasitized eggs attacked by A. sinicus corrected for preference is: and the number of T. mitsukurii attacked is The hyperparasitism rates of the two egg parasitoids are Tj = Tj Na∕ As ⌢ H and Tm = Tm Na∕ As ⌢ H.

Parasitism by tachinid flies
Tachinid species attacking invasive pentatomids in invaded areas are known to form new associations as occurred in North America with T. pennipes attacking large BMSB nymphs and adults, albeit at low rates (Joshi et al. 2019).
Adult mortality due to tachinid parasitism was shown to be density-dependent (Liljesthröm and Bernstein 1990), and hence we model the low attack rates on BMSB adults using a weak density-dependent functional response model. The attack rate (μ adult-para , Eq. 12i) as a proportion is a function of host density with the rate increasing with physiological time above a temperature threshold (0 < Δt(T) = T-12.75 °C).

Predation of eggs and nymphs
Healthy and parasitized BMSB eggs and nymphs are attacked by native parasitoids and predators, and this mortality is also included as a density-dependent mortality rate (μ pred ) on prey density and physiological time (0 < Δt = T-12.75 °C) (Eq. 12ii). (10) Tm Na = As Na ⋅ As Absent sound data, immature BMSB and egg parasitoid life stages are assumed attacked equally.

Weather data
Daily weather from a reanalysis of weather observations (current climate) and from weather scenarios of climate change were used to run the mechanistic PBDM models. The weather data for current climate are from AgMERRA (Ruane et al. 2015), a global baseline forcing dataset for the period 1980-2010 from the Agricultural Model Intercomparison and Improvement Project (AgMIP, http:// www. agmip. org/). The data are a reanalysis of weather observations combined with observational datasets from in situ observation networks and satellites. The AgMERRA weather data (https:// data. giss. nasa. gov/ impac ts/ agmip cf/) have a ~ 25 km spatial resolution for each of the 17,791 lattice cells for the Euro-Mediterranean region (and 15,843 lattice cells for the USA, Mexico, and Central America). A subset of daily weather for the period January 1, 2000, to December 31, 2010, was used to characterize the current period.
Climate change weather data for the European-Mediterranean region for the period 2040-2050 were downscaled from coarser ~ 200 km resolution global climate model data to a ~ 30 km resolution using the PROTHEUS regional climate model (Artale et al. 2010). PROTHEUS is a coupled atmosphere-ocean regional model that allows simulation of local extremes of weather via the inclusion of a fine-scale representation of topography and the influence of the Mediterranean Sea (Artale et al. 2010). This enables increasing the spatial resolution and accommodating poikilotherm sensitivity to fine-scale climate and local topography effects. Specifically, we used the A1B regional climate change scenario that is toward the middle of the IPCC (2014) range of greenhouse gas (GHG) forcing scenarios (Giorgi and Bi 2005). The uncertainty associated with climate model predictions forced using the A1B scenario is low for the Mediterranean region relative to the rest of the globe (Gualdi et al. 2013). A subset of daily weather for the period January 1, 2040, to December 31, 2050, for 10,598 lattice cells was used for the Euro-Mediterranean region.
Daily climate change simulations for North and Central America weather are from the NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) dataset for the period 2055-2065 (Thrasher et al. 2012) at ~ 25 km resolution (https:// www. nccs. nasa. gov/ servi ces/ data-colle ctions/ land-based-produ cts/ nex-gddp) and were derived from global climate simulations of the Max Planck Institute Earth System Model low resolution (MPI-ESM-LR) model (Taylor et al. 2012) using a statistical downscaling technique (Thrasher et al. 2012). Global MPI-ESM-LR climate simulations were forced by the Representative Concentration Pathway 8.5 (RCP 8.5) scenario (Riahi et al. 2011) that corresponds to a range of warming similar to A1B (Rogelj et al. 2012). An evaluation of historical simulations of North American climate using continental metrics of bias relative to observed weather (Sheffield et al. 2013) showed that MPI-ESM-LR is the top ranked among the core set of 17 global climate models considered. Based on this downscaled climate change scenario, a daily weather dataset was developed for the period January 1, 2055, to December 31, 2065, for 20,355 lattice cells in North and Central America.

GIS mapping and marginal analysis
Using low nominal initial densities for each species, the PBDMs were run continuously for each of the lattice cells for specified time periods using daily lattice cell weather data. The daily age structured dynamics are predicted by the models for each lattice cell, but only yearly georeferenced summary data for all variables were written to year-specific text files. At the end of multi-year runs, the mean, standard deviation, and coefficient of variation were computed for each lattice cell omitting the first year when the model is assumed equilibrating. No calibrations of the models were made to fit common-wisdom notions of the geographic distribution and abundance of a species. The open-source GIS GRASS (Neteler et al. 2012) originally developed by the United States Army Corps of Engineers was used for geospatial data management, analysis, and mapping using bicubic spline interpolation on a 3-km raster grid to match the resolution of the underlying digital elevation model. The interpolated PBDM raster maps were overlaid on base map layers made using GRASS GIS with Natural Earth free vector (e.g., administrative boundaries) and raster (e.g., shaded relief) map data available in the public domain (https:// www. natur alear thdata. com/). The digital elevation model used in GIS computations (e.g., for restricting mapping of PBDM raster maps below a specified elevation) is the public domain NOAA Global Land One-km Base Elevation (GLOBE; https:// www. ngdc. noaa. gov/ mgg/ topo/ globe. html). GRASS is supported and further developed by the GRASS Development Team (2022, see http:// grass. osgeo. org). GRASS versions 6.4.4 and 8.2 were used on Windows and Apple computers, respectively.
Multiple linear regression (MLR) models (Eq. 13) with unknown error term ε were used to summarize the simulation data and to explore the impact of weather and computed variables (x i ) and the J interactions of x i (i.e., z j ) on dependent variables such as total annual number of BMSB host stages and of parasitized hosts.
To estimate the marginal effects of x i , the partial derivative of Eq. 13 with respects to x i is computed using the averages of remaining independent variables nymphs x i , i = 1, ..., I .

Results
To illustrate the richness of model output, the daily dynamics of BMSB, T. japonicus, T. mitsukurii and A. sinicus life stages are illustrated for the period January 1, 2005, to BMSB nymphs = f (x 1 , .....,  Daily maximum and minimum temperature (°C), day length (hrs), and daily mm rainfall (mm) are illustrated in Fig. 5a, the dynamics of the egg, nymph, preovigenic and reproductive summer adult stages in Fig. 5b, and the transition dynamics of diapause adults in fall to spring period when surviving preovigenic adults emerge from overwintering sites and become reproductive adults are illustrates in Fig. 5c. Panel Fig. 5d,e,f shows the dynamics of immature and adult stages of T. japonicus, T. mitsukurii and A. sinicus, respectively. For comparative purposes, 1,000-density levels for an arbitrary standard unit of the searched environment are indicated by horizontal dashed lines.
The annual cycle of BMSB near Padua has two egg peaks, one from adults emerging from diapause and the second from summer generation adults. Depending on the physiological time length of the season, a further generation could develop in warmer areas. In fall, maturing nymphs enter the preovigenic stage and become diapausing adults when photoperiod decreases below 12.7 h dl, and break diapause as dl increases above 13.5 h (Fig. 5c).

Effects of T. japonicus and T. mitsukurii search rates
The effects on BMSB nymphal densities of different search capacities of the two Trissolcus egg parasitoids was investigated using search rate scenarios Tj α = Tm α [0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2] for 2005-2010 weather near Padua given the background action of native biotic mortality factors (eqs. 12i, 12ii). The results yield decreasing nymph densities with increasing α and suggest that Tj α = Tm α > 0.5 would result in excellent suppression (Fig. 6). However, field parasitism rates in Asia, Europe and North America are usually lower in the range 50-80% (Zhang et al. 2017;Milnes and Beers 2019), suggesting field search rates of Tj α = Tm α in the range 0.2-0.3.

Estimating the prospective effects of T. japonicus, T. mitsukurii, A. sinicus and Trichopoda spp.
To examine prospectively the impact of the four parasitoids, 25 combinations of the dynamics of the four species with two egg parasitoid search rates Tj α = Tm α = 0.2 or 0.4 with As α = 0.1 were simulated using near Padua weather for the period 2005-2010. Note A. sinicus must always be run in combination with either or both egg parasitoids, while Trichopoda spp. may occur alone in the simulation, and the case of no natural enemies is the null scenario. Background predation was included in all runs (Eq. 12ii). In the MLR analysis, only 2006-2010 simulation data were used, and only presence or absence values (i.e., 1 or 0) for T. japonicus (Tj + ), T. mitsukurii (Tm + ), A. sinicus (As + ) and Trichopoda spp.(Tp + ) indicated by right superscript symbol + , search rates ( Tj α = Tm α) and season length (dd) were used as independent variables, while the dependent variable in the different analyses was average annual total number of each life stage. Only significant independent variables (p ≤ 0.05) were retained in the MLR model. Simple applications of marginal analysis are used to estimate the impact of the presenceabsence of each natural enemy.

Total eggs
The number of BMSB eggs surviving (Eq. 14i) is decreased by the action of the egg parasitoid (Tj + , Tm + ) and the tachinid Tp + . The large negative effect of increasing the search parameter of the egg parasitoids was expected as illustrated in Fig. 6, while the negative effect of a longer season length reflects the longer period for the action of the natural enemies. The regression coefficient for the hyperparasitoid A. sinicus is positive indicating it increases BMSB egg levels by reducing egg parasitoid density. Surprisingly, despite the low attack rate, the action of Tp + is about the same as the combined action of the two egg parasitoids. Removing season length in this (and all subsequent) MLRs yields the following relationships.
Using simulation data when both parasitoids are present in the system, the ratio of total T. mitsukurii parasitzed egg to those parasitized by T. japonicus is ~ 0.34 (R 2 = 0.45).

Total nymphs
The impact of the three primary parasitoids on total nymphs is roughly the same, with the action of the hyperparasitoid A. sinicus being positive but not significant and hence not retained in the model.

Total adults
The effects of the egg parasitoids on total adults are similar, with the impact of the tachinid fly being about the same as that of the two egg parasitoids combined. The action of A. sinicus is positive and on average reduces the combined action of T. japonicus and T. mitsukurii ~ 22%.
The total number of diapausing adults per year regessed on the total number of adults per year is y = 0.59x, R 2 = 0.81. The impact of Tp + acting alone on adult density (see Eq. 12i) is estimated using the data for Tp + = [0, 1] (i.e., absent or present).
Using average Tp + = 0.5 yields an average 19.9% seasonlong adult mortality and an average ~ 0.25% mortality rate per day during an 80-day BMSB adult lifespan. If the developmental time of Trichopoda spp. larvae is 8-10 days, then a point estimate of parasitism would be ~ 2.27% = 0.25% × 9, which is at the lower end of field estimates of 2.38% to 4.5% parasitism rates observed by Joshi et al. (2019). We estimate roughly the average impact of Tp + in our model as follows: suppose a BMSB female lays an average of ~ 210 eggs, hence the product of sex ratio x correction for death at midlife x average fecundity x the regression coefficient from Eq. 14v x average Tp + (i.e., 0.5 × 0.5 × 210 × 335.5 × 0.5) suggests that on average 8,794 BMSB eggs would not be produced during the season due to the action of Trichopoda spp. This value is a similar magnitude to the average effects of Tp + on the reduction of total annual BMSB eggs when all natural enemies are interacting (Eq. 14ii).

Natural enemy interactions
We first explore interspecific competition of T. mitsukurii and Trichopoda spp. on the number of T. japonicus parasitized eggs ( E Tj) (Eq. 15i).
The negative effects of Tm + and Tp + on Tj are similar, and though not significant As + hyperparasitism increased T. japonicus parasitism by decreasing competition from T. mitsukurii. A similar relationship occurs when T. mitsukurii parasitized eggs ( E Tm) is the dependent variable with a greater effect of T. japonicus competition and lesser effects of Tp + and As + (eqns. 15i vs 15ii). The low simulated hyperparasitism by As + reflects the low search rate As α = 0.1 suggested by field data (Mele et al. 2022). In the model, the average total season hyperparasitism rate on T. japonicus is ~ 23% and ~ 26% for T. mitsukurii.

Regional analysis
The model was run across all lattice cells in the European-Mediterranean region for the 2000-2010 period. Absent non-native natural enemies, Fig. 7a shows the prospective distribution of average total annual BMSB nymphs given the action of native natural enemies, and Fig. 7b maps the associated coefficients of variation (CV) as a percentage. Note the geographic distribution of the CVs is much wider than that for average nymph density indicating the additional marginal areas where BMSB could occur but at low highly variable levels. The efficacy of the two egg parasitoids depends on weather and their ability to find hosts. At the regional level, the latter is examined by comparing the results for search parameters Tj α = Tm α = 0.2 in Fig. 7c with Tj α = Tm α = 0.6 in Fig. 7d. The CVs for Fig. 7d are mapped in Fig. 7e. Visually, the densities do not appear to vary much between the Tj α = Tm α = 0.2 and 0.6 scenarios ( Fig. 7c and d), with the maximum densities in both scenarios being approximately 30% that absent egg parasitoids and tachinid parasitoids (Fig. 7a). However, subtracting the lattice cell values in Fig. 7d from corresponding ones in Fig. 7c shows large decreases in central Spain, North Africa, Italy and elsewhere (Fig. 8). Compared to areas of Italy with high BMSB nymph densities absent egg parasitoids and tachinid parasitoids, the prospective reductions in BMSB densities in many areas at α = 0.6 is ~ 85%. The action of the egg parasitoids and tachinid parasitoids increases the distribution of high CVs across the region (Fig. 7b vs. e).
Using MLR, we examine the effects on BMSB of season length (degree days, dd) and annual rainfall (mm) and primary parasitoid search rates ( Tj α = Tm α = [0.2, 0.4, 0.6]) across the region (Eq. 16). On average, BMSB nymph densities increase with annual rainfall and season length, and as shown for site-specific analyses decrease with increasing egg parasitoid search efficiency (α).

Climate change
Using climate model weather for 2040-2050, nymph densities in the absence of egg parasitoids and tachinids increase across the region (Figs. 9a, b vs. 7a, b). Including all the natural enemies with high egg parasitoid search ( Tj α = Tm α = 0.6 and As α = 0.1), prospective nymph densities decrease ~ 80% or more across the region (Fig. 9c) with increased variability (i.e., CVs) in many areas (Fig. 9d).

Discussion
Using a simple stage-structured matrix models, Abram et al. (2020)  Prospective average geographic distribution of annual BMSB nymph densities and coefficients of variation (CV) as a percent for 2040-2050 climate change daily weather. a Average nymph densities absent natural enemies and b CV (%). c Average nymph densities including the action of all natural enemies with egg parasitoid search parameters Tj α = Tm α = 0.6 and. As α = 0.1, and d CV (%). Nymphal densities are truncated at 2754 in Fig. 9c (red symbol) bug (Nezara viridula De Geer, Hemiptera: Pentatomidae) population growth; (ii) that evidence the introduced egg parasitoid T. basalis was responsible for complete control of N. viridula was weak despite ~ 50-90% parasitism of egg masses; (iii) that the apparent impact of T. basalis on N. viridula in Hawaii, Indonesia and the mainland USA (California) was variable or very limited, and sometimes negligible compared to the much larger impact of egg predators such as ants.
Our heuristic study used biologically rich weather-driven physiologically based distributed maturation time demographic models (PBDMs) designed to evaluate the potential of the parasitoids of BMSB egg and adult stages for the control of BMSB (Fig. 1) at local and regional levels. PBDMs capture the biology of the interacting species making their predictions independent of time, place, and detection records, and fall under the ambit of time varying life tables (Gutierrez 1996). This bioeconomic modeling paradigm (Regev et al. 1998) has a long history in assessing the geographic distribution and relative population dynamics of invasive species (Gutierrez and Baumgärtner 1984;e.g., Gutierrez 1996;Gutierrez and Ponti 2013a;Ponti et al. 2021; see supplementary information).
PBDM models (i.e., system) capture the biology of the invasive polyphagous brown marmorated stink bug (H. halys, BMSB) as attacked by two exotic scelionid egg parasitoids (Trissolcus japonicus and T. mitsukurii) that are in turn attacked by a pteromalid hyperparasitoid (Acroclisoides sinicus). As of 2019, adventive populations of T. japonicus and T. mitsukurii were not widely distributed in Europe (Zapponi et al. 2021), but releases are occurring, and the two species are expected to become sympatric across BMSB invaded areas. Our analysis anticipates this outcome. The biological interactions of the two parasitoids were studied in the laboratory by Costi et al. (2022) and Giovannini et al. (2022a) who found that indirect competition between females of the two species led to more progeny produced by the first arriving female of either species and in most cases neither species recognized previously parasitized host eggs, resulting in a wastage of eggs due to multiple parasitism. Multiple and super-parasitism and search are components of the Gilbert-Fraser parasitoid functional response models used in our study.
Host finding by T. japonicus and T. mitsukurii is increased by semiochemicals emitted by BMSB (Bertoldi et al. 2019;Boyle et al. 2020;Malek et al. 2021;Rondoni et al. 2022). Field parasitism of eggs in BMSB native areas ranged from 60 to 80% (Zhang et al. 2017) that in our model suggested an egg parasitoid search rate of 0.25 < α < 0.35 (Fig. 6) estimated via simulation. The simulation model predicts eggparasitism alone would be insufficient to keep the BMSB population below economic levels (see Abram et al. 2020), but would help reduce levels beyond the action provided by autochthonous natural enemies (Costi et al. 2019;Moraglio et al. 2020;Scaccini et al. 2020;Zapponi et al. 2021). Further, local conditions may promote higher levels of control, as observed by the higher impact of parasitoids in fruit orchards close to semi-natural areas compared to more distant ones (Mele et al. 2022).
Hyperparasitoids may also use semiochemical cues in host finding (Cusumano et al. 2020) that could reduce the effectiveness of biological control agents. Theoretical models suggest that some hyperparasitoids may dampen extreme host-parasitoid population oscillations and contribute to the stability of multitrophic systems restraining the efficacy of biological control, but such situations may be rare and, in most cases, invading hyperparasitoids are posited more likely to disrupt stable host-parasitoid interactions (Holt and Hochberg 1998). Giovannini et al. (2021) reared A. sinicus on both T. japonicus and T. mitsukurii, and though not designed to test for host preference, ~ 2.7 times as many A. sinicus were produced under similar conditions on T. mitsukurii than on T. japonicus. Laboratory studies by Mele et al. (2021) found > 20:1 preference for T. mitsukurii compared to T. japonicus using a colony of A. sinicus derived from and maintained on T. mitsukurii parasitized BMSB eggs. However, growing evidence indicates that genetic variation in behavior and conditioning during the adult stage in insects can contribute to preference for the host on which the insect developed (Barron 2001). In parasitoid wasps, conditioning during "critical periods" during the adult stage may have a profound impact on host selection behavior (van Emden et al. 1996;Bjorksten and Hoffmann 1998;Barron 2001). Hence, the strong laboratory host preferences demonstrated in A. sinicus by Mele et al. (2021) may not be applicable to the field where hyperparasitism levels are low (Mele et al. 2022). In the model, we assumed equal preference for T. japonicus and T. mitsukurii, but the capacity to incorporate unequal preference is a component. Average season total hyperparasitism predicted by the model on T. japonicus is ~ 23% and ~ 26% for T. mitsukurii.
Also included in the model was the putative low parasitism by polyphagous tachinid flies (i.e., Trichopoda spp.) that parasitize large nymph and adult stages of many Hemiptera including BMSB. Observed field parasitism of BMSB large nymphs and adults by the tachinid Trichopoda spp. is < 5% (e.g., Joshi et al. 2019), but the impact may be far greater than merely the death of individual BMSB adults, as the costs in time and energy to produce adults is near maximum and the preclusion of their potential egg production can have large demographic consequences that PBDMs are designed to assess. The impact of other extant native parasitoids and predators attacking BMSB egg and nymphal stages is low e.g., Costi et al. 2019; see also Rot et al. 2021) but was not explored.
Overall, the simulation results using near Padua, Italy weather suggests that T. japonicus is more effective than T. mitsukurii, the joint control by both egg parasitoids is greater than either acting alone, and the hyperparasitoid does not greatly increase BMSB densities. The potential impact of adult parasitism by tachinids could be equal to that of both egg parasitoids.

A regional perspective
BMSB densities in many invaded areas are high despite attack by new associations of native egg parasitoids and predators. Under these conditions, the prospective geographic distribution and relative abundance of BMSB in the Palearctic European−Mediterranean Region for the 2000-2010 period is composed of two areas: the smaller area where viable stable populations can occur (Fig. 7a), and the larger area with high coefficients of variation where low but highly variable populations may develop (Fig. 7b). On a regional scale, large decreases in nymphal densities are predicted with increasing egg parasitoid search rates ( Tj α = Tm α = 0.2 → Tj α = Tm α = 0.6) with low rates of hyperparasitism, and the additional action of the tachinid parasitoid and native natural enemies (Eq. 16). Furthermore, the enhanced action of the natural enemies increases CVs across the region under current (Fig. 7b vs. e) and even more under 2040-2050 climate change weather (Fig. 9a,b vs c,d). Using data from Fig. 7d, a regression of total T. mitsukurii parasitized eggs (left superscript E) to T. japonicus parasitized eggs across years and lattice cells across the region ( E Tm = 4.779 + 0.294 E Tj, R 2 = 0.94) suggests that given the same search rate and A. sinicus preference, T. japonicus could attack ~threefold more BMSB eggs than T. mitsukurii. This result reflects the lower developmental threshold and higher fecundity of T. japonicus. Unfortunately, consistent high levels of egg parasitism required for control are not recorded in native or invaded areas.
Research on BMSB has focused on easily observable egg parasitoids and less on the action of tachinid parasitoids of adults and predation by native species, likely because the latter are more difficult to assess. However, Abram et al. (2020) identified three tachinid species present in different areas that formed new associations with the invasive N. viridula: T. pennipes in North America, T. giacomelli (Blanchard) in South America, and T. pilipes (F.) in the West Indies. Further, Joshi et al. (2019) underscored the underappreciated potential contribution of low parasitism rates by polyphagous tachinid parasitoids to BMSB control. Our model suggests that despite low adult parasitism rates, the impact of the tachinid fly would be approximately that of the combined action of the two egg parasitoids. From a bioeconomic perspective (Gutierrez and Regev 2005), tachinid parasitism of large nymph and adult BMSB results in a maximum per capita loss in time and energy investments to the BMSB population in addition to the loss of eggs not deposited. In contrast, the loss of an egg results in little loss in time and energy. Hence, a recommendation that emerges is that new biotypes and species that attack the adult stage (e.g., P. latifascia, Chen et al. 2020) should be sought in BMSB's native range. Because extant populations of T. japonicus and T. mitsukurii are adventive, they may be viewed as new biotype associations on invading BMSB populations; hence, given the large effect of increased search rate as shown here, additional biotypes of the egg parasitoids should also be sought. A classic example of the potential of new biotype introductions is the biological control of walnut aphid (Chromaphis juglandicola (Kaltenbach)) in California by an Iranian biotype of the parasitoid Trioxys pallidus Haliday after an introduced French biotype failed to control the pest (van den Bosch et al. 1979;see Lozier et al. 2008).

Climate change
From a biological point of view, climate change is simply other weather patterns that may alter the interactions, relative abundance, and geographic distribution of species. The prospective effects of climate change on BMSB nymph abundance given the action of all natural enemies with search parameters Tj α = Tm α= 0.6 and As α= 0.1 (i.e., good biological control) were estimated by comparing regional 2000-2010 simulation results with 2040-2050 results (i.e., Figs. 7d vs. 9c). With climate warming, the results show that the distribution and abundance of BMSB is expected to decrease as some areas of the Mediterranean basin become less favorable due to increased temperature directly, and indirectly due to longer seasons in favorable areas that would enhance the action of the introduced natural enemies.

Can we believe model results?
Like any ecological model, PBDMs are incomplete as they are not one to one descriptors of nature. Gutiérrez Illan et al.
(2022) used the species distribution MaxEnt algorithm to predict the distribution of BMSB under current and climate change weather in the USA, concluding that occurrence was most affected by winter precipitation and proximity to populated areas, that abundance was influenced strongly by evapotranspiration and solar photoperiod, and that environmental factors that mediated the geographical distribution of BMSB are different from those driving abundance. They concluded that linking models of establishment (occurrence) and population dynamics (abundance) offers a more effective way to forecast the spread and impact of BMSB and other invasive species than simply occurrence-based models. We heartily agree. Laboratory studies suggest that desiccation in hot and dry climates may limit BMSB's geographical range (Scaccini et al. 2019;Khadka et al. 2020;Fisher et al. 2021;Grisafi et al. 2021;Stahl et al. 2021) supporting the Gutiérrez Illan et al. (2022) MaxEnt finding that evapotranspiration was a key factor limiting BMSB (see supplemental materials for applications).
A well-parameterized PBDM enables prospective examination of the biological interactions of multiple species as driven by weather at the local and regional levels not readily observable in the field. This point was made by Toko et al. (2019) who opined in review of the highly successful biological control program against cassava mealybug (Phenococcus manihoti Matile-Ferrero, Hemiptera: Pseudococcidae) in sub-Saharan Africa that the system model (Gutierrez et al. 1993) demonstrated that the encyrtid parasitoid Apanogyrus lopezi Howard reduced mealybug populations by 90%, while indigenous lady beetles cut occasional peak mealybug populations by 20%, a result that no field study could have demonstrated. A well-parameterized PBDM was used to analyze the biological control of the spotted alfalfa aphid (Therioaphis maculata Buckton, Hemiptera: Aphididae) across the disparate ecological zones of California by the combined time and place action of parasitoids, disease, and coccinellid predators (Gutierrez and Ponti 2013a). Tritrophic PBDM analyses of the biological control of invasive pest species, include vine mealybug (Planococcus ficus (Signoret); Gutierrez et al. 2008), glassy-winged sharpshooter (Homalodisca vitripennis (Germar); Gutierrez et al. 2011), Asian citrus psyllid (Diaphorina citri Kuwayama; Gutierrez and Ponti 2013b), coffee berry borer (Hypothenemus hampei (Ferrari); Cure et al. 2020) and others. In each case, the underlying bases for the success or failure of the biological control effort were informed by the models.
The tritrophic BMSB system model objectively captured the known global change biology and identified as key elements the need for high field search rates and host preference in the egg parasitoids T. japonicus and T. mitsukurii. It further elucidated the putative action of tachinid parasitoids of BMSB adults as an important component for successful region wide biological control of BMSB. High egg parasitism rates contribute to reducing the impact of BMSB regionally, but alone would appear to be insufficient to control its population to non-pest levels (see Abram et al. 2020). This finding has parallels in the PBDM analysis of the ongoing biological control of the invasive yellow starthistle (Centaurea solstitialis L., YST) in Western North America where partial success was achieved by the introduction of several capitula/seed-feeding insects (Gutierrez et al. 2017). A critical element identified for the biological control of YST was the need for natural enemies that kill mature whole plants and/or greatly reduce seed production by attacked survivors (e.g., the rosette weevil, Ceratapion basicorne (Illiger)). The area invaded by BMSB is increasing, and like YST we must await the system reaching dynamical equilibria before we know where along the spectrum of biological control scenarios it settles and why. In the interim, the PBDMs for BMSB can be easily updated as new information becomes available (e.g., evapotranspiration effects, cf. Grisafi et al. 2021) and can be used to evaluate prospectively new natural enemies for introduction and for providing insights and direction for moving efficiently toward the desirable goal of regional biological control of BMSB.

Author contributions
All authors reviewed and approved the manuscript.

Distributed maturation time population dynamics models
Physiologically based demographic models (PBDMs) are time-varying life tables, the theoretical basis of which was reviewed by Gutierrez et al. (1994) and Mills and Gutierrez (1999). It is commonly supposed that PBDMs have large numbers of parameters making them difficult to develop, but as demonstrated here and in prior studies, this is not the case, and further, the same underlying model(s) can be used across species and trophic levels.
Only an overview of the time-invariant distributed-maturation time demographic model used in our study (Manetsch 1976) is presented noting that other demographic models with distributed maturation time could also be used ( Vansickle 1977;Di Cola et al. 1999;e.g., Buffoni and Pasquali 2007). Species have different life stages (left superscript s), and the same discrete dynamics model is used for all of them. The model for the ith age class of a life stage s with i = 1, 2,…, s k age classes (Eq. 17; Gutierrez 1996;Severini et al. 2005) can be viewed as composed of s k dynamics equations for each life stage. The forcing variable is temperature (T), with time (t) being a day (d) that from the perspective of ectotherm species is of variable length in physiological time units (i.e., s Δ x (T(t)) in degree days (dd)), or proportional develoment ( s r(T(t))) that may differ for each stage (and species) having different mean developmental times s Δ and thresholds. The state variable s N i (t) is the density of the ith age class, and s i (t) is the proportional age-specific net loss rate due to temperature, net immigration, and other factors during s Δ x (T(t)) (Gutierrez 1996). Following the notation of Di Cola et al. (1999, page 523), the ith age class of stage s was modeled as follows: The total density in life stage s is s N(t) = ∑ k i=1 s N i (t). Ignoring stage notation, new eggs enter the first age class of the egg stage for say BMSB (i = 1), flow occurs between age classes and between stages (egg to nymph to preova and mature adults), and surviving adults exit as deaths at maximum age (i = A k). Absent mortality, the theoretical distribution of cohort developmental times is estimated by Erlang parameter k = Δ 2 ∕var , where var is the variance of developmental time Δ. In our study, we used k = 25 for the egg and nymphal stages, while for the BMSB adult stage A k = 80 is equal the average longevity in days at 25 °C. One can use proportional development or physiological time units (dd, degree days) more commonly interpreted by biologists. For example, if the average BMSB nymphal stage is 446 dd, the standard deviation is 89.2 dd for k = 25. Furthermore, developmental times may vary with nutrition and other factors (Gutierrez 1996) and given appropriate data, this aspect can be easily accommodated using the time varying form of the model (Vansickle 1977). Lastly, because of non-linearities and time varying nature, the model can only be evaluated numerically (Wang and Gutierrez 1980). The numerical solution for Eq. 18 can be found in Abkin and Wolf (1976), as implemented here in Gutierrez (1996, pages 157-159), or as Pascal code from APG. from the NEX-GDDP dataset, prepared by the Climate Analytics Group and NASA Ames Research Center using the NASA Earth Exchange, and distributed by the NASA Center for Climate Simulation (NCCS). The data for the climate weather scenario for the Euro-Mediterranean region were developed and provided by the Laboratorio Modellistica Climatica e Impatti at ENEA, Rome, Italy. Model development was by CasasGlobal.org, Kensington, CA, USA. The work of CREA was supported by the Italian Ministry of Agriculture Food and Forestry Policies (grant projects "Salvaolivi" DM-0033437-21/12/2017 and "Proteggo 1.4" DISR-05-0001837-04/01/2022). A.P., D.S. and A.M. were supported by funding provided by Regione Veneto-Unità Organizzativa Fitosanitario. The study was supported by project TEBAKA (project (17) d s N i dt = s k ⋅ s Δ x s Δ s N i−1 (t) − s N i (t) − s i (t) s N i (t).
(18) d dt s Δ s n i (t) s k = s n i−1 (t) − s n i (t) − s i (t) s n i (t) s Δ s k .
ID: ARS01_00815) co-funded by the European Union-ERDF and ESF, "Pon Ricerca e Innovazione 2014-2020." Funding Open access funding provided by Ente per le Nuove Tecnologie, l'Energia e l'Ambiente within the CRUI-CARE Agreement. CASAS NGO funded model development. The work of CREA was supported by the Italian Ministry of Agriculture Food and Forestry Policies (grant projects "Salvaolivi" DM-0033437-21/12/2017 and "Proteggo 1.4" DISR-05-0001837-04/01/2022. A.P., D.S. and A.M. were supported by funding provided by Regione Veneto-Unità Organizzativa Fitosanitario. The study was supported by project TEBAKA (project ID: ARS01_00815) co-funded by the European Union-ERDF and ESF, "Pon Ricerca e Innovazione 2014-2020." Data availability All biological data used in the analysis are publicly available or were sourced from tables, text or estimated from figures in the cited literature. All data from the different sources are shown in the fitted biodemographic functions plotted in Figs. 2, 3, and 4, with references cited in the figure legends and included in the reference list (Nielsen et al. 2008;Cira et al. 2016Cira et al. , 2018Baek et al. 2017;Costi et al. 2017;Sabbatini-Peverieri et al. 2020;Ciancio et al. 2021). Base geodata layers used to generate maps are available in the public domain as Natural Earth vector and raster data (https:// www. natur alear thdata. com/) and GLOBE digital elevation model data (NOAA Global Land One-km Base Elevation, https:// www. ngdc. noaa. gov/ mgg/ topo/ globe. html).

Conflict of interest
The authors have no relevant financial or nonfinancial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.