Key message

  • We developed a reaction–diffusion model to predict P. japonica's spatio-temporal dynamics.

  • The parameterisation is based on data from the Phytosanitary Service of Lombardy from 2015 to 2021.

  • The rate of expansion from the initial point of invasion varies along different directions.

  • We used a linear model to describe the influence of land use on the rate of expansion of the pest.

  • We described the spatio-temporal dynamics patterns and predicted future trends of pest expansion.

Introduction

The Japanese beetle (Popillia japonica Newman) is a highly polyphagous pest of the Scarabaeidae family. The beetle is native to Asia and, outside its area of origin, the species is present in North America and continental Europe. In North America, the species is a major pest of turf and crops (Potter and Held 2002). The first report of the Japanese beetle in continental Europe dates to July 2014 in northwestern Italy, in the Ticino River natural park, located between the Piedmont and Lombardy regions (Pavesi 2014). The species most likely established in that area a few years before 2014, and at the time of the first detection the population was actively expanding. Currently, the species has been able to invade more than 16,000 km2 of territory, covering vast areas in Lombardy and Piedmont and reaching some municipalities in northwest Emilia Romagna and in Valle d’Aosta (EPPO 2022b; Gotta et al. 2023; Europhyt Outbreak No. 574, Update 07/2022–12-05). Since June 2017, the pest has been present in Switzerland in an area located at the border with Italy, where it is currently under a containment strategy (EPPO 2021). Most recently (26th of July, 2023), the pest has established for the first time north of the Alps, near the airport of Zurich (https://www.Popillia.eu/blog/first-outbreak-of-japanese-beetle-north-of-the-alps). In Europe, the species has been included in the A2 list of the European Plant Protection Organization as a quarantine pest (EPPO 2022a) and, more recently, it has been added to the list of priority pests (European Commission, 2019).

To develop approaches for the rational containment of an invading population, it is important to understand the processes of population growth and spread as they are fundamental in describing and predicting the spatio-temporal dynamics of pest expansion in a naïve environment. The understanding of these population processes leverages models that summarise the available data and simulate the growth and spread of the pest during an ongoing invasion phase. Despite their importance, current models describing spatio-temporal dynamics of P. japonica are quite limited.

Population growth represents a major driver influencing the spread of an invasive species. The increase in the number of individuals over time enables the attainment of the propagule pressure necessary to colonize new patches where environmental conditions are suitable for local establishment (Nathan et al. 2003; Liebhold and Tobin 2008). Data are available on the role of environmental variables on the main processes linked to P. japonica population growth, i.e. individual fertility, development, and survival (or mortality). For instance, temperature (Ludwig 1928; 1930; 1932; Payne 1928; Fleming 1972), air relative humidity (Payne 1928; Ludwig 1936; 1945; Ludwig and Landsman 1937; Bellucci 1939), soil humidity (Fleming 1972; Regniere et al. 1981; Allsopp et al. 1992), and soil texture (Regniere et al. 1981; Allsopp et al. 1992) play major roles in influencing the life-history of the species. Some studies have also investigated the ecological niche of P. japonica (Hirzel and Le Lay 2008; Simonetto et al. 2022; Zhu et al. 2023), and identified the range of suitable conditions for the species which might also favour population growth (Hirzel and Le Lay 2008). Despite the availability of detailed information on pest biology and ecology, a quantitative description of the pattern of temporal population dynamics with an estimation of the intrinsic population growth rate are scarce in the literature. However, the analysis of the ongoing spread of P. japonica in Lombardy allowed the identification of a clear logistic growth pattern (Gotta et al. 2023). Results obtained by the GESPO project (funded by Lombardy Region) have led to the estimation of a function describing the population growth of P. japonica in Lombardy. These results are available at the project’s website (https://sites.google.com/unibs.it/gespo-project-en/outcomes/estimating-population-abundance-and-growth).

The sum of all individual’s dispersal events combined with population growth and successful establishment in novel areas result in the progression of the front of an invasion, i.e. the population spread (Liebhold and Tobin 2008; Lewis et al. 2016; Nathan et al. 2012). In relation to individual dispersal, we can further distinguish between short-distance dispersal (SDD), i.e. the movement of individuals, commonly over a short range within their environment, and long-distance dispersal (LDD), i.e. the movement of individuals over a wide range, commonly aided by humans (through the movement of infested materials, hitchhiking events, etc.) (Liebhold and Tobin 2008; Gilioli et al. 2013; Tobin and Robinet 2022). The diffusive process of P. japonica observed in northern Italy seems to be primarily triggered by individual SDD (Borgogno-Mondino et al. 2022), resulting in a continuous spread at the population level. Mark–recapture studies reported rather variable average individual dispersals, ranging between 50 m/day (Lacey et al. 1995) and 2.4 km/day (Lessio et al. 2021). Studies from the US reported a maximum daily flight distance covered by the pest ranging between 400 m (Holmes and Barrett 1997) and 700 m (Hamilton 2003). Several variables influence P. japonica individuals’ dispersal: time of the day (Paoli et al. 2023), wind speed (Lacey et al. 1994), air humidity (Fleming 1972), cloud coverage (Lacey et al. 1994), land management and pesticide treatments (Bohlen and Barrett 1990; Lacey et al. 1995; Holmes and Barrett 1997). Additionally, the gregarious behaviour of the adults (Kowles and Switzer 2012) plays an important role in influencing individual dispersal. This is evident as the insect tends to follow hormonal cues left by a pioneer female or volatile molecules emitted by plants suffering foliar damage (Potter and Held 2002) or preferred as hosts (Fleming 1972). Regarding wind and its contribution to the dispersal of the pest, Fox (1927) suggested that it plays a fundamental role in the westward expansion of the pest in the US. On the other hand, Borgogno-Mondino et al. (2022) indicated that wind plays a marginal role in the population expansion in the Po Valley, primarily because it is usually low or absent during the peak flight activity hours.

The successful establishment of dispersing propagules contributes to the spread of a species (Nathan et al. 2012). The population spread can be measured as a rate of expansion, either by the distance reached by the pest from the initial invasion point (Lutscher 2019), or by the area occupied by the species (Shigesada et al. 1995). Most of the studies pertaining to the spread of P. japonica come from the US. Hadley and Smith (1923) quantified the growth of the infested area from the first year of detection (period from 1916 to 1922), while Fox (1927) reported distances reached by the pest in 1926 from the initial point of invasion, emphasising the varying expansion rates along different directions. More recently, Allsopp (1996) and Lewis et al. (2016) investigated the changes in the rate of radial expansion, highlighting an exponential increase in the initial phase, followed by a linear increase in the following two decades. For the Italian case, Borgogno-Mondino et al. (2022) proposed an iterative model based on a geostatistical approach.

In our work we relied on a mechanistic reaction–diffusion model to integrate population growth and diffusion to describe and predict the spatio-temporal dynamics of P. japonica populations (Roques et al. 2008; Gilioli et al. 2013). Biologically meaningful parameters were estimated using the monitoring data of the adult stage of P. japonica collected by the Phytosanitary Service of Lombardy Region between 2015 and 2021. The reaction component (i.e. population growth) was described using a logistic Beverton–Holt function (Beverton and Holt 1957; Liang et al. 2017). The diffusive component was modelled by a Gaussian kernel. An integrodifference equation model was derived to describe the rate of expansion, i.e. the speed at which the leading edge of a population wave moves over time, indicated as the asymptotic spreading speed in Lutscher 2019. To investigate the potential influence of the landscape on the rate of expansion of P. japonica (Fox 1927), the diffusive component was estimated on 13 different trajectories originating from the area of first detection of the pest in Italy and a linear regression model was applied to interpret the effect of land cover and use on the spread rate of P. japonica.

Materials and methods

Study area

The study area (Fig. 1) corresponds to the area infested by P. japonica in the Lombardy Region (northern Italy) from 2015 to 2021. The southern portion of the area comprises a part of the Po Valley, an alluvial plain formed by the Po River, characterised by arable lands, rice paddies, broadleaf forests, and natural and semi-natural vegetation. The Milan metropolitan area extends in the central part of the infested region and is characterised by a highly anthropized and fragmented landscape. Further to the north, the infested area appears to be characterised by an abundance of hills, mountains, and fluvio-glacial valleys occupied by broadleaf mixed and coniferous forests, which is typical of the pre-alpine and alpine landscapes.

Fig. 1
figure 1

Map showing the study area. Red numbered lines departing from the presumed initial point of invasion represent the thirteen trajectories considered in this paper. Grey hexagons represent the cells providing P. japonica’s adult abundance data used in this paper to determine the rate of expansion along the trajectories. Lines coloured with a gradient of greys represent the monitored progression of the front of invasion from 2015 (lightest grey line) to 2021 (black line)

Full size image

The climate of the study area is temperate continental, with cold winters and warm summers and abundant (700–1400 mm/y) rainfall in spring and autumn. The study area was divided by the Regional Phytosanitary Service of Lombardy by superimposing a mesh of hexagonal spatial units (hereinafter, cells) with 1.42 km sides and an area of 5.42 km2 in order to characterise the expansion of the pest.

Data on population abundance and spread

To study the spread of P. japonica, we referred to adult population abundance data collected by the Regional Phytosanitary Service of Lombardy from 2015 to 2021. We classified cells as “infested” if they contained at least one detected adult individual. Adult abundance data were collected during the whole flight period using Trécé traps (USDA-APHIS). Traps were deployed between the end of May and the first days of June, corresponding to the beginning of the flight period of the species. The traps were baited with both the pheromone Japonilure and a mixture containing 3:7:3 blend of geraniol, eugenol and phenyl propionate. Traps were emptied twice a week during the peak adult flight period and once a week at both the beginning and end of the flight period. The monitoring season concluded between the end of September and the initial ten days of October, coinciding with the removal of traps from the fields. From 2015 to 2019, traps were positioned for mass trapping in infested cells. In 2020 and 2021, traps were employed solely for pest monitoring, with a primary emphasis on the perimeter of the affected zone. Nevertheless, several cells that had been infested since the first year were monitored to track the growth dynamics of the species. To investigate the population growth of P. japonica considering the infestation age of the cell (i.e. the number of years between the first year of infestation of the cell and the year of sampling), we excluded the cells that were sampled during the first year of the monitoring activities (i.e. 2015) as we were not able, for those cells, to define the year of first infestation. For each cell and for each sampling year, we estimated the daily average abundance of P. japonica adults during the flight period. The procedure involved regularisation of the trap capture data by an Erlang distribution (Curry and Feldman 1987). A description of the adopted procedure is available on the GESPO project website. For the sake of clarity, here is a concise description of the methodology: for each cell and each sampling year, the infestation age was calculated as the difference between the year of sampling and the first year in which that cell was declared infested, plus 1. This calculation resulted in a range from 1 (indicating the first year of infestation) to 6 (pertaining to cells sampled in 2021 that were infested since 2016). The final step involved estimating the average daily adult abundance in cells with the same infestation age, which was then used to determine the intrinsic population growth rate \({R}_{0}\) of P. japonica. To investigate the rate of expansion of the pest we relied on a subset of cells sampled by the Phytosanitary Service of Lombardy. This subset comprised 88 cells that intersected the 13 selected trajectories, as illustrated in Fig. 1. However, these 88 cells were sampled over multiple years, and those closer to the initial invasion point intersected more than one trajectory. This resulted in more sampling points per trajectory than the raw count of 88 cells suggests, totalling 224 cells (the sum of sampling points detailed in Table 1), ranging from a minimum of 8 sampling points for trajectory 4 to a maximum of 28 for trajectory 11. The average number of sampling points per trajectory is 17. Additionally, as each cell was set with multiple traps, the rounded average number of traps per cell per trajectory is also provided in Table 1.

Table 1 Values of diffusivity (D), rate of expansion (c*), Root Mean Square Error (RMSE) between observed and predicted relative abundance values, number of sampling points, average number of traps per cell, and maximum distance (in km) reached by the P. japonica in 2030 for each trajectory considered in the paper
Full size table

Modelling population growth and spread

We adopted a reaction–diffusion model (Kot et al. 1996; Kondo and Miura 2010) to describe the population growth in a given cell and the population diffusion among cells at discrete time steps (1 year). The reaction component of the model simulates population growth. The diffusion component describes the outcome of individual dispersal at the population level and it is achieved with a convolution operation, which combines two functions to express their transformation over time. The output of the model is the estimation of the adult population over the spatial dimension at the time t + 1 starting from the adult abundance at the time t. The model has been solved mathematically by an integrodifference equation with Gaussian kernel. From the numerical point of view, the population abundance has been approximated by continuous piece-wise high-order polynomial functions and the integrals have been computed by composite Gauss–Legendre quadrature formulas (Quarteroni et al. 2014). We defined the space of analysis as\(\Omega =\left[0,L\right]\), with \(L>0\) being the maximum diffusion distance possible in the model, and the period of analysis as\([0,{t}_{f}]\), with \({t}_{f}>0\) being the last year of simulation. The value \(\mathcal{N}(t,x)\) is the number of individuals located at the point \(x\in \Omega\) at time\(t\in [0,{t}_{f}]\). We discretised the period of analysis \([0,{t}_{f}]\) in \(s\in {\mathbb{N}}\) intervals\(, s\ne 0\), as follows: the time step is \(\tau =\frac{{t}_{f}}{s}\) and the partition of \([0,{t}_{f}]\) is\({t}_{0}=0,{t}_{1}=\tau ,\dots ,{t}_{j}=\tau *j,\dots ,{t}_{s}=\tau *s={t}_{f}\). Let \({\mathcal{N}}_{0}:=\mathcal{N}\left(0,\bullet \right)\) be the initial distribution of individuals in \(\Omega\) at time\({t}_{0}=0\). Then, for every discrete time\({t}_{j}\), we looked for a function.

$${\mathcal{N}}_{j+1}\left(x\right)={\int }_{\Omega }\kappa \left(x-y\right)g\left({\mathcal{N}}_{j}\left(y\right)\right)dy,$$

where the kernel \(\kappa \left(x-y\right)\) represents the diffusion process (from \(y\) to \(x\)) that only depends on the distance between the positions \(x\) and \(y\), while the function \(g\) is the reaction component of the model that describes population growth over time.

For the function \(g\) we used the Beverton–Holt function

$$g\left({\mathcal{N}}_{j}\right)=\frac{{R}_{0}K{\mathcal{N}}_{j}}{K+\left({R}_{0}-1\right){\mathcal{N}}_{j}} .$$

The parameter \({R}_{0}>0\) is the intrinsic growth rate; while, \(K\) is the population carrying capacity of the environment. Both the values of the intrinsic growth rate (i.e. the proliferation rate per generation) \({R}_{0}=2.70233\) and the local carrying capacity \(K=233\) adults per trap per day used in the present study were obtained from the results presented on the GESPO project website. Subsequently, the abundance data used in the present study were normalised to the carrying capacity.

For the kernel \(\kappa \left(x-y\right)\) we chose the gaussian probability density function

$$\kappa \left(x-y\right)=\frac{1}{\sqrt{2\pi {\sigma }^{2}}}{e}^{\frac{-{\left|x-y\right|}^{2}}{\left(2{\sigma }^{2}\right)}}$$

where \({\sigma }^{2}\in {R}^{+}\) is the variance of the normal distribution. \(D={\sigma }^{2}/2\) is the diffusivity coefficient, depending on the variance of the normal distribution and is an indirect measure of the capacity of the pest to move across the landscape.

We estimated the parameter \({\sigma }^{2}\) of the reaction–diffusion model in order that the numerical solution fits, in the least square sense, the mean abundance per year of infestation from 2015 to 2021. Then we projected the spread of P. japonica until \({t}_{f}=2030\). The computed values of \({\sigma }^{2}\) can be extracted from Table 1 bearing in mind that \(D={\sigma }^{2}/2\). To estimate the growth component, we fitted the Beverton–Holt function with the population abundance mean values calculated for the cells of the same age. To estimate the rate of expansion we used the mean abundance per year of infestation and fitted it to the reaction–diffusion model solution in order to predict a fitted value of \(D\), necessary to derive the rate of expansion of the pest. To estimate the rate of expansion of the pest, we adopted the approach proposed by Lutscher (2019) based on the asymptotic spreading speed scenario (i.e. the speed at which the leading edge of a population wave moves over time): we defined the rate of expansion \({c}^{*}\) as follows:

$${c}^{*}=\sqrt{2D log \left({R}_{0}\right).}$$

\({c}^{*}\) corresponds to the value at which any observer travelling at a speed higher than \({c}^{*}\), eventually finds itself beyond the population invasion front. Conversely, if an observer travels at speeds lower than \({c}^{*}\), it eventually becomes surrounded by the population (Lutscher 2019). This model can be applied with any initial distribution, provided that the moment-generating function of the kernel \(\kappa\) exists (true for the gaussian kernel) and that the growth function \(g\) is sublinear (i.e. \(g\left(\mathcal{N}\right)\le {g}{\prime}\left(0\right)\mathcal{N}\)). In our case, this latter assumption was verified as we assumed that the population was not subject to the Allee effect (e.g. the density-dependent effect suffered from establishing populations during the process of invasion) (Liebhold and Tobin 2008). All the fitting procedures were performed using the function lsqcurvefit of MATLAB (The Mathworks Inc 2022).

We assessed the model simulation performance by computing the Root Mean Square Error (RMSE) between the observed values (coloured dots in Fig. 2) and the corresponding predicted values (solid lines) of the relative abundance per year of infestation for each trajectory. As the calculated abundances are relative to the carrying capacity (and hence reach a maximum value of 100%), the RMSE represents the percentage of error in the model.

Fig. 2
figure 2

Spatio-temporal population dynamics of P. japonica along three emblematic trajectories: a trajectory number 11; b trajectory number 10; c trajectory number 9. The x axis represents the distance in km from the initial invasion point. The y axis represents the normalised population density (i.e. the percentage of the carrying capacity). Coloured dots represent the data used for the calibration of the model. Coloured lines represent the adult abundance curves. Solid lines represent the years from 2015 until 2021, for which abundance data were available, while dashed lines are predictions for the years from 2022 to 2030

Full size image

Effect of land use on the spread of P. japonica

To evaluate the influence of the landscape on the rate of expansion of P. japonica, we analysed the spread of the species in the study area along different radial trajectories (Fig. 1) originating from the initial point of establishment in Italy. The exact location of this point is currently not known. Thus, in our work we assumed the invasion starting point as the centroid of the infested area (encompassing both the Lombardy and Piedmont Regions) in 2014, which includes the municipalities of Pombia, Oleggio, Cameri, Galliate Marano Ticino, Bellinzago Novarese in Piedmont and Ferno, Lonate Pozzolo, Nosate, Turbigo and Vizzola Ticino in Lombardy. Starting from the centroid of the initial infested area, we hypothesised that the pest spread in all directions with a speed assumed to be dependent on land use. To test this hypothesis, we considered 13 trajectories in order to cover the entire expansion area in Lombardy and intercept the highest number of cells for which adult abundance data were available. The trajectories span over an angle of 160° and are relatively equidistant from one another with a mean angle around 13° (Fig. 1).

Data on land cover refer to all the cells intersected by each trajectory, starting from the initial invasion point, and finishing at the farthest occupied cell. Data were obtained from the third level of the regional Use of Agricultural Land and Forestry cartography (Destinazione d’Uso dei Suoli Agricoli e Forestali–DUSAF 6) (Regione Lombardia 2018). Based on the work of Borgogno-Mondino et al. (2022) we classified the land use classes in two categories: suitable and unsuitable. We deemed as suitable for the pest all the urban green and recreational areas, meadows, arable land, permanent cultures, broadleaf forests and mixed (broadleaf and coniferous) forests. The inclusion of broadleaf and mixed forests among the suitable land uses is based on a comparison between the silvicultural map of Lombardy (https://www.geoportale.regione.lombardia.it) and the EPPO host species list (https://gd.eppo.int/taxon/POPIJA/hosts), revealing a significant overlap in species. Coniferous forests (Fox 1927), heavily urbanised areas without green spaces, production areas (i.e. factories and commercial and artisanal productive facilities) and infrastructures were considered unsuitable. For each cell, we computed a land cover suitability index, calculated as the percentage of area covered by suitable land cover classes. Then, we averaged the suitability index of all the cells belonging to each trajectory. To investigate the role of land use on the spread of the pest along the different trajectories we relied on a linear model (lm function in R Core Team 2023) linking the rate of expansion calculated along each trajectory to its land cover suitability index.

Results

Pattern of population dynamics in space and time

The diffusion component for the Gaussian kernel was computed for the 13 different trajectories (see Table 1). The diffusivity \(D\) ranged from 0.05 to 0.48, with a mean of 0.19 ± 0.15 SD. The estimated \({c}^{*}\) values varied greatly between trajectories and resulted in a three-fold increase from the slowest (4.5 km/y for the 10th trajectory in the south eastly direction) to the quickest trajectories (13.7 km/y for the 4th trajectory in the north-east direction; 13.8 km/y for the 13th trajectory in the southerly direction). The mean and standard deviation of the values of \({c}^{*}\) were 8.2 km/y and 3.1 km/y, respectively. The RMSE ranged from 5.59 to 13.45, with a mean of 9.52 (please, refer to Table 1). The model predicted the pest to spread from 81 to 207 km by the year 2030 from the point of initial establishment to the distance corresponding to a normalised population density of 1% of the carrying capacity. Figure 2 shows the population abundance normalised to the carrying capacity K along three notable trajectories (to demonstrate a range of results generated by the model) as functions of the distance from the centroid and over time. The plots consider both the period in which we have observations (from 2015 to 2021) and a projection for the period 2022–2030. For the first period the simulated population abundance was comparable with the observations. The projections up to the year 2030 were performed considering a constant diffusivity D, therefore they disregard any variation in the parameter D due to landscape features beyond the last occupied cell along each trajectory.

The carrying capacity at the initial invasion point was reached by the 9th year of simulation (2023). Trajectory number 11 (Fig. 2a) showed one of the fastest dynamics, predicting that P. japonica could reach distances of over 169 km from the point of initial invasion by the end of the simulation (i.e. 2030). Figure 2b (trajectory number 10) shows a rather slow rate of expansion: the pest reached a distance of 81 km from the initial invasion point in 16 years of simulation. A case of intermediate rate of expansion is represented in Fig. 2c (trajectory number 9). After 16 years the model predicts the invasion front will reach a distance of 112 km from the starting point.

Effect of land cover on the rate of expansion

The land cover types varied greatly among the 13 trajectories (see Table 2): urbanized areas covered from 7.78% to 51.58% (with a mean value of 34.47 ± 14.19 SD) of the total area of the cells intersected by a specific trajectory; arable lands covered from 9.86% to 68.02% (mean 33.27 ± 19.34 SD) of the area intersected by the trajectories, showing a clear increasing trend from the trajectories north to south; permanent cultures (e.g. orchards) generally occupied a negligible fraction of the area of interest, although with notable presence (up to 9.18%) for the trajectories facing south; permanent meadows, instead, occupied around 5% (mean 4.52 ± 1.91 SD) of the area in all the trajectories; woodlands (without distinction on the species composition) represented around 50% of the occupied area for the northward trajectories and around 10% for those pointing eastward and southward (mean 23.10 ± 16.35 SD).

Table 2 Percentage of the main land cover classes (urban, arable land, permanent cultures, meadows and woodlands) and suitability index for each trajectory considered in the paper
Full size table

The land suitability index (see Table 2), computed as the percentage of suitable areas of the cells intersected by a trajectory, varied from 50 to 89% (mean 65% ± 13% SD). The three most suitable trajectories were those facing south, towards the middle of the Po Valley and yielded indices above 80%. The remaining trajectories had lower values of the suitability index, from a low of 50% to a maximum of 65%.

The linear model confirmed that the suitability of the land use for P. japonica impacts the rate of expansion (Fig. 3). The regression coefficient is equal to 0.258 km/y, indicating an increase in the rate of expansion of approximately 260 m/y per unit (percentage) of land cover suitability. The adjusted R-squared for the regression is 0.33, indicating a moderate goodness of fit of the model to the data.

Fig. 3
figure 3

Linear model (blue line) explaining the rates of expansion (asymptotic speed) in kilometres travelled by P. japonica per year given the land suitability index (expressed as percentage). Black points are the data for each of the 13 trajectories selected. The grey-shaded area represents the confidence interval

Full size image

Discussion

In this paper we presented a model that simulates the spatio-temporal population dynamics of P. japonica during an ongoing invasion process based on solid mechanistic foundations and biologically meaningful parameters. The discrete time Beverton–Holt model used for simulating the population growth (reaction component) includes two important parameters such as the carrying capacity (estimated as \(K=233\) adults per trap per day) of the environment, and the intrinsic growth rate of the population (\({R}_{0}\) = 2.70233). The nonlinear Beverton–Holt growth function accounts for the density-dependent regulation of the population (Lewis et al. 2016). These values have been estimated within a specific environmental context (northern Italy); therefore, we anticipate potential variations in both the growth rate and carrying capacity based on local conditions related to climate, soil, and plant communities characterising different areas. Furthermore, the relatively short period covered in our studies prevents the ability to highlight any reaction of the receiving biotic community (e.g. natural controls) that might result in a reduction of the environmental carrying capacity. Preliminary observations performed by authors of this paper (Ciampitti, Cavagna and Bianchi) seem to confirm the existence of a biotic reaction (e.g. parasites, predators, etc.) in the receiving community leading to a decrease in pest abundance a few years after reaching a local maximum. The candidate role of biological control agents in regulating Popillia japonica population abundance has been reported also in other countries as it is for example in the case of Ovavesicula popilliae (Smitley et al. 2022), although most probably the populations are regulated by a combination of biotic and abiotic factors (Potter and Held 2002).

In order to estimate the speed of invasion of the pest, we adopted the approach proposed by Lutscher (2019) for nonlinear growth functions without an Allee effect. This approach is more general than other approaches: the point-release scenario and the travelling front scenario, both require a linear growth function, and either a Dirac delta initial distribution (the former) or an initial developed front (the latter). Our selection of Lutscher’s approach was primarily driven by the fact that our model incorporates the nonlinear growth Beverton–Holt function, a Gaussian diffusive kernel, and an initial Gaussian distribution that is not necessarily fully developed. This choice stems from the demonstration provided by Lutscher (2019) that the asymptotic speed of invasion of the nonlinear approach equals the minimum speed of the travelling front scenario and the minimum invasion speed of the point-release scenario (both easily computable), thus we adopted that approach for estimating the speed of invasion of P. japonica.

The high variability in the rates of expansion reported in our study (between 4.5 and 13.8 km/y, with a mean of 8.2 km/y) is in line with the results of Fox (1927). Considering the first 10–15 years of the pest’s invasion in the US, Fox (1927) observed high variability in the mean rate of expansion along the cardinal, intercardinal and secondary cardinal directions from Riverton, New Jersey (i.e. the point of first introduction in the US). The author discussed the possible causes of such variability, addressing wind conditions, orography, human-mediated long-distance dispersal, and land cover. Considering the average rate of expansion, our results are remarkably close to those reported in Borgogno-Mondino et al. (2022) and in various studies from the US. Concerning the Italian case, Borgogno-Mondino et al. (2022) proposed an iterative model based on a geostatistical approach to describe the spread of P. japonica. They obtained variograms showing spatial autocorrelation between 7.5 and 15 km. Their analysis represents a first attempt to integrate population growth, individual dispersal, and natural death using a geostatistical modelling approach. Regarding the spread of P. japonica in the US, rates of expansion derivable from Hadley and Smith (1923) vary from less than 1 km/y to10 km/y, while those reported in Allsopp (1996) are of 7.7 km/y for the second decade from the beginning of the invasion and 11.9 km/y for the third decade. Both the above papers, along with Lewis et al. (2016), mention an increase in the rate of expansion of the pest in the first years of invasion, becoming stable after the seventh (Lewis et al. 2016) or the tenth (Allsopp 1996) year. The increasing rate of the range expansion in the initial phases of invasion in the US poses the question whether it is acceptable to represent the expansion in Lombardy Region as a linear process, since the establishment occurred recently. The assumption of a linear rate of expansion appears reasonable considering that the first report of the presence of the insect in northern Italy (in 2014) most probably does not coincide with the date of the actual establishment of the pest in the region (EFSA 2019). Thus, the population could have gone unnoticed during the initial expansion phase.

A more comprehensive model of the spread of P. japonica should include a component simulating the LDD of the pest, but the lack of data prevents the estimation of reliable distributions describing the probability and pattern of LDD events. Human-assisted LDD events are mainly due to hitch-hiking on transported commodities but seem to play a secondary role in P. japonica spread dynamics in northern Italy and continental Europe. However, incursions of individuals that did not originate established populations in Italy have been reported in Parma and Piacenza in 2020 and in 2021, near the main airport in Sardinia in 2021, in the Valle d’Aosta Region in 2021 (EPPO 2022b), and in the Friuli Venezia Giulia Region in 2021 (Bassi et al. 2022). Incursions without establishment have also been reported for other European countries, such as the Netherlands in 2018 (EPPO 2019) and Germany in 2021 and 2022 (EPPO 2022c). By now, only one LDD event is known to have led to an established population in Europe, i.e. the recently reported established population near the airport of Zurich. The limited amount of incursion reports of the pest might be due to the effectiveness of the measures adopted for the confinement of the pest, such as limitations imposed on the transportation of propagative material from the infested areas and the inspection of high-risk areas such as plant nurseries, airports, and camping areas (Regione Lombardia 2023). Moreover, the above-mentioned detections all occurred nearby airports, railroads and motorways. Such conditions seem also to confirm that wind-mediated LDD events are uncommon, as already stated in Borgogno-Mondino et al. (2022). However, we expect that the expanding range of the pest could elevate the risks linked to LDD events, consequently heightening the significance of these events in influencing the spatio-temporal dynamics of P. japonica in continental Europe in the coming years.

Investigating the role of landscape on the rate of expansion of invasive pest insects is crucial for planning appropriate management strategies. The land suitability index presented in this study summarised the effects of land cover on the rate of expansion of P. japonica populations in a single, intuitive value. Trajectories 4, 11, 12 and 13 are the fastest. The first one points north-east, while the last three all point south-southeast towards the Po Valley in a landscape typically devoted to intensive, irrigated agriculture, which offers plenty of feeding and nesting sites for the species and with no obstacles to the movement of the pest. Moreover, trajectories 11, 12 and 13 show a high percentage of permanent cultures when compared to the other trajectories, suggesting that this land use might play a role in determining the rate of expansion of the pest. Trajectory 4, instead, has a low suitability index and a low percentage of land use dedicated to permanent cultures. Nevertheless, it yields one of the highest rates of expansion. The reason behind this counterintuitive result could lie in the specific landscape context of trajectory 4, where the lake of Como and the Prealps act as barriers to the undirected expansion of the pest. This could result in the convergence of most dispersing individuals along the aforementioned trajectory. The trajectories with the lowest rates of expansion are number 2, 5 (north–west trajectory) and 10 (south–west trajectory). All of these are characterised by a medium to high level of urbanised areas with no vegetation (the percentage of urban land cover ranges between 38 and 42%). For trajectories 2 and 5 orography could have had a role in determining the dynamics of expansion of the pest; these trajectories point towards the Alps, which offer a physical barrier (both in terms of elevation and in terms of decreasing temperatures with height) to the dispersal of individuals. Taken together, the results of the regression model linking the rates of expansion with the suitability indexes suggest that other variables might be relevant in determining the rate of expansion of the pest. For example, Simonetto et al. (2022), demonstrated the influence of environmental variables other than land use on P. japonica population dynamics. Included among them are soil organic carbon, soil texture, soil and air temperature, rainfall, and relative humidity.

In conclusion, our model was able to accurately predict the spatio-temporal dynamics of P. japonica populations and allowed the estimation of the spread of the pest, including the role of land use in the Lombardy region. The mathematical structure of the model proposed is highly flexible and allows room for fine-tuning the parameters and modification of the diffusion kernel as new knowledge and data become available. Our results provide fundamental insights on the growth pattern and spread of P. japonica in continental Europe. Knowledge on the landscape-specific growth and rate of expansion and the factors influencing them allows the use of this model as a tool supporting the design and the implementation of appropriate containment and, where possible, eradication actions, especially in high-risk areas.

Author contributions

Conceptualisation contributed by GG, GS, AS, PG, ADF, and AnBa; data curation contributed by GS, MC, BC, and AlBi; formal analysis contributed by PG; investigation contributed by MC, BC, and AlBi; methodology contributed by GG, PG, AS, GS, and ADF; software contributed by PG; supervision contributed by GG; writing–original draft contributed by ADF, GS, and PG; writing–review and editing contributed by GG, GS, AS, MC, BC, AlBi, AnBa, NM, ADF, and PG