Non-linear physiological responses to climate change: the case of Ceratitis capitata distribution and abundance in Europe

Understanding how climate change might influence the distribution and abundance of crop pests is fundamental for the development and the implementation of pest management strategies. Here we present and apply a modelling framework assessing the non-linear physiological responses of the life-history strategies of the Mediterranean fruit fly (Ceratitis capitata, Wiedemann) to temperature. The model is used to explore how climate change might influence the distribution and abundance of this pest in Europe. We estimated the change in the distribution, abundance and activity of this species under current (year 2020) and future (years 2030 and 2050) climatic scenarios. The effects of climate change on the distribution, abundance and activity of C. capitata are heterogeneous both in time and in space. A northward expansion of the species, an increase in the altitudinal limit marking the presence of the species, and an overall increase in population abundance is expected in areas that might become more suitable under a changing climate. On the contrary, stable or reduced population abundances can be expected in areas where climate change leads to equally suitable or less suitable conditions. This heterogeneity reflects the contribution of both spatial variability in the predicted climatic patterns and non-linearity in the responses of the species’ life-history strategies to temperature.

(year 2020) and future (years 2030 and 2050) climatic scenarios. The effects of climate change on the distribution, abundance and activity of C. capitata are heterogeneous both in time and in space. A northward expansion of the species, an increase in the altitudinal limit marking the presence of the species, and an overall increase in population abundance is expected in areas that might become more suitable under a changing climate. On the contrary, stable or reduced population abundances can be expected in areas where climate change leads to equally suitable or less suitable conditions. This heterogeneity reflects the contribution of both spatial variability in the predicted climatic patterns and non-linearity in the responses of the species' life-history strategies to temperature.

Introduction
Invasion history and impacts of Ceratitis capitata in Europe The Mediterranean fruit fly (or medfly), Ceratitis capitata (Wiedemann) (Diptera: Teprhitidae) is a phytophagous insect pest considered one of the most important threats for horticulture and the fruit industry (Liquido et al. 1990;Papadopoulos et al. 2001;De Meyer et al. 2002;Weldon et al. 2018). The medfly has been reported on more than 300 host plants including citrus, stone fruits, pome fruits, tomatoes, figs and others (Morales et al. 2004;Meats and Smallridge 2007). The damage on fruit is caused by adult oviposition, which also facilitates the transmission of bacteria and fungi, and by larval trophic activity (Cayol et al. 1994). The need to apply pre-and postharvest control measures and limitations to export products from infested areas to medfly-free areas results in economic impacts (Hulme 2009). The medfly has been classified as an EPPO A2 quarantine pest (Rossler and Chen 1994). The species originated from eastern Sub-Saharan Africa (White and Elson-Harris 1992;Gasperi et al. 2002;Malacrida et al. 2007). Facilitated by human-mediated transportation and trade (De Meyer et al. 2002;Bonizzoni et al. 2004;Malacrida et al. 2007), the medfly colonised the whole African continent, the Mediterranean Basin (Robinson and Hooper 1989), the Middle East, western Australia (Bonizzoni et al. 2004) and Latin America (Harris and Lee 1989;Headrick and Goeden 1996;Vera et al. 2002).

Ecology of Ceratitis capitata
The biological traits of C. capitata explaining its successful invasion and establishment in new areas are related to its high polyphagy (Malacrida et al. 2007), short generation times (Duyck and Quilici 2002;Grout and Stoltz 2007), capacity for long distance flights (Meats and Smallridge 2007), and adaptability to different climatic conditions (Nyamukondiwa et al. 2010). Responses of C. capitata to climate are of particular interest. Recent research demonstrated the capacity of C. capitata to tolerate a wide range of temperature conditions both in laboratory (Nyamukondiwa et al. 2010;Weldon et al. 2011Weldon et al. , 2018 and field conditions (Nyamukondiwa et al. 2013), explaining its success in colonising temperate areas (Malacrida et al. 1992;De Meyer et al. 2002). Furthermore, transient populations might survive in refuge areas, allowing C. capitata to reinvade new regions as soon as the seasonal weather conditions become suitable for the species (Carey 1991). The climatic conditions that limit the establishment of C. capitata are still under debate. The species was traditionally considered limited to the 41st parallel north (Robinson and Hooper 1989;Israely et al. 2004) and to habitats where temperatures drop below 10°C (Mwatawala et al. 2015). However, Papadopoulos et al. (1996Papadopoulos et al. ( , 1998Papadopoulos et al. ( , 2001 reported the capacity of the species to overwinter as larva in cold winter conditions and survive in subfreezing temperatures. Stable C. capitata populations have been reported in southern France (Cayol and Causse 1993), northern Italy (Rigamonti et al. 2002;Rigamonti 2004Rigamonti , 2005Zanoni 2018;Zanoni et al. 2019), and Austria (Egartner et al. 2017) suggesting the capacity of the species to complete its life-cycle in areas far above its traditional northward distribution limit. Climate change is expected to play a primary role in altering the potential distribution and performance of the species worldwide (Gutierrez and Ponti 2011;Sultana et al. 2020), and in the species' capacity to reinvade new areas by local transient populations (Lux 2018).
Requirements for the development of modelling tools supporting the management of Ceratitis capitata Given the severe impacts to agriculture and the high invasive capacity of C. capitata, it is fundamental to produce reliable risk scenarios that account for an appropriate spatial and temporal resolution in relation to the scale of management. Models are suitable tools for exploring how the main drivers of risk (i.e., pest distribution and abundance) (EFSA PLH Panel 2018) are affected by climate variability and change both in time and space (Hill et al. 2016b).
Several models estimating the potential distribution and the risks linked to C. capitata have been proposed.
The vast majority are species distribution models (SDMs), which estimate the nature of the relationship between the species' current distribution (i.e., presence/absence) and environmental variables (i.e., temperature, rainfall, land use etc.) (Vera et al. 2002;Gevrey and Worner 2006;De Meyer et al. 2007;Li et al. 2009;Szyniszewska and Tatem 2014;Godefroid et al. 2015;Kaya et al. 2017). These models are powerful tools that, under specific assumptions (Soberon and Nakamura 2009;Wiens et al. 2009;Warren 2012), allow estimating the current and the projected distribution of a species in the form of habitat suitability indices. In particular, these models are useful for estimating the potential distribution of pest species in the case that life-history strategy parameters at the individual level are not known and only occurrence data are available (Srivastava et al. 2021). One of the main shortcomings of SDMs is that they cannot provide information on potential population abundance, which is one of the most important elements for assessing the impacts of a pest in a certain area (Gilioli et al. 2017a;EFSA PLH Panel 2018). Additionally, these models lack of a mechanistic description of the species' life-history strategies (Soberon and Nakamura 2009) not allowing a description of the non-linear physiological responses to environmental drivers at the individual level and the consequences on species' population dynamics (Gutierrez and Ponti 2013;Gilioli et al. 2014).
When data on individual physiological responses to temperature and other environmental variables (e.g. humidity, rainfall, land use) are available, processbased models allow a fine description of the biological processes at the basis of population dynamics, both in space and time (Gutierrez and Ponti 2011;Régnière et al. 2012b;Gilioli et al. 2016;Maino et al. 2016). Despite the intensive efforts required for the development of process-based models (Kriticos et al. 2012) the biological realism at the basis of these models may allow obtaining more accurate predictions of the potential distribution and abundance (Gutierrez and Ponti 2013;Gilioli et al. 2014;Pasquali et al. 2015;Ponti et al. 2015).
Here we propose a process-based model to investigate the potential role of climate on the life-history strategies of C. capitata, and how physiological responses at the individual level influence the distribution, the abundance and the activity of the pest in Europe.
Many potential drivers (e.g. humidity, rainfall, land use etc.) can influence the population dynamics and thus the potential establishment of a pest species. Given the major role of temperature on poikilotherm organisms (Logan et al. 2003;David Logan et al. 2005;Tamburini et al. 2013;Sultana et al. 2017), we used climatic variables to predict the potential establishment, abundance and activity of C. capitata. The climatic scenarios used in the present study describe the local temperature patterns at a high spatiotemporal resolution, allowing for the estimation of the role of temperature on species' individual physiology and population dynamics. This enables the exploration of current and future scenarios linked to pest distribution and impacts (Kriticos et al. 2015;EFSA PLH Panel 2018). In addition to temperature, we included the role of density-dependent factors influencing larval and adult population abundance due to the effects of intra-specific competition, as well as a mortality component that accounts for an extrinsic mortality factor combining the effects of inter-specific competition and predation due to natural enemies (Agnew et al. 2002;Tamburini et al. 2013). Finally, we ensured a complete procedure for model calibration and validation, which is fundamental to guarantee reliable pest risk scenarios (Venette et al. 2010;Venette 2015). The model presented provides population projections in terms of stage-specific pest distribution, abundance and activity that are essential elements for assessing the risks linked to a pest. Results consider in particular larvae and adults because they directly affect crop productivity and quality. We compare the obtained simulation results with data and models available in literature. Finally, we thoroughly discuss the effects of non-linear responses to climate change on the current and the future distribution, abundance and activity of C. capitata.

Data on population abundance and distribution
Population abundance data used for model calibration refer to time series adult trap catch data for C. capitata collected from five independent studies (Fig. 1).
Studies were conducted on fruit orchards or on wild plants in Latium (Italy) in 2007 (Sciarretta and Trematerra 2011), in Tarragona (Spain) in 2007(Martínez-Ferrer et al. 2010, in Thessaloniki (Greece) in 1991 (Papadopoulos et al. 2001), in Dubrovnik (Croatia) in 2002 (Bjeliš et al. 2007), and in Kibbutz Zova (Israel) between 1994(Israely et al. 2004. Information on the occurrence of the species in Trentino Alto Adige/Südtirol (the northern distribution limit of C. capitata in Italy) is used for validating the model's capacity to predict species' establishment along an altitudinal gradient (see Sect. Model validation). In the area of Riva del Garda, at the foothills of the Alps, the species has been sporadically reported since 1990 with only marginal impacts on fruit orchards. However, in 2016, the area was characterized by a high presence of the species showing a transition towards more stable populations (Zanoni 2018). In two other locations in the same region, Bressanone and Bolzano, further up north from Riva del Garda and into the Alps, the species occurrence was not reported, therefore these areas can be considered not suitable for its establishment.
Occurrence (i.e., presence) data collected from the Belgian Biodiversity Platform (http://www. biodiversity.be) and from the Global Biodiversity Information Facility dataset (https://www.gbif.org/) were used to test model's ability to predict the current distribution of C. capitata in Europe (see Sect. Model validation). The occurrence of the species does not necessarily mean that established populations are present in a certain location. Therefore, we selected occurrence data only for those countries where the species was considered present, based on the information reported in the European and Mediterranean Plant Protection Organization (EPPO) Global database (https://gd.eppo.int/).

Point-based temperature data
Temperature data used for model calibration and for model validation based on predicting the altitudinal limits of the species were obtained from the Global Surface Summary of Day Product created by the US National Climatic Data Centre (https://www.ncdc. noaa.gov/). Daily minimum and maximum air temperatures were extracted from the nearest station to each location representing calibration or validation data. Hourly temperature data were then calculated from daily maximum and minimum air temperature using the algorithm described in Woodhead (1979) and applied in Gilioli et al. (2014).

Climate scenarios
Climate scenarios were obtained from the Coordinated Regional Downscaling Experiment (CORDEX) and refer to Regional Climate Model (RCM) simulations for the European Domain (from 27th to 72th parallel north and from the 22nd meridian west to -45th meridian east) (Jacob et al. 2014). Scenarios provide tri-hourly bias-adjusted air temperature at a 0.11°9 0.11°spatial resolution. The scenarios refer to the Coupled Model Intercomparison Project Phase 5 (CMIP5) and are based on Representative Concentration Pathways (RCPs) (Van Vuuren et al. 2011), i.e., the greenhouse gases emission scenarios corresponding to a determined radiative forcing up to the year 2100. We used two RCPs available for the European domain: the RCP4.5 (referring to a stabilization of radiative forcing from 2150 onwards at 4.5 W/m 2 ) and the RCP8.5 (referring to rising radiative forcing crossing 8.5 W/m 2 at the end of twenty-first century). Temperature data were re-gridded to a regular 0.1°9 0.1°grid through bilinear interpolation using the command-line software Climate Data Operators (Schulzweida 2019) and averaged over a period of 10 years in order to avoid extreme climatic variations. The climatic scenarios used in the present study refer to 2020 RCP4.5 (average for the years 2016-2025) assumed as the current climate (see supplementary material 1 Table S1 and Figure S1), 2030 RCP8.5 (average for the years 2026-2035) and 2050 RCP8.5 (average for the years 2046-2055).

Model assumptions
The population dynamics model for C. capitata is based on the following assumptions.
Individual life history strategies are represented by development, mortality and fecundity stage-specific rate functions. These functions depend on temperature, which is the main environmental driving force influencing poikilotherm population dynamics (Gutierrez 1996;Régnière et al. 2012a;Rebaudo and Rabhi 2018). Individual physiological responses to temperature are non-linear.
Besides temperature, the fecundity rate function includes the effects of aging on the number of eggs produced by each adult female. This allows simulating a pre-oviposition period (i.e., no eggs production), a peak in the production of eggs, and then a reduction of fecundity towards the end of the life of the adult female.
In addition to temperature, the mortality rate function has two more components. We include a mortality term to account for density-dependent responses simulating the intra-specific competition, and a mortality term accounting for an extrinsic control due to inter-specific interactions (with competitors and predators) both acting on the trophic stages (i.e., larvae and adults) (Carey et al. 1995;Dukas et al. 2001;Diamantidis et al. 2020).
The model simulates the local population dynamics in each node of a regular grid 0.1°9 0.1°(lattice model) covering the European continent. Population dynamics in a node are influenced only by conditions in that specific node (temperature, population abundance and structure). It is assumed that nodes are isolated from one another (no flux of individuals among the nodes) and identical in terms of resources for the pest and biotic control agents.
We denote by x 2 ½0; 1 the individual physiological age, defined as the level of development of an individual within a stage over time Pasquali 2007, 2010). With x ¼ 0 we denote individuals at the beginning of the development within a certain stage, while x ¼ 1 refers to individuals that have completed the development within that stage. Let / i t; x ð Þ be the population abundance of individuals in stage i at time t, with physiological age in ½x; x þ dx, where dx stands for an infinitesimal age increment. Hence, the number of individuals in stage i at time t is We described the population dynamics of C. capitata through a process-based model based on the Kolmogorov equation (Lee et al. 1976;Weiss 1986;Bergh and Getz 1988;Iannelli 1995;Buffoni and Pasquali 2007;Rafikov et al. 2008;Solari and Natiello 2014;Lanzarone et al. 2017). This modelling approach has already been applied to investigate the dynamics and the potential distribution of Bemisia tabaci (Gilioli et al. 2014), Lobesia botrana (Gilioli et al. 2016), Pomacea caniculata (Gilioli et al. 2017b, c), Aedes albopictus (Pasquali et al. 2020) and the phenology of Cydia pomonella (Pasquali et al. 2019). Full mathematical details of the model presented can be found in supplementary material 2.

Model rate functions estimation
The development and the mortality rate functions were defined for each developmental stage. The fertility rate function was defined for the adult stage. The parameters of the rate functions were estimated through the least squares method using the lscurvefit function in MATLAB version R2018a (MATLAB, R2018a, The MathWorks, Inc., MA, USA).Termination tolerance was set on the first-order optimality at 10 -6 . In the following, we will omit for brevity the dependence on the stage i ¼ 1; . . .; 4, specifying in the corresponding tables the stage-specific parameters of each function.

Development rate function
The stage-specific temperature-dependent development rate function m T ð Þ was described by the Brière function (Briere et al. 1999 where T ¼ TðtÞ is the air temperature at time t, r is a scaling parameter, T inf and T sup represent the stagespecific minimum and maximum temperature at which the development occurs, respectively. The parameters of function (1) were estimated from experimental data on temperature-dependent development times (in days) of C. capitata exposed to different constant temperatures. The parameters of function (1) were estimated using the least squares method based on data from Carey (2011), Diamantidis et al. (2011), Duyck and Quilici (2002), Grout and Stoltz (2007), Gutierrez and Ponti (2011), Muñiz and Gil (1986), Ricalde et al. (2012), Rössler (1975), Shoukry and Hafez (1979) and Vargas et al. (1997Vargas et al. ( , 2000. Parameters of the development rate functions and related figures for each stage are reported in supplementary material 3 (Table S2 and Figure S2).

Mortality rate function
The stage-specific temperature-dependent finite mortality rate M(T) (percentage of mortality in a stage at a given temperature T) was expressed by where T l inf and T l sup are the abscissa of the intersections among the polynomial function defined in (2) and the constant function k ¼ 0:9 representing the hypothesised maximum stage-specific mortality. The parameters of function (2) were estimated from experimental data on C. capitata survival at different temperatures using the least squares method based on data from Diamantidis et al. (2011), Duyck and Quilici (2002), Gutierrez and Ponti (2011) and Vargas et al. (2000). As no suitable data were available on the stagespecific mortality of the adult stage, we assumed that it was equal to the stage-specific proportional mortality of pupae. The parameters of function (2) and related figures for each stage are reported in supplementary material 3 (Table S3 and Figure S3).
Following the method proposed in Gilioli et al. (2016) and Pasquali et al. (2020) we derived the temperature-dependent instantaneous mortality rate l T ð Þ from the finite mortality rate in (2) where the parameters c Lj and c Rj , j ¼ 1; 2; 3 of the outer branches of l were inferred to obtain the desired slope of the curve and a sufficiently regular connection with the middle branch (l of class C 1 ). The stagespecific parameters of the functions (3) and related figures for each stage are reported in supplementary material 3 (Table S4 and Figure S4). For the egg and pupal stages (i ¼ 1; 3) the stagespecific mortality rate m i T ð Þ appearing in supplementary material 2 was defined as For the larval and the adult stages (i ¼ 2; 4) we also considered a density-dependent mortality term (intraspecific competition) and a mortality term due to interspecific interactions (competition and predation). For these two stages the mortality rate was modified as follows where N i max refers to the maximum stage-specific abundance that might be sustained with respect to the available trophic resources, while the term d i is an extrinsic biotic mortality due to the effects of competitors and natural enemies. The parameters in (4) were estimated from population dynamics data following the procedure presented in Sect. Model calibration.

Model calibration
The model calibration procedure consisted of estimating the two biotic mortality terms in the mortality rate function (4) applied to larvae and adults. Parameters of the biotic mortality terms in function (4) were obtained through the minimisation of the squared distance between the estimated and the observed population dynamics of C. capitata adults (least square method) considering time-series trap catches data collected in five different locations in the Mediterranean Basin (Italy, Spain, Portugal, Croatia and Israel). The calibration procedure allowed estimating parameters that represent population variability in the physiological responses under different environmental conditions. Yearly temperature datasets at hourly time resolution were used as input data for model calibration.
Denoting by N 4 j t i ; d 2 ; d 4 ; N 2 max ; N 4 max À Á the adult abundance at location j at time t i obtained by integrating the solution (see supplementary material 2) with respect to physiological age, for parameters d 2 ; d 4 ; N 2 max ; N 4 max (see (4)), we defined the function where A j t i ð Þ is the observed adult abundance at location j at time t i and R j is the number of available data for location j. Then, we found the optimal parameters d 2 ; d 4 ; N 2 max ; N 4 max by minimising Q, namely The minimization was performed using the MATLAB function fmincon, with the default tolerance of 10 -5 .

Model validation
Model validation was based on testing the model's capacity to predict i) the altitudinal limit of the species distribution in a transect crossing the Alps in northern Italy, and ii) the current distribution of C. capitata in Europe, taking into account the countries where the species is established. For the altitudinal limit, we tested the model in predicting establishment of the species in Riva del Garda, Bressanone and Bolzano (Trentino Alto Adige, Italy) in 2016 (Zanoni 2018) using local yearly temperature datasets at hourly time resolution. The model's ability to predict the current distribution of C. capitata was tested by comparing data on C. capitata occurrence in Europe with the predicted distribution of the species obtained by model simulation under current climate scenario.

Model simulations and climate scenarios
We used model simulations to estimate the current and the future distribution, abundance and activity of C. capitata. Simulation results were obtained by numerically solving the model in each point of the 0.1°9 0.1°grid covering most of Europe, but including some of Central Asia and North Africa (Fig. 1). The mathematical model, that is a partial differential equation, was discretised in both time and physiological age by standard schemes and the simulations were performed using the scientific software MATLAB.
The model outputs presented in this study are: (i) the yearly mean adult and larval abundance considering only the abundances higher or equal to stagespecific abundance thresholds, and (ii) the yearly adult and larval activity calculated as the number of days in which the adult and the larval abundance is higher or equal to stage-specific abundance thresholds. The stage-specific abundance thresholds were set to one individual for the adult stage and 30 individuals for the larval stage. The area of potential distribution of adults and larvae are defined by the points of the grid in which the yearly mean adult and larval abundance is higher or equal to the stage-specific abundance thresholds. For each node of the grid, the calculated yearly mean adult and larval population abundance refers to a sampling unit that corresponds to the range covered by a single trap for C. capitata adults baited with trimedlure. Manoukis et al. (2015) reported that a single trap effectively catches adult individuals approximately on a 7-14 m radius. Thus, in our work, we assumed that the spatial unit A 0 covered by each grid's node is equal to the area of influence of a single trap (7-14 m radius). For any area of size A, adult abundance varies linearly with the ratio A=A 0 . In our study, we assume a constant trap efficacy; therefore, to obtain the real population density estimates in terms of individuals per area unit, it is necessary to introduce a correction factor accounting for trap efficacy.
For all the simulations, the initial condition on the 1 st of January was set to five larvae uniformly distributed within their physiological age interval ½0; 1. In order to obtain stable population dynamics as well as model outputs independent from the initial conditions, we ran the model for ten consecutive years. Model outputs were calculated in the last year of the simulation. Model outputs obtained for the current climate scenario (2020) were then compared with those obtained for the 2030 and 2050 scenarios (Sect. Distribution, abundance and activity of C. capitata under climate change scenarios).

Model calibration and validation
The parameters obtained throughout the calibration procedure representing the biotic components of the mortality function (4) d i and N i max for larvae and adults were used in the present study (Table 1).
The validation procedure reveals that the model is able to interpret the altitudinal limit of the distribution in a transect crossing the Alps in northern Italy (see supplementary material 4 Figure S6). The model successfully predicts the presence of stable populations in Riva del Garda (Trento province, Italy) (Zanoni 2018). A sharp drop to zero in the population abundance is predicted along the altitudinal gradient in Bolzano and Bressanone (Bolzano province, Italy). Stable populations of C. capitata were not reported in the aforementioned locations, due to weather conditions that prevent the survival of the species in winter. The vast majority of the occurrence data for C. capitata falls within the area of distribution predicted by the model (see Fig. 1), confirming the model's capacity to accurately predict the overall distribution and the northernmost distribution limit of C. capitata in Europe.
Distribution, abundance and activity of C. capitata under current climate scenario The simulated distribution and the yearly mean abundance of C. capitata adult and larvae in Europe for 2020 is reported in Fig. 1. The predicted current distribution of C. capitata covers the Iberian Peninsula, Morocco, Algeria, Tunisia, Greece, part of Turkey, Syria, the coastal area of the Balkans, the Danube area of Romania and Bulgaria, Italy and southwestern France. Mountainous areas are not suitable for the species, as is clearly shown in Spain, Italy and the Balkans. Established populations of C. capitata can potentially reach the 48th parallel north in France with some sporadic populations in the central part of the country where climatic conditions are more suitable. In southern France and northern Italy, the distribution is limited to the 46th parallel north. This latitude also marks the distribution limit in eastern Europe where the distribution is not continuous over the Balkans, but it is influenced by local climatic conditions. The yearly average number of adults ranges between 1 (the threshold set to define the possibility of establishment) and 30 individuals per spatial unit (the maximum level of abundance reported in southern Italy). Yearly mean larval abundance ranges between 30 (the threshold set to define the possibility of establishment) and the maximum value of 2545 individuals per spatial unit reported in Morocco. In Europe, a higher average number of adults (30 individuals) and larvae (above 2500 individuals) are reported in Andalucía (southern Spain), Sicily (southern Italy) and Cyprus. Yearly mean population abundance decreases along a south-north gradient. Higher population abundance is predicted in the warmer areas of southern Europe and on the eastern coasts of the Mediterranean Basin (Turkey and the Middle East). Population abundance drops in the northernmost distribution limits of the species and towards mountainous areas.
In supplementary material 5 ( Figure S7) we report a graph showing the relation between the average yearly temperature and the average yearly adult and larval abundance under current (year 2020) climatic scenario. The graph highlights the non-linear response of population abundance to the yearly average temperature. The model predicts that populations of C. capitata are not able to establish in areas characterised by a yearly average temperature below 12.4°C. Above this threshold, population abundance clearly increases with temperature. The rate of increase progressively diminishes when average temperatures are higher than 15-16°C.
The simulated adult and larval activity of C. capitata in Europe for 2020 is presented in Fig. 2. Low activity is recorded towards the northern distribution limit of the species or in the inland areas, while moving south and along the Mediterranean coast, an increased activity of both adults and larvae is predicted. The model predicts a continuous activity of adults in the coastal and in the inland areas of Morocco, Algeria, and Tunisia, and the southern coastal areas of Spain, Portugal, Italy, Greece, Turkey, Cyprus and Middle East. The area where larval populations are present throughout the year is more extended, and covers the Mediterranean coasts of southern Europe, the inland areas of Spain, Portugal, Turkey, and the Middle East.
Distribution, abundance and activity of C. capitata under climate change scenarios We ran the model using the 2030 and 2050 climate scenarios and compared the results with the 2020 scenario. The results (Figs. 3 and 4) are presented in term of variation of the output variables (population abundance and activity of C. capitata) with respect to the 2020 scenario. A positive variation means an increase in the output variable (higher abundance or longer period of activity with respect to the 2020 scenario); a negative variation shows a decrease in the output variable (lower abundance or shorter period of activity with respect to the 2020 scenario).
Results of the comparison between years 2030 and 2020 are shown in Fig. 3. In 2030, the potential distribution of the species increases in the inland areas of France, and reaches southern Germany, Switzerland, Austria, and Hungary. A range of variation of -4 (-13%) and ? 14 (? 47%) for adult individuals and -234 (-9%) and ? 700 (? 27%) for larval individuals is predicted in 2030 compared to 2020. The variation of population activity ranges between -64 and ? 110 days for adults and between -152 and ? 220 days for larvae.
Results of the comparison between the year 2050 and 2020 are shown in Fig. 4. In 2050, the model predicts established populations reaching northern France, central Germany, Poland, and the Netherlands. The comparison between the 2020 and 2050 scenarios shows a variation of population abundance ranging between -4 (-13%) and ? 18 (? 60%) for adults and between -129 (-5%) and ? 888 (? 35%) for larvae. The variation of population activity ranges between -50 and ? 173 days for adults and -28 and ? 279 days for larvae.
Overall, scenario comparisons show both an increased population abundance and activity towards the species' northern distribution limit and towards the inner areas of the Mediterranean basin. The decrease in population abundance is particularly evident in continental Turkey, the Iberian Peninsula, Morocco,  Figures a2 and b2 show the activity (expressed as the number of days in which the adult population abundance is equal or greater than one individual and the larval population abundance is equal or greater than 30 individuals per spatial unit) of adults and larvae, respectively Tunisia and Algeria. The decrease in population activity is more pronounced in 2030 (especially, in central Spain and Portugal).
In Fig. 5, we represent the change in the average population abundance of adult individuals of C. capitata with respect to the average yearly current temperature and the temperature variation under climate change (yearly average variation between 2020 and 2050). The graph highlights the complex and non-linear population responses respect to climatic variations; the main observed patterns are discussed below.

Discussion
Predicted current distribution, abundance and activity of C. capitata Model simulations show that the current northernmost distribution limits of C. capitata reaches the 46 th parallel north in Italy and in eastern Europe while sporadic populations are able to reach the 48 th parallel north in France. The predicted current distribution of C. capitata complies with the most updated data related to the distribution of the species in Europe (http://www.biodiversity.be; https://www.gbif.org/). This serves as an indicator on the validity of the predictions of the model presented. In our model, areas outside the predicted current distribution limits should be considered as climatically unsuitable for the (expressed as the number of days in which the adult population abundance is equal or greater than one individual and the larval population abundance is equal or greater than 30 individuals per spatial unit) of adults and larvae, respectively establishment of C. capitata even though an accidental introduction of the species might result in the local presence of transient populations surviving in refuge areas or under particularly warm winter conditions. In contrast with other studies (Li et al. 2009; Szyniszewska and Tatem 2014) our model does not predict suitable habitats for the establishment of C. capitata in the UK, Ireland and in vast areas of northern France, Belgium, Netherlands and Poland. Our results seems more in line with current knowledge, attesting the northernmost presence of stable C. capitata populations in southern France, northern Italy and Austria (Cayol and Causse 1993;Rigamonti et al. 2002;Rigamonti 2004Rigamonti , 2005Egartner et al. 2017;Zanoni 2018;Zanoni et al. 2019). Other results similar to our predictions are found in Vera et al. (2002) and, limited to Italy, in Gutierrez and Ponti (2011).
In relation to population abundance and activity, overall the model predicts a south-north gradient with higher populations observed in southern Europe and along the Mediterranean coastal areas. Within the area of potential establishment, climate prevents the establishment of the species moving towards mountainous areas. This result is in line with those provided by other authors (Vera et al. 2002;Li et al. 2009;in Gutierrez and Ponti 2011;Szyniszewska and Tatem 2014). The continuous presence of adult individuals within the year is predicted for the southern Mediterranean coasts and the Middle East, while the area interested by the continuous presence of larvae is more extended and reaches the southern coasts of France. This result suggests a high capacity of C. capitata to overwinter as larva, especially in Mediterranean areas, while adults are not able to survive cold winter conditions (Papadopoulos et al. 1996;Peñarrubia-María et al. 2012). Even if this hypothesis is still under debate (Israely et al. 2004), results from other studies suggest that the species might be able to overwinter in areas characterised by temperate winters or subfreezing temperatures (Papadopoulos et al. 1996(Papadopoulos et al. , 1998(Papadopoulos et al. , 2001Rigamonti 2004;Escudero-Colomar et al. 2008).
Our predictions on the potential abundance and activity of the species should be carefully considered due to a series of limitations affecting our modelling approach. First, only the role of temperature on the demographic processes is taken into account in obtaining the predicted population abundance. Indeed, temperature represents one of the major driving force influencing species' responses at the individual and population levels (Gutierrez 1996;Régnière et al. 2012a;Rebaudo and Rabhi 2018). The realism of population abundance predictions might benefit from the inclusion of the effects of abiotic factors other than temperature as long as information on the biological responses to those factors are available. The modelling framework here proposed can be easily adapted to include the influence of other abiotic (e.g. relative humidity) or biotic (e.g. host plants phenology) drivers. The density-dependent control and the extrinsic mortality factor are introduced without considering local information on resource availability and abundance of natural enemies and competitors. Another limit of our approach is that model calibration was performed using adult trap catches that represents the standard protocol for monitoring C. capitata. Timeseries abundance data referring to the other life-stages (e.g. eggs, larvae and pupae) are rare in literature. The availability of high resolution and harmonised data on both adult and larval stages would allow finer model predictions on population density (i.e., individuals per area unit) for the larval stage, which is the most relevant for assessing crop impacts. Despite these limitations, the outputs predicted by our model still provide biologically meaningful information for assessing the pest population pressure in both space and time. This information is highly relevant for the quantitative estimation of the risk posed by this species to host plants (EFSA PLH Panel 2018).
Potential effects of climate change on distribution, abundance and activity of C. capitata The effects of climate change on model outputs, namely species distribution, abundance and activity, are not homogeneous in the whole territory under investigation and no simple trends or gradients can be derived (Walther et al. 2002;Svobodová et al. 2014;Battisti and Larsson 2015;Hill et al. 2016a). Model projection allows identifying patterns of variation in Europe that describe the joint effects of the spatial heterogeneity in the expected change in temperature and the non-linear temperature-dependent response of the species' life-history traits. The results of our scenario comparison show five different patterns in the investigated area.
A northward expansion of the area of potential distribution of the species. The area of distribution increases of ? 12% in 2030 and ? 42% in 2050 respect to the 2020 scenario. The northward expansion of the species is evident in central Europe. Therefore, areas that are currently not suitable for the establishment of the species due to low temperatures might become more suitable under a changing climate.
A rise in the altitudinal limit marking the presence of established populations. This is particularly evident along the border between Spain and France (Pyrenees), in Italy (Apennines) and along the borders between Italy, France, Switzerland, and Austria (Alps).
Increase in population abundance of C. capitata. The percent area predicted to have an increase in population abundance is 62% (adults) and 64% (larvae) in 2030, and 79% (adults) and 86% (larvae) in 2050. The increase in the abundance and activity is prominent towards the northward distribution limit and along the altitudinal limit of distribution of the species. This increase can be explained in terms of species' physiological responses to temperature variation. In these areas, temperature variations act to the almost linear and positive (slope [ 0) area of the rate functions expressing the species' life history traits (development, fecundity and survival).
Stability in the population abundance of C. capitata. The percent area predicted to have stable population abundance is 21% (adults) and 13% (larvae) in 2030, and 7% (adults) and 4% (larvae) in 2050. Stability of the populations is evident in southern Portugal, the coastal areas of Italy and Greece, limited areas of the Middle East, Morocco and Algeria. In these areas, temperature variations are likely within the range of optimal temperature for the rate functions. Thus, negative and positive effects on species' life history result in generally stable populations.
A decrease of population abundance of C. capitata. The percent area predicted to have a decrease in population abundance is 17% (adults) and 23% (larvae) in 2030, and 14% (adults) and 10% (larvae) in 2050. This pattern is particularly relevant in vast areas of Morocco, Algeria and Tunisia, the Middle East, Turkey, southern Spain and limited areas in southern Italy. This decrease in population performance is probably due to the effects of temperature increasing to the point where the rate functions are non-linear and negative (slope \ 0).

Non-linear population responses of C. capitata to climate change
The non-linearity of population responses to climate change derives from the interaction between individual's physiological responses and the spatial and temporal heterogeneity of current and future climate. The simulated changes in the average population abundance of adult individuals of C. capitata with respect to the average yearly current temperature and the temperature variation under climate change (yearly average variation between 2020 and 2050) allow the identification of the following.
Higher positive variation in the average adult population abundance is reported for the areas characterised by lower yearly average current temperature and where higher positive temperature variation is expected in 2050.
Adult population abundance variation shows a nonlinear decrease as the yearly average current temperature rises.
In the areas characterised by yearly average current temperature approximately above 17°C and a temperature variation below 1.5°C, the variation in the abundance of C. capitata might even be negative.
In the areas characterised by temperature higher than 22°C, and with positive temperature variations above 1.5°C due to climate change, we observe a positive and mild increase in the abundance of C. capitata.

Conclusions
In this paper, we present the results of a process-based, stage-structured model providing quantitative information on the distribution, abundance and activity of C. capitata under current and future climatic scenarios. The model presented provides future pest projections that can inform decision-makers in the allocation of effort and resources for pest management in the areas considered at major risk of invasion of C. capitata (Hill et al. 2016a;Weldon et al. 2018). The model allows for the investigation of non-linear physiological responses to environmental forcing drivers (and their variations both in time and in space) at the individual and population levels. The capacity to describe and interpret the heterogeneity in population responses to climate change are fundamental elements in guiding the assessment and the management of future risks linked to C. capitata at both the regional and the local spatial scales (Hill et al. 2016a;Weldon et al. 2018). The model can also be applied at the local scale for estimating high temporal and spatial resolution outputs in terms of within-year pest population dynamics and phenology (Gilioli et al. 2021). Additionally, this output can be used to facilitate the scheduling of monitoring or pest treatment activities at the local scale (Manoukis and Hoffman 2014;Kean and Stringer 2019;Rossi et al. 2019). Simulations at high spatial and temporal resolution are a valid support for pest management in areas where the species is well established, but they might also be applied to investigate transient populations of C. capitata and assess the risks of establishment of the species in new areas (Papadopoulos et al. 2013;Kean and Stringer 2019). this paper. We are grateful to Professor Andrea Sciarretta (University of Molise, Campobasso, Italy) for having provided relevant information on the biology and population dynamics of C. capitata and on the main trapping systems used to monitor the species. We kindly thank Dr. Fabrizio Ferrari (TerrAria srl, Milan, Italy), Dr. Simone Gabriele Parisi (Milan, Italy) and Professor Benjamin Zaitchik (Johns Hopkins University, Baltimore, USA) for the support provided in the elaboration of the climatic scenarios used in the present study. We kindly thank Dr. Matilde Zubani (University of Brescia, Brescia, Italy) for her valuable contribution on the revision of the paper concepts and structure. We thank Dr. Irene Tosato (University of Brescia, Brescia, Italy) for her work of proofreading and English revision. We express our gratitude to Mr. Gabriele Merli (University of Pavia, Pavia, Italy) for his advice and kind helpfulness and to Professor Luca Pavarino (University of Pavia, Pavia, Italy), in charge of the new high-performance computing cluster (HPC) of the University of Pavia. We also acknowledge the World Climate Research Programme's Working Group on Regional Climate, and the Working Group on Coupled Modelling, former coordinating body of CORDEX and responsible panel for CMIP5. We also acknowledge the Earth System Grid Federation infrastructure an international effort led by the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison, the European Network for Earth System Modelling and other partners in the Global Organisation for Earth System Science Portals (GO-ESSP). We acknowledge the E-OBS dataset from the EU-FP6 project UERRA (https://www.uerra.eu) and the Copernicus Climate Change Service, and the data providers in the ECA&D project (https://www.ecad.eu).
Authors' contribution GG, GSp, SP and GSc contributed to the conceptualisation and the design of the work. MC and PG analysed the model and performed the numerical simulations. AW and ARD performed distributional and climatic data acquisition and interpretation. GG and GSp drafted the manuscript.
Funding Open access funding provided by Università degli Studi di Brescia within the CRUI-CARE Agreement. This research has been funded by the Julius Kühn-Institut (JKI), Institute for National and International Plant Health (Braunschweig, Germany) under the cooperation agreement ''Assessment of the impacts of plant pests under climate change''. This research has been supported by Fondazione Cariplo (Italy) and Regione Lombardia (Italy) under the project: ''La salute della persona: lo sviluppo e la valorizzazione della conoscenza per la prevenzione, la diagnosi precoce e le terapie personalizzate''. Grant Emblematici Maggiori 2015-1080.

Declarations
Conflict of interests The authors declare no conflict of interest.
Data availability Raw data used in the present study are available upon request.
Code availability Custom software simulating C. capitata population dynamics is available upon request.
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://creativecommons.org/licenses/by/4.0/.