Key message

  • There are high concerns linked to the potential establishment and impacts of Spodoptera frugiperda in Europe

  • We developed a physiologically-based model to quantitatively assess the risks linked to S. frugiperda

  • Risk of establishment of S. frugiperda was predicted in the Mediterranean coastal areas of Southern Europe

  • Risks linked to transient populations (i.e. populations migrating from suitable areas) were predicted in Southern and Central Europe

Introduction

The fall armyworm, Spodoptera frugiperda (J.E. Smith) (Lepidoptera: Noctuidae) is a phytophagous pest considered a major threat to agricultural production and food security (Early et al. 2018; Tambo et al. 2021), especially in developing countries (Devi 2018; FAO 2020; Suby et al. 2020; Koffi et al. 2020). The species is known to feed on more than 350 host plants including economically valuable crops such as maize, rice, soybean, sorghum, wheat, barley, and cotton (de Freitas Bueno et al. 2011; Hardke et al. 2015; Montezano et al. 2018). Impacts on crops are caused mainly by late instar larvae (Overton et al. 2021) feeding on stems, branches, leaves, and reproductive structures of the host, and causing direct yield loss, defoliation, and general weakness of the plant (Harrison 1984; Vilarinho et al. 2011). The larval trophic activity might favour plant infection caused by fungi (Farias et al. 2014). Yield loss to maize ranges from 11 to 67% (Hruska and Gould 1997; Day et al. 2017; Kumela et al. 2019; Baudron et al. 2019). Reported average losses for other crops are 26% for sorghum, 24% for sweet corn, 13% for bermudagrass, and 5% for rice. The control measures for protecting crops affected by the species and restrictions on the trade of potentially infested products cause further economic and social costs (Overton et al. 2021). The species is native to tropical and subtropical America where it is considered a prevalent pest for maize, soybean, cotton, and other major crops (Nagoshi et al. 2007; Baudron et al. 2019; Koffi et al. 2020). Human-mediated transportation and trades (Cock et al. 2017), the high migratory capacity (the species might fly up to 100 km per night) (Rose et al. 1975; Westbrook et al. 2016), and the high prolificacy (more than 1500 eggs laid per moth) (Luginbill, 1928) of the species facilitated the dispersal of the pest in non-native areas. In 2016, the species was accidentally introduced in Central and Western Africa (Goergen et al. 2016) where it was able to spread in vast areas of sub-Saharan and North Africa (Day et al. 2017; Cock et al. 2017; EPPO 2020a). Since 2018, the species invaded vast areas of the Middle East (EPPO 2019a; EPPO 2020b, c, d), South Asia (EPPO 2018; Sharanabasappa et al. 2019), South-Eastern Asia (EPPO 2019b, c; EPPO 2020e; Sartiami et al. 2020; Zaimi et al. 2021), East Asia (EPPO 2019d; Suh et al. 2021), North-Eastern Asia (EPPO 2019e) and Oceania (EPPO 2020f). In Europe, the species is currently (FAO 2021) present in the Canary Islands (EPPO 2021). S. frugiperda is on the EPPO A2 list of quarantine pests, and on the European Commission list of priority pests (EU 2019) due to the risk of introduction, establishment, and consequences of this pest to Europe. Fresh plant products imported from Latin America represent the main pathway of entry of the species into the EU (EFSA PLH Panel et al. 2017; 2018a; EFSA et al. 2020). Another pathway of introduction is represented by the possibility of eggs and adults to entry as hitchhikers on international flights (Early et al. 2018). The high migratory ability of the species causes concerns about the potential impacts of transient populations moving from hotspots to new areas during the favourable season (EFSA PLH Panel et al. 2018a; Timilsena et al. 2022). A realistic threat is the introduction of individuals from North Africa to Europe due to natural or wind-mediated dispersal (Westbrook et al. 2016; Early et al. 2018).

Given the potential threats of S. frugiperda to European agriculture, it is fundamental to quantitatively estimate the risk of establishment and the potential impacts linked to the species. This information is fundamental for planning and implementing surveillance and inspections to reduce the likelihood of introduction and establishment of the pest in Europe (EFSA PLH Panel et al. 2018a; EFSA et al. 2020). So far, many species distribution models have been developed for predicting the potential habitat suitability for S. frugiperda (Ramirez-Cabral et al. 2017; Du Plessis et al. 2018; Early et al. 2018; Liu et al. 2020; Baloch et al. 2020; Fan et al. 2020; Zacarias 2020; Huang et al. 2020; Tepa-Yotto et al. 2021; Ramasamy et al. 2021). However, there is high uncertainty on the risks linked to the establishment of this pest in Europe. For instance, no risk of establishment but only risk linked to transient populations was predicted by Du Plessis et al. (2018). On the contrary, suitable areas in Southern Europe were identified by Early et al. (2018) and EFSA PLH Panel et al. (2018a). Other authors identified risks in areas of Central Europe (Zacarias 2020) or further north, up to Ireland (Liu et al. 2020; Ramasamy et al. 2021) and Southern Norway (Tepa-Yotto et al. 2021), although with low habitat suitability indices. This high uncertainty reflects the need to establish sound criteria and reliable models for obtaining a realistic assessment of the risk linked to a pest (Ponti et al. 2015).

In this work, we aimed at providing a solid and quantitative assessment of the risks linked to S. frugiperda in Europe through the application of a physiologically-based (i.e. mechanistic) modelling approach. This approach allows for the faithful description of important aspects of the biology of the species (Sparks 1979), such as the nonlinear responses to temperature and the influence of relevant abiotic drivers (density-dependent factors, mortality due to biotic agents) on the individual physiology, population distribution and dynamics (Régnière et al. 2012a; Gutierrez and Ponti 2013). The model was used to respond to the following assessment questions (AQ), which are highly relevant for estimating the risks linked to S. frugiperda in Europe (EFSA PLH Panel et al. 2018b). AQ 1—Is the model able to predict the pattern of population dynamics and the limits of establishment in the area of current distribution? (current distribution and dynamics); AQ 2—Can the species establish in Europe? If yes, what is the area of potential establishment of the species? (establishment in Europe); AQ 3—What is the population dynamics of the species in the areas of potential distribution in Europe? (population dynamics in Europe); AQ 4—Can the species originate transient populations in Europe? If yes, can population abundance in transient populations represent a risk for cultivated plants? (dynamics of transient populations).

Materials and methods

The model

In this work, we developed a physiologically-based model using a system of Kolmogorov partial differential equations to simulate the stage-specific population dynamics of S. frugiperda considering the two dimensions of time \(t\) and physiological age \(x\) (Buffoni and Pasquali 2007; Rafikov et al. 2008; Solari and Natiello 2014; Lanzarone et al. 2017) (full mathematical details of the model are present in Section S1 of supplementary material 1). We assumed that a population of S. frugiperda is composed by four stages \(i\), namely egg (\(i = 1\)), larva (\(i = 2\)), pupa (\(i = 3\)), and adult (\(i = 4\)). Physiological age in the \(i\)-th stage, \(x^{i} \in \left[ {0,1} \right]\) represents the level of development of an individual in the stage \(i\) (Buffoni and Pasquali 2007). With \(x^{i} = 0\), we represent an individual at the beginning of \(i\)-th stage, while with \(x^{i} = 1\) we represent an individual at the end of the \(i\)-th stage. The term \(\phi^{i} \left( {t, \;x} \right)\) represents the number of individuals in stage \(i\) at time \(t\) with physiological age \(\left[ {x,\;x + {\text{d}}x} \right]\). The overall number of individuals in stage \(i\) at time \(t\) is calculated as \(N^{i} \left( t \right) = \mathop \smallint \limits_{0}^{1} \phi^{i} \left( {t, x} \right){\text{d}}x\). The population abundance in the stage \(i,\;N^{i} \left( t \right)\), is defined by the number of individuals in a spatial unit as defined in ‘Definition of the spatial unit’ section. In the present work, we considered the predicted adult population abundance (reported as the average number of moths per trap per week) as a descriptor of the potential impacts of S. frugiperda. The simulations were performed using MATLAB version R2018a (MATLAB, R2018a, The MathWorks, Inc., MA, USA). We assumed the population dynamics of S. frugiperda was dependent on the species’ life-history strategies. These were described at the individual level by stage-specific development, mortality, and fecundity rate functions. Since temperature is considered one of the main variables influencing the physiology of poikilotherms (Gutierrez 1996; Régnière et al. 2012b; Gilioli et al. 2021a), the effects of the time-dependent temperature profile \(T\left( t \right)\) affecting the species’ life-history strategies were considered in the model (Barfield et al. 1978; Silva et al. 2017; Du Plessis et al. 2020).

Development rate function

We defined \(\nu^{i} \left( {T\left( t \right)} \right)\) as the temperature-dependent development rate function of individuals in stage \(i\) as a function of temperature \(T\left( t \right)\). For the stages \(i = 1,2,3\) we used the development rate functions that are defined in Gilioli et al. (2021b). For the stage \(i = 4\), respect to Gilioli et al. (2021b), we increased the life-span of the adults by reducing the development rate function \(\nu^{4} \left( {T\left( t \right)} \right)\) by a fixed factor of 2.5 to obtain more realistic adult survival curves (He et al. 2021a; Zhang et al. 2021). The methodology used for estimating parameters of the development rate function \(\nu^{i} \left( {T\left( t \right)} \right)\) is presented in Section S1.1 of supplementary material 1.

Mortality rate function

As in Gilioli et al. (2021b), we assumed the mortality rate function \(m^{i} \left( t \right)\) for the stages \(i = 1,\; 3, \;4\) depending on temperature according to the following law

$$m^{i} \left( t \right) = \mu^{i} \left( {T\left( t \right)} \right), \;i = 1,\; 3,\; 4$$

with \(\mu^{i} \left( {T\left( t \right)} \right)\) being the temperature-dependent instantaneous mortality acting on individuals within each stage at time \(t\) (see Section S1.2 of supplementary material 1). The mortality of larvae is affected by multiple factors, such as weather conditions (Varella et al. 2015), the attack of biotic agents (e.g. predators, parasites and pathogens) (Escribano et al. 2000; Zanuncio et al. 2008), and density-dependent factors (e.g. cannibalistic behaviour) (Chapman 1999; Chapman et al. 2000; Andow et al. 2015; He et al. 2021b). To account for these factors, in the present work, the mortality rate function for the larval stage \(m^{2} \left( t \right)\) is expressed as follows

$$m^{2} \left( t \right) = \mu^{2} \left( {T\left( t \right)} \right)\left( {1 + \alpha \left( {\frac{{N^{2} \left( t \right)}}{\gamma }} \right)^{2} } \right) + \beta$$

with \(\left( {1 + \alpha \left( {\frac{{N^{2} \left( t \right)}}{\gamma }} \right)^{2} } \right)\) representing a density-dependent component simulating the intraspecific competition (e.g. cannibalistic behaviour), \(\alpha > 0\) representing a multiplicative term, and \(\beta > 0\) representing a biotic component simulating the role of predators, parasites and pathogens.

The parameter \(\gamma > 0\) represents the larval carrying capacity based on resources availability. The parameter \(\gamma\) was set to 3000 which corresponds to the larval abundance at the carrying capacity in the spatial unit considered in the present study (see ‘Definition of the spatial unit’ section for details). The methodology used for estimating the temperature-dependent component of the mortality rate function \(\mu^{i} \left( {T\left( t \right)} \right)\) is presented in Section S1.2 of supplementary material 1. Parameters \(\alpha\) and \(\beta\) were estimated through the calibration procedure (‘Model calibration’ section).

Fecundity rate function

For the adult stage, we defined the fecundity rate function \(F^{1} \left( t \right)\) representing the production of eggs by adult females (Johnson 1987). As in Gilioli et al. (2021b), the fecundity rate function depends on female age and temperature. In the present work, we further introduced a density-dependent regulation term to account for the role of intraspecific competition in egg production due to limitations in the per-capita food supply (Leather 2018). The fecundity rate used in the present study is

$$F^{1} \left( t \right) = g\left( {T\left( t \right)} \right)\left( {1 - \frac{{N^{4} \left( t \right)}}{{S + N^{4} \left( t \right)}}} \right)\mathop \smallint \limits_{0}^{1} \phi^{4} \left( {t,x} \right)h\left( x \right){\text{d}}x$$

with \(g\left( {T\left( t \right)} \right)\) describing the temperature-dependent component, \(\left( {1 - \frac{{N^{4} \left( t \right)}}{{S + N^{4} \left( t \right)}}} \right)\) describing the density-dependent component, \(\phi^{4} \left( {t,x} \right)\) being the number of adult individuals at time \(t\) and physiological age \(x\), and \(h\left( x \right)\) describing the physiological age-dependent component. The terms \(g\left( {T\left( t \right)} \right)\) and \(h\left( x \right)\) were taken from Gilioli et al. (2021b) (see Section S1.3 of supplementary material 1 for details). The term \(S\) is a half-saturation term in the density-dependent regulation of female fecundity. Based on the assumption that the adult abundance at the carrying capacity in the spatial unit is \(N_{K}^{4} = 320\) adult individuals per week (see ‘Definition of the spatial unit’ section for details), we set \(S = 0.5 N_{K}^{4} = 160\). With this set up, the density-dependent term \(\left( {1 - \frac{{N^{4} \left( t \right)}}{{S + N^{4} \left( t \right)}}} \right)\) is almost 1.00 (negligible density-dependent effects) for adult population abundances lower than 10 individuals, and almost 0.35 (relevant density-dependent effects) for adult population abundances approaching 320 individuals per trap per week.

Model calibration

The calibration procedure consisted in estimating the parameter \(\alpha_{j}\) and \(\beta_{j}\) that will be used for the definition of the parameters \(\alpha\) of the density-dependent mortality term and the biotic mortality term \(\beta\) included in the mortality rate function of the larval stage \(m^{2} \left( t \right)\). Parameters \(\alpha_{j}\) and \(\beta_{j}\) were estimated by minimising the mean squared distance between the simulated and the observed adult abundance for each of the 21 observation datasets \(j\) representing the calibration dataset (see ‘Data on pest population dynamics’ section). The minimisation was performed for each of the 21 observation datasets \(j\) through solving the following function

$$Q_{j} \left( {\alpha_{j} ,\beta_{j} } \right) = \mathop \sum \limits_{j = 1}^{21} \frac{1}{{R_{j} }}\mathop \sum \limits_{i = 1}^{{R_{j} }} \left| {N_{j}^{4} \left( {t_{i} ;\alpha_{j} ,\beta_{j} } \right) - A_{j} \left( {t_{i} } \right)} \right|^{2}$$

The term \(A_{j} \left( {t_{i} } \right)\) represents the observed adult abundance in the dataset \(j\) at the time \(t_{i}\) corresponding to the time at which adult abundance was sampled. The term \(R_{j}\) represents the number of sampled data available for each dataset \(j\). With \(N_{j}^{4} \left( {t_{i} ;\alpha_{j} ,\beta_{j} } \right)\), we define the adult abundance in the dataset \(j\) at time \(t_{i}\), obtained by solving the Kolmogorov equations with the parameters \(\alpha = \alpha_{j}\) and \(\beta = \beta_{j}\) keeping fixed the other parameters. The optimal parameters \(\hat{a}_{j}\) and \(\hat{\beta }_{j}\) were the minimisers of the \(Q_{j}\), i.e. they allow for the minimum difference between simulated and observed adult population abundance

$$Q_{j} \left( {\hat{\alpha }_{j} ,\hat{\beta }_{j} } \right) = \min_{{\alpha_{j} ,\beta_{j} }} Q_{j} \left( {\alpha_{j} ,\beta_{j} } \right)$$

For the minimisation procedure, we used the MATLAB function fmincon with step tolerance equal to 10–5 for the stopping test.

Simulation design

The population dynamics model of S. frugiperda was used to explore the four assessment questions reported in ‘Introduction’ section. To account for the uncertainty linked to the estimates of parameters \(\alpha\) and \(\beta\), the model was implemented considering three assessment scenarios (see ‘Generation of assessment scenarios’ section).

Assessment question 1— Current distribution and dynamics

The capacity of the model to predict the local population dynamics of S. frugiperda was tested by comparing simulated and observed adult population abundance using data obtained in three locations selected along a latitudinal gradient in the area of current distribution in North America (see ‘Data on pest population dynamics’ section). From south to north, we considered a highly suitable location (Miami Dade County, Florida), a location at the edge of the area of establishment (Alachua County, Florida), and a location that is currently known to be reached only by migrating populations (Tift County, Georgia) (Westbrook et al. 2016; Garcia et al. 2018). The population dynamics were simulated using the temperature profile of the current climate in the tested locations as input data (see ‘Temperature data’ section). Initial conditions were set to 5 pupae uniformly distributed in their physiological age (from 0 to 1) on the 1st of January. The model was implemented for four consecutive years, repeating the same yearly temperature profile, to obtain stable population dynamic patterns and model outputs that were independent of the initial conditions. We assumed that no migration of individuals was possible from and to each location in which the model was implemented. The assessing variables considered were the yearly average number of moths per trap per week, the number of generations per year, and the maximum adult population abundance reached over the last year of simulation.

Assessment question 2— Establishment in Europe

For assessing the potential distribution and abundance of S. frugiperda in Europe, we implemented the model in a spatial grid of 0.1° × 0.1° representing the European territory (see ‘Temperature data’ section). In each node of the grid, the population dynamics was assessed using the same initial conditions defined in AQ 1 and the temperature profile of that specific node (Gilioli et al. 2014, 2021c; Pasquali et al. 2020). The species was considered established in a node if, at the first time-step of January 1st of the last year of simulation, the adult abundance was higher than an adult abundance threshold \(\left( {A_{0} = 0.01} \right)\). The threshold \(A_{0}\) was set by considering the average of the minimum population abundance reached by the species at the northernmost edge of the area of establishment in a set of locations in North America, including the location of Alachua County (Florida) tested in AQ 1. Species’ potential distribution was estimated over the last year of simulation. The area of potential establishment of S. frugiperda in Europe was given by the set of grid nodes where the species was considered established.

Assessment question 3— Population dynamics in Europe

The local population dynamics of S. frugiperda in Europe was assessed by implementing the model in 3 locations using the initial conditions explained in AQ 1. Locations were chosen based on the simulated S. frugiperda potential dynamics in Europe obtained by answering AQ 2. Based on the model’s result, a highly suitable location was selected in Cyprus, and two less suitable locations were selected, in Southern France and on the Atlantic coast of Portugal. We considered the same assessing variables presented in AQ 1.

Assessment question 4— Dynamics of transient populations

Transient populations are analysed in a hypothetical scenario in which migrating adults arrive in a location characterised by temporary suitable conditions (e.g. warm temperature conditions during spring or summer), but where the species is not able to survive during fall or winter. To assess the dynamics of transient populations, we simulated the introduction of an inoculum characterised by five adult individuals uniformly distributed in their physiological age (from 0 to 1) in four maize production areas in Europe, outside the predicted area of establishment: Rădoiești (Romania, 44th parallel north), Ghedi (Italy, 45th parallel north), Ouarville (France, 48th parallel north), and Engelsberg (Germany, 48th parallel north). The dynamics of transient populations was assessed considering three different Days of the Year (DOY) for the introduction of the inoculum: April 1 (90th DOY), June 1 (150th DOY), and August 1 (210th DOY). The model was implemented from the date of introduction of the inoculum to the end of the year, using as temperature profile the current climate in the tested location (see ‘Temperature data’ section). The assessing variables considered were the average number of moths per trap per week, the number of generations, and the maximum adult population abundance over the simulation period. We assumed that the inoculum was not able to originate a transient population if the predicted adult population abundance reached values below or equal to the adult abundance threshold \(A_{0}\) during the simulation period.

Generation of assessment scenarios

Considering the range of distribution of the parameters \(\alpha_{j}\) and \(\beta_{j}\) estimated through the calibration procedure (see ‘Model calibration’ section), we calculated the 10th, the 50th, and the 90th quantiles of the distributions for the definition of parameters \(\alpha\) and \(\beta\). To account for variability in the population dynamics, we generated 9 different assessment scenarios, combining the quantiles of \(\alpha\) and \(\beta\). In the present study, we consider the worst-case assessment scenario where the species has lower mortality \(\left( {\alpha = 10{\text{th; }}\beta = 10th} \right)\), The median-case assessment scenario, obtained considering the medians of parameters distribution \(\left( {\alpha = 50th;\;\beta = 50th} \right)\), and the best-case assessment scenario, where S. frugiperda mortality is high \(\left( {\alpha = 90th;\;\beta = 90th} \right)\). The values of parameters related to the three investigated scenarios are reported in Table 1. The population dynamics of S. frugiperda in the current area of establishment in North America (AQ 1) and the dynamics of transient populations in Europe (AQ 4) were predicted considering the median-case assessment scenario. The best-case, the median-case, and the worst-case assessment scenarios were considered for predicting the potential distribution of S. frugiperda in Europe (AQ 2) and the population dynamics of the pest within the predicted area of establishment (AQ 3).

Table 1 Estimates of parameters \(\alpha\) and \(\beta\) linked to larval mortality for the best-case, the median-case, and the worst-case assessment scenario discussed in the present study

Data

Data on pest population dynamics

Data on pest population dynamics were used for estimating parameters in the function describing the larval mortality (see ‘Model calibration’ section) and to test the model’s capacity to predict the population dynamics patterns and the establishment of S. frugiperda in North America (AQ 1). Population dynamics data used for calibration purposes (hereinafter, calibration dataset) refer to 21 time-series adult trap catches data collected in the area of establishment in Central and North America from 1982 to 2019, selected for their quality in terms of completeness of the time-series and realism of population trends (see supplementary material 2) (Silvain and Ti-A-Hing 1985; Pair et al. 1986; Nagoshi and Meagher 2004; Meagher and Nagoshi 2004; Rojas et al. 2004; Salas-Araiza et al. 2018; Salazar-Blanco et al. 2020). The calibration dataset covers latitudes between 4.85 and 28.76 parallel north and it includes, 1 time-series data collected in French Guyana (Matoury), 1 in Costa Rica (Guanacaste), 3 in Mexico (Manzano and Irapuato), and 16 in the United States (from Southern to Northern Florida). Additional population dynamics data were used for answering the AQ 1 and refer to time-series adult trap catches collected in three locations: Miami Dade County (Florida, 25th parallel north), Alachua County (Florida, 29th parallel north), and Tift County (Georgia, 31st parallel north) (Pair et al. 1986; Meagher and Nagoshi 2004; Garcia et al. 2019).

Definition of the spatial unit

The simulated adult abundance variable used in our model \(N^{4} \left( t \right)\) refers to the number of adult individuals caught in a trap per week. To consistently allow the comparison between observed and simulated adult population abundance, the temporal unit of the population dynamics data used in the present study was referring to weekly adult trap catches. Since a pheromone-baited trap can effectively catch insects within a range of two hectares (Tingle and Mitchell 1979), the spatial unit for the definition of the adult population abundance was considered two hectares in the present study. Our model required the estimation of the larval carrying capacity \({ }\gamma\). Considering the whole calibration dataset, we first calculated the average maximum observed adult abundance (284 individuals per trap per week). Based on this result we assumed a conservative value representing the carrying capacity of the adults in the spatial unit \(N_{K}^{4} = 320\). The relation between the seasonal fluctuations of adults (captured using pheromone-baited traps) and larvae (captured using sweep nets) of S. frugiperda was investigated for three consecutive years (1981–1983) by Silvain and Ti-A-Hing (1985). From their work, we extracted 10 datasets and calculated the average amount of larvae produced by a single adult (i.e. the ratio between larval and adult abundance at the peaks of the population) \(P = 9.34\). Based on this result, we calculated the carrying capacity of larvae \(\gamma = N_{K}^{4} P = 2989\) which was rounded to \(\gamma = 3000\) in the present study.

Data on species physiology

The development \(\nu^{i} \left( {T\left( t \right)} \right)\), mortality \(m^{i} \left( t \right)\) and fecundity \(F^{1} \left( t \right)\) rate functions were estimated considering data available in the literature on stage-specific responses of S. frugiperda exposed to different constant temperature conditions. Data referring to the average stage-specific duration in days were used for estimating the development rate function \(\nu^{i} \left( {T\left( t \right)} \right)\) (Barfield et al. 1978; Simmons 1993; Oeh et al. 2001; Busato et al. 2005; Milano et al. 2008; Barros et al. 2010; Ríos-Díez and Saldamando-Benjumea 2011; Garcia et al. 2018). Data referring to the stage-specific percentage survival were used for estimating the temperature-dependent component \(\mu^{i} \left( {T\left( t \right)} \right)\) of the mortality rate function \(m^{i} \left( t \right)\) (Barfield et al. 1978; Pashley et al. 1995; Murúa and Virla 2004; Busato et al. 2005; Milano et al. 2008; Barros et al. 2010; Garcia et al. 2018). Data referring to the temperature-dependent average total fecundity, average daily fecundity, and average duration in days of the oviposition period were used for estimating the temperature- \(g\left( {T\left( t \right)} \right)\) and the physiological age-dependent \(h\left( x \right)\) components of the fecundity rate function \(F^{1} \left( t \right)\) (Barfield et al. 1978; Pashley et al. 1995; Oeh et al. 2001; Milano et al. 2008; Barros et al. 2010; Garcia et al. 2018).

Temperature data

Yearly temperature data used as inputs during model calibration refer to the 5th generation of European ReAnalysis (ERA5-Land), reporting hourly air temperature data at a 0.1° × 0.1° spatial resolution (Muñoz Sabater 2019). Bilinear interpolation was used to obtain temperature data for each location of the calibration dataset. The current climatic scenario used to respond to the assessment questions was extracted from the Coordinated Regional Downscaling Experiment (CORDEX) (Jacob et al. 2014) and refers to the Coupled Model Intercomparison Project Phase 5 (CMIP5). The scenario is based on Representative Concentration Pathways (RCPs) which consider the greenhouse gases emissions up to the year 2100 (van Vuuren et al. 2011). The climatic scenario provides tri-hourly temperature data on 0.11° × 0.11° spatial resolution for the European domain over a period ranging between 2016 and 2025. Temperature data were regridded through bilinear interpolation to a regular 0.1° × 0.1° grid using Climate Data Operators command lines (Schulzweida 2019). We then averaged tri-hourly data over the whole decade (2016–2025) of the scenario to obtain an annual average temperature profile, which was assumed as the current climate (see Section S2 of supplementary material 1).

Results

Below are presented the answers to the four assessment questions, based on the results of the model.

AQ 1— Predicted population dynamics and limits of establishment of Spodoptera frugiperda in areas of current distribution

The graphical results of the model implemented along a south–north latitudinal gradient under the median-case assessment scenario are presented in Fig. 1. In the area of Miami Dade (Florida), the model predicted seven peaks (i.e. generations) per year; the predicted yearly average number of moths per trap per week was around 64 individuals and the maximum adult population abundance was around 165 individuals reached on the 6th generation. In the area of Alachua (Florida), the model predicted two generations per year; the yearly average adult abundance was around 17 individuals, and the maximum adult population abundance was around 98 individuals reached on the 2nd generation. Adult population abundance reached values lower than the adult population threshold \(A_{0}\) over the simulation period in the area of Tift (Georgia). Thus, the potential establishment of the pest was considered not possible in the above-mentioned area.

Fig. 1
figure 1

Observed (red asterisks) and simulated (blue line) adult population dynamics of Spodoptera frugiperda within the area of current distribution in North America: (left) Southern Florida, Miami Dade County, and (right) Northern Florida, Alachua County

AQ 2— Risk of establishment and potential distribution of Spodoptera frugiperda in Europe

Figures 2 and 3 show the risks of establishment of S. frugiperda in Europe under the three assessment scenarios. In the median-case scenario, risk of establishment was predicted in the southern coastal areas of the Mediterranean basin (Cyprus, Syria, Lebanon, Southern Turkey, Southern Italy, Southern, and Western Spain). A lower risk of establishment was expected in the Atlantic coasts of Portugal, and sporadic locations on the west coast of Sardinia. In the median-case assessment scenario, the measured area of potential establishment was 0.26% of the whole area under assessment. The area of establishment decreased by 89% (0.03% of the total assessed area) in the best-case assessment scenario and increased by 116% (0.57% of the total assessed area) in the worst-case assessment scenario. The northernmost latitudinal limit marking the presence of S. frugiperda populations was the 38th parallel north (Eastern Spain), the 43rd parallel north (Southern France), and the 44th parallel north (Northern Italy) in the best-case, median-case, and worst-case assessment scenarios, respectively.

Fig. 2
figure 2

Heat map showing the predicted distribution and the yearly average abundance of adult individuals of Spodoptera frugiperda per trap per week under the median-case assessment scenario

Fig. 3
figure 3

Heat maps showing the predicted distribution and the yearly average abundance of adult individuals of Spodoptera frugiperda per trap per week under the best-case A and worst-case B assessment scenarios

AQ 3— Predicted population dynamics of Spodoptera frugiperda in areas of potential distribution in Europe

Estimated population abundance within the area of potential establishment in Europe was highly variable depending on the assessment scenarios. The predicted yearly average number of moths per trap per week (± standard deviation) in the spatial unit was 5 (± 4) in the best-case, 17 (± 5) in the median-case, and 139 (± 22) in the worst-case assessment scenario. More details on the yearly population dynamic patterns of S. frugiperda are provided by the results of the local implementation of the model in areas with different suitability for the species in Europe (Fig. 4). The results of the model implemented in a highly suitable area (Cyprus) showed low population abundances at the beginning of the year due to low temperatures. Approaching the spring season, a rise in the adult population abundance was predicted according to temperature increase. Four adult population peaks (i.e. generations) were predicted around the 186th, 225th, 264th, and 314th DOY with the maximum adult population abundance reached on the third generation. Predicted adult population abundances during the peaks ranged between 90 and 130 individuals. After the fourth generation, a decline in the abundance of adults was observed, due to temperature drops in fall. The model implemented in the northernmost edge of the predicted establishment area in Europe (La Seyne-sur-Mer, Southern France) resulted in a single generation at around the 250th DOY. Adult abundance reached 90 individuals and then a sharp drop in population abundance was observed due to cold temperatures in fall. The model implemented in the Atlantic coast of Portugal (Alcochete) showed a single generation at around the 297 DOY with low adult population abundance during the peak (18 individuals).

Fig. 4
figure 4

Simulated (blue line) population dynamics of adults of Spodoptera frugiperda along a south–north latitudinal gradient within the area of potential establishment in Europe, in A Cyprus, B Southern France, and C Atlantic coast of Portugal

AQ 4— Risks linked to transient populations of Spodoptera frugiperda in Europe

Simulations of the inoculum in different periods of the year (considering the median-case assessment scenario) clearly showed that the species might be able to establish transient populations outside the predicted establishment area in Europe (Fig. 5). The results of the model implemented in Rădoiești (Southern Romania, 44th parallel north) and Ghedi (Northern Italy, 45th parallel north) showed risks linked to transient populations in all three introduction periods. A single generation was predicted for introductions at the 90th DOY, and at the 210th DOY and two generations were expected for introductions at the 150th DOY when weather conditions can be particularly suitable for the species. The predicted average number of moths per trap per week ranged between 19 and 20 adults with peaks at around 70 individuals (introduction at the 90th DOY) to 43–52 individuals (introductions at the 150th and 210th DOY). The model was implemented in areas further north in Europe (48th parallel north) in Ouarville (Northern France) and Engelsberg (Southern Germany). Introductions at the 90th DOY did not represent a risk of transient populations (\(N^{4} < A_{0}\) during the simulation period) due to the unsuitable environmental conditions affecting species’ survival. Introductions occurring during warmer periods in late spring (150th DOY) and summer (210th DOY) allowed the species to originate transient populations. However, only low yearly average adult population abundances (1–2 individuals) were predicted, thus representing a low risk linked to transient populations.

Fig. 5
figure 5

Population dynamics of adults of Spodoptera frugiperda simulating an inoculum of 5 adults at the 90th, (green line), 150th (blue line), and 210th (red line) day of the year in (upper left corner) Rădoiești (Southern Romania), (upper right corner) Ghedi (Southern Italy), (lower left corner) Ouarville (Northern France) and (lower right corner) Engelsberg (Germany)

Discussion

The model presented was able to satisfactorily predict the population dynamics, the variability in the number of generations, and the limits in the area of establishment of S. frugiperda along a latitudinal gradient within the area of current distribution in North America. The model implemented under the median-case assessment scenario predicted up to seven generations per year in an area where the species is well established (Florida, Miami Dade County) and only two generations per year in a location situated at the northernmost edge of the establishment for the species (Florida, Alachua County). These results are in agreement with the available observations reporting a high population abundance and population dynamics characterised by continuous generations throughout the year in tropical areas of Central America characterised by warmer temperature conditions (Sparks 1979; Busato et al. 2005), and around six generations in warm areas of North America (Luginbill 1928). According to the results of the model, the number of generations progressively decreases moving towards the northern areas of distribution of the species (Johnson 1987; Ramirez-Cabral et al. 2017; Schlemmer 2018). Correctly, the model predicted no establishment in areas (Tift County, Georgia) considered reached only by migratory populations. These results highlight the prominent role of climate in influencing the distribution and the dynamics of S. frugiperda (Capinera 2002; Garcia et al. 2018, 2019).

The predicted distribution and dynamics of S. frugiperda in Europe clearly highlighted risks of establishment of the species, especially in the coastal areas of the Mediterranean basin due to more favourable climatic conditions. In particular, higher average adult population abundances were predicted in the coastal areas of Southern Spain, Southern Italy, Greece, Cyprus, Southern Turkey, and Lebanon. Our results are in partial disagreement with the work of Du Plessis et al. (2018), which reported a low risk of establishment, but mainly risk linked to transient populations, associated with the pest in Europe. In agreement with our results, suitable areas for the establishment of S. frugiperda in the Mediterranean coasts of Europe have been reported in EFSA PLH Panel et al. (2018a), Early et al. (2018), Liu et al. (2020), Baloch et al. (2020), Zacarias (2020), and Tepa-Yotto et al. (2021). The predicted northernmost limit that might be reached by S. frugiperda was between the 38th and the 44th parallel north based on the different scenarios under assessment. These results are in partial disagreement with Liu et al. (2020), Baloch et al. (2020), Fan et al. (2020), and Tepa-Yotto et al. (2021). These authors reported areas potentially suitable for the establishment of the species further north, reaching the United Kingdom and Southern Sweden (although with low habitat suitability indices). Current knowledge of the biology of S. frugiperda seems to justify our predictions, especially in the light of the prominent role of climate in shaping the area of distribution of the species (Early et al. 2018). In particular, it is reported that the species does not enter diapause and suffers from cold weather conditions (Capinera 2002; Nagoshi et al. 2012). This might prevent the establishment of S. frugiperda in cold areas (EFSA PLH Panel et al. 2018a).

The population dynamics pattern of S. frugiperda predicted by the model in a highly suitable location in Europe showed that a rise in the adult population abundance (up to 90–130 adult individuals per trap per week) may occur during the early summer period, with up to four generations per year. A single generation per year and lower population abundances (around 18–90 adult individuals per trap per week) were predicted in less suitable locations in Europe. This result is in agreement with EFSA PLH Panel et al. (2018a) which reported up to four generations per year in the most suitable areas in Southern Europe. The population pressure expected in a suitable European location might represent a risk for local crop production.

Given the high migratory ability of the species, transient populations might represent a threat in areas outside the area of potential establishment of the species. The results of the model showed that in Europe, high risk due to transient populations can be expected in areas up to the 45th parallel north, with adult population abundances (20–70 adult individuals per trap per week) that might cause impacts on local crop production. Lower risks due to transient populations in Europe can be expected in areas up to the 48th parallel north, where unsuitable climatic conditions hinder the survival of the inoculum.

The physiologically-based modelling approach used in the present study requires the definition of biologically meaningful parameters describing the life-history of S. frugiperda. Parameters related to the temperature-dependent components influencing development, mortality, and fecundity were easily estimated using the large amount of data available in the literature. Conversely, the lack of data linked to density-dependent effects on adult fecundity and larval mortality forced us to make reasonable assumptions in the mathematical description of these components. Parameters representing density-dependent and biotic regulation affecting larval mortality have been calibrated using time-series population dynamics data. Parameters might be further fine-tuned if more data become available. In our model, temperature is the only abiotic variable influencing biological processes. If relevant, the model might be easily extended to include the influence of other environmental variables, such as relative humidity. The model does not consider the influence of different host plant species on the life-history (Chen et al. 2022), the dynamics, and the distribution of S. frugiperda (Baloch et al. 2020). Another source of uncertainty not considered in the model is represented by variability in the physiological responses associated with different strains of the species that might reach Europe (Sarr et al. 2021). The scenario-based approach we implemented seeks to cover part of the issues linked to parameters’ estimates and model’s limitations.

Conclusion

In this work, we present the results of a physiologically-based model applied to S. frugiperda to (i) predict the species’ population dynamics and abundance, (ii) assess the risk of establishment of the species in Europe, and (iii) predict the risk linked to transient populations in Europe. To the best of our knowledge, this is the first physiologically-based model simulating the life-history of S. frugiperda used for investigating the potential dynamics and distribution of the species in Europe. The physiologically-based modelling approach allowed us to simulate the influence of biotic (density-dependent effects and mortality due to biotic agents) and abiotic (e.g. temperature) variables on the life-history strategies of a pest (Soberon and Nakamura 2009; Gutierrez and Ponti 2013). This approach provides realistic predictions that are independent of data on the current distribution of the species that might be incomplete and/or biased (Wiens et al. 2009).

The model presented can provide fundamental elements for supporting the management of the pest considering different spatio-temporal scales and management contexts (Sperandio 2021). In case S. frugiperda is still absent from a territory, preventive measures should be taken, namely pest risk analysis (PRA), update of phytosanitary regulations (including the potential for response measures), inspection and diagnostics, and surveillance (EU 2018; FAO 2021). The quantitative outputs of the model (e.g. average population abundance, population dynamics, and potential distribution of the pest) provide fundamental information for the analysis of the risks posed by S. frugiperda in Europe. Risk maps generated by the model can be used to guide the implementation of detection surveys in the identification of high-risk areas that might be particularly suitable for the establishment of the pest. Similarly, the definition of the frequency and the intensity of the inspection measures can be guided by coupling risk maps on the potential establishment of the pest and information on trade routes and movement of people (EFSA PLH Panel et al. 2018b; FAO 2021). The model also provides relevant information on the potential impacts caused by transient populations that might represent a risk for local crop production. In case the species becomes established in mainland Europe, predictions on the population phenology and dynamics of S. frugiperda can be used for the timely implementation of control actions aimed at reducing pest population pressure and thus reducing the impacts on local crops (Rossi et al. 2019).

Authors’ contribution

GG and GS conceptualised the work. GG, GS, AS, and PG performed simulations. GG, GS, AS, and PG, interpreted outputs. GG, GS, and MC acquired and interpreted data. All authors drafted the manuscript.