Biological invasion risk assessment of Tuta absoluta: mechanistic versus correlative methods

The capacity to assess invasion risk from potential crop pests before invasion of new regions globally would be invaluable, but this requires the ability to predict accurately their potential geographic range and relative abundance in novel areas. This may be unachievable using de facto standard correlative methods as shown for the South American tomato pinworm Tuta absoluta, a serious insect pest of tomato native to South America. Its global invasive potential was not identified until after rapid invasion of Europe, followed by Africa and parts of Asia where it has become a major food security problem on solanaceous crops. Early prospective assessment of its potential range is possible using physiologically based demographic modeling that would have identified knowledge gaps in T. absoluta biology at low temperatures. Physiologically based demographic models (PBDMs) realistically capture the weather-driven biology in a mechanistic way allowing evaluation of invasive risk in novel areas and climes including climate change. PBDMs explain the biological bases for the geographic distribution, are generally applicable to species of any taxa, are not limited to terrestrial ecosystems, and hence can be extended to support ecological risk modeling in aquatic ecosystems. PBDMs address a lack of unified general methods for assessing and managing invasive species that has limited invasion biology from becoming a more predictive science.

followed by Africa and parts of Asia where it has become a major food security problem on solanaceous crops. Early prospective assessment of its potential range is possible using physiologically based demographic modeling that would have identified knowledge gaps in T. absoluta biology at low temperatures. Physiologically based demographic models (PBDMs) realistically capture the weather-driven biology in a mechanistic way allowing evaluation of invasive risk in novel areas and climes including climate change. PBDMs explain the biological bases for the geographic distribution, are generally applicable to species of any taxa, are not limited to terrestrial

Introduction
Predicting the potential geographic distribution and relative abundance of invasive species prior to invasion is key to the assessment of risk and management (Bradshaw et al. 2016;Lenzner et al. 2019). A bellwether failure of pre-invasion assessment is the South American tomato pinworm Tuta absoluta (Meyrick) (Lepidoptera: Gelechiidae), a serious native insect pest of tomato in South America (Desneux et al. 2010). It was not identified as a serious threat by the European Union (EU), the USA, or other tomato-growing areas, until it invaded Spain in 2006, from where it spread rapidly across Europe, Africa, and Asia (Pratt et al. 2017;Biondi et al. 2018;Han et al. 2019). This oligophagous leaf miner and its primary tomato host are native to South America, specifically Peru where both originated (Blanca et al. 2015;Biondi et al. 2018). The species was described from individuals collected in Huancayo in the central Andean highlands of Peru (Biondi et al. 2018), a location at -12.065°latitude, -75.205°longitude, and 3245 m altitude with a cold semi-arid steppe climate. The pest develops on solanaceous species, the most suitable ones being tomato, potato, and European black nightshade (Solanum nigrum) (Biondi et al. 2018), but it is also recorded from other cultivated Solanaceae (e.g., eggplant) (Cherif et al. 2019;Sylla et al. 2019). Tomato is widely grown in the western Palearctic (Europe and Mediterranean Basin) and North America (Supplementary Information, SI, Figs. S1 and S2).
Given its host range of agriculturally important crops and its status as a key pest of tomato in South America, it is surprising that large tomato producing entities with well-organized plant quarantine infrastructure such as the EU and the USA failed to consider the invasion risk the moth posed (Biondi et al. 2018;Desneux et al. 2010). In 2006, moths from a single Chilean population invaded Spain (Guillemaud et al. 2015), and in ten years expanded its distribution to 60% of the global tomato growing area from northern Europe to South Africa, and from West Africa to South Asia. This area includes six of ten top tomato-growing countries in the Palearctic (i.e., India, Turkey, Egypt, Iran, Italy, Spain). The moth has recently invaded China (the top producer) (Zhang et al. 2019;Li et al. 2020) but not the USA (fourth) and Mexico (ninth) (Biondi and Desneux 2019) although it is present in Haiti (Verheggen and Fontus 2019). To understand this invasion, we need to know why this species was able switch from being a local pest of tomato in South America to become a global food security threat in climes beyond its recorded distribution. Accurate risk analysis of T. absoluta before 2006 would have triggered global quarantine measures, preventing its unprecedented rapid invasion.
Three post hoc attempts were made between 2010 and 2019 to estimate its potential global distribution as affected by extant weather (Desneux et al. 2010;Han et al. 2019) and climate change (Santana et al. 2019) using the ecological niche modeling software CLI-MEX (Sutherst and Maywald 1985). CLIMEX and other correlative methods are de facto standards that use presence-absence data to find climatic correlates as indices of favorability to predict the prospective distribution of ectotherm species (Elith 2017). Such studies may not contain records from the full climatic range of the species, and further lack mechanistic biological underpinnings, hence as occurred for T. absoluta may fail to predict and explain the full invasive range in novel areas such as the Palearctic.
Here, a mechanistic physiologically-based demographic model (PBDM) (Gutierrez 1996;Ponti et al. 2015a, b) that captures the weather-driven biology of the moth is developed to explain and predict prospectively its invasiveness in the western Palearctic and in uninvaded areas such as the USA and Mexico under extant and climate change scenarios.

Materials and methods
The PBDM for T. absoluta is mathematically the same as the model of the daily dynamics developed for the light brown apple moth Epiphyas postvittana (Walk.) (Lepidoptera: Tortricidae) (Gutierrez et al. 2010), and consists of two components: • A general age-structured distributed maturation time model of the daily dynamics of the moth (Eq. 1); and • Weather-driven biodemographic functions (BDFs) that parameterize the age-structured biology T. absoluta, and capture the time varying effects of temperature and relative humidity on vital rates (sensu Hughes and Gilbert 1968;Gutierrez 1992) (see Fig. 1).
Biological models will always be incomplete, and the field and laboratory data used to parameterize them may have differences due to experimental methods and other factors. Hence the use of data from the literature to parameterize the BDFs must be assessed for internal consistency (data outliers) and for consistency with other data sets for the same factor (see below). Well parameterized PBDMs have given robust predictions about the time-place varying densities of species, and in the aggregate across a lattice of locations with different weather; results that are mapped to predict their geographic distribution and relative abundance (e.g., Gutierrez and Ponti 2013). In sequence we review the dynamics model and the imbedded biodemographic functions that capture the weather-driven biology of T. absoluta.    1 Thermal biology of Tuta absoluta: development, reproduction, and mortality as a function of temperature. a Developmental rate of the egg-pupal stages as a function of temperature (data from Barrientos et al. 1998;Krechemer and Foerster 2015;Martins et al. 2016;Campos et al. 2020). b Agespecific fecundity (data from Marcano 1995). c Normalized gross fecundity as a function of temperature (data from Marcano 1995;Krechemer and Foerster 2015). d Oviposition per female as a function of age and temperature (Eq. 8). e Temperaturedependent mortality (data from Van Damme et al. 2015;Krechemer and Foerster 2015;Martins et al. 2016;Kahrer et al. 2019) with the line indicating the fitted polynomial function used in the model The dynamics model Cohorts of individuals in the population experience different weather and hence have different patterns of time-varying means and variance for developmental times, and this is easily modeled using a weatherdriven distributed maturation time model (Eq. 1) based on Manetsch (1976) and Vansickle (1977). We note that other dynamics models could also be used (see Gutierrez 1996;Di Cola et al. 1999;Buffoni and Pasquali 2007). The four life stages (left superscript s = E, L, P, A for the egg, larval, pupal, and adult stages respectively) are modeled by a system of s k equations (Eq. 1) that describes the dynamics of the ith of s k age class. The distribution of developmental times in the absence of mortality is described by an Erlang distribution with parameter k equal to the number of within stage age classes required to approximate the observed mean ( s D) and variance ( s var) of developmental time of cohort members. Following the notation of Di Cola et al. (1999, p. 523), the ith age class of stage s is modeled as follows: In terms of flux, s n i ðtÞ ¼ s N i ðtÞ s m i ðtÞ where s m i ðtÞ ¼ Note that new eggs flow into the first age class of the egg stage with the outflow from the kth age class of the previous stage flowing into the first age class of the next stage (e.g., egg to larval, etc.). Ignoring the stage superscript, N i is the number density of the ith age class, dt is the change in time (e.g., a day), D is the expected mean stage developmental time, D x is the daily increment of physiological age in degree days for all life stages that varies with temperature, and l i (t) is the proportional age-specific net loss rate as modified by temperature, age, net migration, and mortality due to natural enemies. The Euler integration scheme implemented for the discretized model is found in Gutierrez (1996, page 157).

Developmental rate
The time step in our PBDM is a day of variable length in physiological time units. Shorter time steps could be used, but the adage holds that ''one needn't measure something with a micrometer if one must cut it with an axe''.
Development may also be modified by nutrition, dormancy, and other factors, but are not included in the model (Gutierrez 1992(Gutierrez , 1996. The developmental rate per day at temperature T is the reciprocal of developmental time in days (d): r(T) = 1/d(T). Developmental times and rates at different temperatures ( Fig. 1a) are for the total egg, larval, and pupal period (D E-P ). We could have modeled the developmental rate of each life stage separately (Marchioro et al. 2017), but accurate data for the leaf miner is difficult to obtain, and daily observations for life stages with shorter developmental times tend to have larger error rates. The data for the E-P period from different studies in the literature are largely consistent (Fig. 1a), and a non-linear model (Eq. 3) was iteratively fit to the data (Gutierrez and Ponti 2013) (v 2 = 0.000312, reduced v 2 = 0.000015; computed using the lmfit Python package version 1.0.2, Newville et al. 2014).
The parameters of Eq. 3 are defined as follows: T(t) is the mean daily temperature on day t, the coefficient a = 0.0024 is the initial slope of the function (Bieri et al. 1983), h L = 7.9°C is the lower thermal threshold for development, h U = 34.95°C is approximately the highest temperature before the function begins to decline to zero, and b = 3.95 is a constant. A similar function was proposed by Briére et al. (1999). We assume that h L is the same for all life stages.
Temperatures experienced by T. absoluta by the egg and adult stages are ambient, while larvae and pupae in leaf mines experience temperatures that differ from ambient, and we must account for the difference. For example, Baumgärtner and Severini (1987) developed a microenvironment simulator based on a daily balance equation of a wet body to estimate temperatures experienced by apple leaf miner moth (Phyllonorycter blancardella) larvae and pupae. Pincebourde and Casas (2006) developed a heat budget model for P. blancardella, and found that mine temperatures were higher than leaf surface temperatures and varied with irradiance. Their heat budget model with irradiance levels ranging from 0 to 1300 lmol m -2 s -1 , temperatures ranging from 11 to 36°C, and relative humidity ranging from 20 to 90%, enabled them to map changes in leaf mine temperatures. Unfortunately, due to lack of appropriate data, we were unable to implement their heat budget model explicitly in our PBDM, but we could estimate the effects of temperature and irradiance from their simulation data (i.e., Fig. 8 in Pincebourde and Casas 2006). A conversion of 1 W m -2 & 4.6 lmol m -2s -1 was made to accommodate the units in our weather data (Sager and McFarlane 1997), and a linear multiple regression model was fit to the data to estimate the change in temperature in mines relative to ambient (DT m ) at time t, temperature T (i.e.,°C), relative humidity (0 \ RH \ 1), and irradiance (Wm -2 ): The daily mean temperature (T _ ) experienced by larvae and pupae in leaf mines is and substituted for T in Eq. 3. Increases in leaf surface temperatures relative to mine temperatures across all irradiation values are * 0:65DT m . Our PBDM is a discrete time model and daily increments of degree days (D x ) are computed as the proportions of the egg to new adult developmental time completed during day t at temperature T. Using the stage notation in Eq. 1, the thermal constant for substage s D in degree days (dd above h L ) is computed in the linear range of favorable temperatures as s D = s d(T)(T-h L ). Specifically, average E D for the egg stage is 82.3dd, L D = 218.8dd, P D = 122.9dd, and A D = 189.2dd is adult female longevity. The preoviposition period under nonlimiting host conditions is approximately 30.4dd (Bentancourt et al. 1996;Foerster 2015, 2017). A life stages completes development when P D x (t) = s D.

Oviposition
Oviposition occurs on the leaf surface and is assumed to be at ambient temperatures. Per capita oviposition was estimated by several authors, but we used the data from Marcano (1995) (Fig. 1b) to estimate the maximum age-dependent oviposition profile at the optimum temperature of 25°C that was fit with the leftbiased Bieri et al. (1983) function, where adult age is in days (d).
At T opt , D x is 17.1dd per day, and hence d in Eq. 6 can be converted to adult physiological age (i.e., days 9 17.1dd). Reversing the process, we can estimate the number of adults ( A N(t,d)) of age d at time t. Fecundity also varies with temperatures ( Fig. 1b), and while studies report different values, when normalized within each data set, a consistent relationship with temperatures arises. Hence, to correct Eq. 6 for the effects of temperature, we model the normalized data from Marcano (1995) and Krechemer and Foerster (2015) with a concave right-skewed function (Eq. 7) ( Fig. 1c) 0 / fec ðTðtÞÞ ¼ 0:0665ðTðtÞ À 7:9Þ 1 þ 2:2 ðTðtÞÀ27:5Þ 1 ð7Þ where 7.9°C is the lower thermal threshold for oviposition and * 27.5°C is a fitted constant where the oviposition rate begins to decline to zero at * 30°C. From this, we generate a three-dimensional fecundity function on age and temperature ( Fig. 1d). Note the lower threshold for oviposition is the same as for development, but the upper threshold is * 3°C lower. Thus, the number of eggs (E) deposited by the population at time t at temperature T with assumed sex ratio sr = 0.5 is As indicated, new eggs enter the first age class of the egg population dynamics model (Eq. 1).

Mortality
There are many sources of mortality, but in the model, we include only the major effects of temperature and a composite mortality function for predation, host finding, and other biotic factors.
Temperature-dependent mortality Studies under simulated greenhouse conditions (Cuthbertson et al. 2013) and in the laboratory (Bentancourt et al. 1996;Barrientos et al. 1998;Van Damme et al. 2015;Krechemer and Foerster 2015;Martins et al. 2016;Ö zgökçe et al. 2016) estimated temperature-dependent mortality rates. However, only Kahrer (2015) and Kahrer et al. (2019) estimated the degree of cold hardening that can occur in the moth allowing it to survive at temperatures well below 7.9°C; conditions that normally occur to varying degrees during the fallwinter-spring period in different parts of T. absoluta's expanding range. Hence, the effects of daily temperature on daily mortality rates across all temperatures (l T (T)) were estimated from literature data (Van Damme et al. 2015;Krechemer and Foerster 2015;Martins et al. 2016;Kahrer et al. 2019), and were fit with a polynomial model with biological limits (Eq. 9, R 2 [ 0.96) (see Fig. 1e). A simpler function could have been used, but we chose accuracy of fit rather than elegance with the extra decimal places given for reproducibility by others. Note that ambient temperatures are used in Eq. 9 for the egg and adult stages, while estimated mine temperatures (i.e., T ¼ T _ ) are used for the larval and pupal stages.
Other sources of mortality Miranda et al. (1998) conducted field life table analyses in Brazil to estimate survivorship of T. absoluta immature stages, but the results are time and place specific. A type II model (Eq. 10, Gutierrez 1996) is used to capture the effects of composite sources of predation/parasitism mortality per dd at temperature T (i.e., D x (t)), and density of each stage ( s N(t)), and to keep the dynamics within reasonable bounds (Eq. 10).
where c = 0.00025 is the assumed apparency rate of moth stages per dd.
To estimate the joint mortality, we write the interaction as survivorship terms for age i as 1 À s l i ðtÞ ¼ ð1 À s l T ðTðtÞÞ Á ð1 À s l _ ðtÞÞ and then solve for s l i ðtÞ ¼ s l T ðTðtÞÞ þ s l _ ðtÞ À s l T ðTðtÞÞ s l _ ðtÞ in Eq. 1.

Weather data
Daily weather data for maximum and minimum temperature, rainfall and solar radiation (MJ/ m 2 /day 2 ) were used to drive PBDM simulations. Weather for present climate simulations was from AgMERRA, the global baseline forcing dataset of the Agricultural Model Inter-comparison and Improvement Project (AgMIP, http://www.agmip.org/) (Ruane et al. 2015) (https://data.giss.nasa.gov/impacts/ agmipcf/). AgMERRA is a reanalysis of weather observations (Rienecker et al. 2011) combined with observational datasets from in situ observational networks and satellites (Ruane et al. 2015). Daily AgMERRA weather at * 25 km spatial resolution was used for the period 1 January 1990 to 31 December 2010 for 17,791 lattice cells across the Euro-Mediterranean region and 15,843 for the USA and Mexico. The climate data summary for Huancayo, Peru is from CLIMATE-DATA.ORG (Climate-Data.org 2019).
Daily weather for climate change simulations in the Euro-Mediterranean region was from the A1B regional climate change scenario with * 30 km resolution for the Euro-Mediterranean region developed by Dell'Aquila et al. (2012) using the regional climate model PROTHEUS (Artale et al. 2010) to refine coarser (* 200 km resolution) global climate simulations (http://www.nextdataproject.it/?q=en/content/ regional-climate-simulations-enea). PROTHEUS is a coupled atmosphere-ocean regional model that allows simulation of local extremes of weather via the inclusion of a fine-scale representation of topography and the influence of the Mediterranean Sea (Artale et al. 2010). The A1B scenario is towards the middle of the IPCC (IPCC 2014) range of greenhouse gas (GHG) forcing scenarios (Giorgi and Bi 2005), and the uncertainty associated with climate model predictions forced using the A1B scenario is low for the Mediterranean region relative to the rest of the globe (Gualdi et al. 2013). Based on this downscaled climate change scenario, a daily weather dataset was developed for the period 1 January 2040 to 31 December 2050 for 10,598 lattice cells in the Euro-Mediterranean region.
Daily climate change simulations for the USA and Mexico were driven by the NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) dataset (Thrasher et al. 2012) at * 25 km resolution, derived from global climate simulations (Taylor et al. 2012) of the Max Planck Institute Earth System Model low resolution (MPI-ESM-LR) model (https://www. nccs.nasa.gov/services/data-collections/land-basedproducts/nex-gddp). Global MPI-ESM-LR climate simulations were forced by the Representative Concentration Pathway 8.5 (RCP 8.5) scenario (Riahi et al. 2011) that corresponds to a range of warming similar to A1B (Rogelj et al. 2012). Sheffield et al. (2013) evaluated historical simulations of North American climate using continental metrics of bias relative to weather observations and showed that MPI-ESM-LR is the top ranked among the core set of 17 global climate models considered. Based on this downscaled climate change scenario, a daily weather dataset was developed for the period 1 January 2045-31 December 2075 for 20,355 lattice cells in the USA and Mexico.
The downscaled climate change scenarios used in this study include high-resolution daily weather data designed for assessing climate change impacts on processes that are sensitive to fine-scale climate and local topography, such as biological processes of poikilotherm organisms like T. absoluta and other invasive pest insects. Downscaling (i.e., increasing the spatial resolution) of global climate model output addresses two primary limitations: the relatively coarse spatial resolution of global climate models (e.g., hundreds of km) and their statistical bias when compared with observations (Thrasher et al. 2012). We used climate scenarios that were tested in previous PBDM analyses Gutierrez et al. 2019) in an effort to contribute towards providing credible information on regional climate changes, impacts, and risks for stakeholders and policy-makers (Eyring et al. 2019).

PBDM simulation and GIS analysis
The PBDM for T. absoluta was run continuously across all lattice cells for the specified periods, with no migration occurring between lattice cells. The starting day for all runs was 1 January of the first year with an initial nominal density of 0.1 m -2 for each life stage at all locations. While the model predicts the daily density of all stages of T. absoluta, only cumulative pupal density per square meter per year was used as the main metric (indices) of favorability. The simulation data were geo-referenced and written to batch files by year for mapping and analysis. Because the model must equilibrate, the first year's results for all simulation runs were not used to compute lattice cell averages, standard deviations, and coefficient of variations (%).
The open source geographic information system (GIS) GRASS (Neteler et al. 2012) originally developed by the United States Army Corps of Engineers was used for geospatial data management, analysis, and mapping using bicubic spline interpolation on a 3-km raster grid. The version of GRASS used is maintained and further developed by the GRASS Development Team (2020) (see http://grass.osgeo. org). Further processing and plotting of raster layers generated by GRASS were done using R (R Core Team 2020) with packages raster (Hijmans 2020) and ggplot2 (Wickham 2016).
Geographic information layers for climate-limited tomato distribution are from the Global Agro-ecological Zones (GAEZ v3.0) dataset (IIASA and FAO 2012) (http://www.fao.org/nr/gaez/en/). We used two GAEZ temperature-constrained potential distribution maps for tomato (irrigated, with intermediate input level, and with a 105 days growth cycle, for average climate : the distribution of temperatesubtropical tomato cultivars and the distribution of tropical-subtropical tomato cultivars in open field conditions. We combined the two distributions, including only the highest favorability class in the maps. The combined tomato distribution (SI, Figs. S1 and S2) was used in this study when mapping the geographic distribution of T. absoluta to correct for the distribution of tomato. Data on global tomato production are from FAOSTAT (see http://faostat.fao. org/).
Using GRASS GIS, CLIMEX maps and data from other studies were digitized from PDFs of published papers, orthorectified, and converted to the Albers Equal Area projection used in this study.

Results
Modeling T. absoluta dynamics inside leaf mines During the year, weather in temperate areas may vary from hot summers to cold winters. Under either extreme, conditions may prove unfavorable for the development, reproduction, and survival of moth life stages. Simulated pupal dynamics as affected by weather are used as a metric of T. absoluta abundance and as a measure of location favorability.
We use weather data for Sacramento, California (USA), an area where field tomato is widely grown, to illustrate the interacting effects of temperature and irradiance on Tuta dynamics. Summers in Sacramento are hot and winters may be cool to cold with occasional frosts. The 2000-2010 simulation for Sacramento in Fig. 2 is shown as a detailed subset for 2009 and 2010 in Fig. 3. Using ambient temperatures, the seasonal dynamics of the moth (pupae m -2 ) in Sacramento are affected by high temperatures, while during the winter period, oviposition stops but low numbers of pupae survive (\ one pupa m -2 ; Fig. 2a, solid line). However, if we include the effects of full irradiance on the outer margins of the canopy, the higher mine temperatures increase larval and pupal mortality during summer (not shown), driving pupal populations to near zero ( Fig. 2a; shaded area). Tuta populations rebound only slightly with cooler fall weather, but pupal survival during the winter period is low.
However, in the field, adult moths select more favorable conditions within the canopy for oviposition and for the development of their progeny (Gomide et al. 2001;Torres et al. 2001;Cocco et al. 2015). Using Beer's Law with a light extinction coefficient of 0.8 and assuming a full canopy (LAI = 4.5), 0.5 light penetration would occur at * 19.15% of canopy depth on fully expanded leaves. Hence, comparing the effects of ambient temperatures without the effect of irradiance and with 0.5 irradiance in the canopy (Fig. 2b), pupal dynamics of the two scenarios are nearly identical. Further, earlier appearance of pupae is predicted in spring at 0.5 irradiance compared to ambient. This is similar to that found by Baumgärtner and Severini (1987) for the apple leaf miner at Chur, Switzerland, where changes in the microenvironment of mines influenced the emergence times, rates, and the duration of the flights with the magnitude dependent on local environmental conditions. This suggests tradeoffs occur between increases in developmental rates and temperature-dependent mortality. Variations of temperature and irradiance-driven dynamics occur throughout the moth's range, and hence for all georeferenced lattice cells in all runs of the model, we use 0.5 irradiance as a standard including the additional effects of relative humidity (see Eq. 3).
The daily weather and the simulated population dynamics for Sacramento in Fig. 3 illustrate some of the complexity of factors affecting the dynamics of T.  Figure 3a and b illustrate the maximum and minimum temperatures and solar radiation used to drive the model, but relative humidity is not illustrated. Figure 3c shows the dynamics of moth life stages as affected by weather (see Fig. 3d; Eq. 9). Mortality occurs due to cold temperatures below the lower thermal threshold (black line) and high temperatures above the optimum temperature for oviposition (gray line) ( Fig. 3d; see Fig. 1e), but the mortality is not severe enough to limit invasion by the moth (see below Sect. 3.3). Across the geographic landscape, each location (lattice cell) and year, weather would produce distinctive dynamical population dynamics, with the yearly total of pupae used as a metric of favorability for the moth.

Invasiveness in the Euro-Mediterranean region
Highest predicted populations of T. absoluta are in coastal areas of southern Spain and North Africa, with populations decreasing into colder areas of Europe ( Fig. 4a, b, c; mean pupae, CV(%), std respectively). The moth's distribution coincides with the known distribution of temperate and tropical field tomato (see SI, Fig. S1) grown under open field conditions in all areas except in hotter regions of North Africa and the Levant (Fig. 4d, e, f). The limits of distribution can be approximated by the annual cumulative daily mortality rates at low (l Ã \7:9 C ¼ P l \7:9 C ) and high (l Ã [ 25 C ¼ P l [ 25 C ) temperatures (Fig. 4g, h respectively). Cold temperatures are limiting in temperate areas, with the geographic range of the moth limited above l Ã \7:9 C * 3.5 (see Fig. 4g). In contrast, high-temperature mortality is limiting in hot-dry  S1). The bottom row shows the cumulative yearly temperature mortality rates at low (l Ã \7:9 C \ 3.5) (g) and at high temperatures (l Ã [ 25 C ) (h), and the total mm rainfall per year (i). The full data interval in (g) is [0.0, 58.0] with values [ 3.5 mapped using the same darkest shade of red pupae m À2 ¼ À 42:558 þ 0:029dda þ 0:023mm rain À 0:399 l Ã \7:9 C À 6:616l Ã The computed t-values are 102.36 for dda [ 7.9°C , 25.24 for mm rain , -16.39 for l Ã \7:9 C , and -89.06 for l Ã [ 25 C . The positive signs of the regression coefficients for dda and mm rain infer the positive effects of a longer favorable season, and the negative signs for cumulative low-and high-temperature mortality infer the limiting effects of cold and hot weather. The effect of high temperature is 16.5-fold greater per unit mortality than that of low temperature. This is due primarily to cold hardiness, and a low upper threshold for reproduction (* 30°C) (see Fig. 1). The summed low and high temperature mortality rates are qualitative metrics of favorability.

Prospective invasiveness in USA and Mexico
Tuta absoluta has not invaded North America, and hence our analysis is prospective of its potential geographic distribution and abundance. Highest populations are predicted in more temperate and tropical areas of Mexico, with moderate populations predicted in the southeastern USA, and in the Central Valley and southern near-coastal areas of California (Fig. 5a, b,  c). The hotter dryer areas of North America are unfavorable, as are higher colder elevations and more northern reaches of the continent (Fig. 5a, b, c). The area of favorability for the moth is again quite similar to the distribution of temperate-tropical tomato cultivation (Fig. 5d vs. a). The limits to its geographic distribution are illustrated using the same constraints of average cumulative daily low (l Ã \7:9 C \3:5; Fig. 5e) and high (l Ã [ 25 C ; Fig. 5f) temperature mortality rates found for the Euro-Mediterranean region. Absent irrigation, rainfall is also limiting in hot-dry areas of southwestern USA and Mexico (Fig. 5g). Note that in southern Texas and adjacent areas of Mexico, high temperature effects on mortality and reproduction are limiting despite abundant rainfall (see Fig. 5d and compare Fig. 5f vs. Fig. 5g).

Geographic distribution with climate warming
In Europe, under climate warming as estimated by climate models, the moth's distribution expands considerably northward and eastward to Eurasia (Fig. 6a) as seen by the recession of the l Ã \7:9 C * 3.5 boundary (Fig. 6b vs. Fig. 4g), while in southern reaches of North Africa and the Levant, the area of favorability shrinks due to high temperatures (l Ã [ 25 C ; Fig. 6c vs. Fig. 4h).
In the USA and Mexico, the area of favorability shrinks profoundly in the southeastern USA and northern Mexico (Fig. 6d vs. Fig. 5a). The l Ã \7:9 C * 3.5 boundary recedes northward (Fig. 6e  vs. Fig. 5e), while the cumulative high-temperature mortality rates increase in the desert areas of the southwestern USA and northwestern Mexico, and in the southeastern USA including Texas and eastern areas of Mexico including the Yucatan Peninsula (Fig. 6f vs. Fig. 5f). Areas of favorability increase in higher elevations of the western US states.

Discussion
Before invading the Euro-Mediterranean region, T. absoluta was not a regulated quarantine pest, and tomato imports were not inspected for this pest at international borders in the EU and USA (Biondi et al. 2018). This occurred despite its known range expansion in South America and severe pest status on tomato in areas of Brazil having warm climates similar to Mediterranean coastal areas that were quickly invaded after it was first recorded in Spain (Desneux et al. 2010). Most Brazilian tomato is grown in three areas (Melo 1992): the Northeast with semi-arid climate (da Silva et al. 2005); the Southeast with tropical climate (Melo 1992); and the Center-West with Cerrado humid subtropical climate (Alvares et al. 2013). Temperatures in these areas of Brazil are on average 8-10°C warmer than the cold semi-arid climate of the Andean highlands (Biondi et al. 2018).
But how best to examine T. absoluta biology in determining its geographic range and relative abundance was an open critical question. Based on records from these warm areas, its establishment at higher latitudes in Europe was considered unlikely (Desneux et al. 2010), but ultimately proved wrong. Among the approaches used to predict the moth's invasiveness were cellular automata that posited the importance of relative humidity and host plant distribution (Guimapi et al. 2016), and the widely used correlative environmental niche model (ENM) software CLIMEX. The CLIMEX software predicts the potential range of the moth using monthly averaged climatic correlates of occurrence in its known range (Desneux et al. 2010;Santana et al. 2019), and hence predictions are influenced by the distribution of these records (Elith 2017). The roots of CLIMEX are the early studies of Fitzpatrick and Nix (1970) who developed growth indices to estimate the climatic limits of Australian grasslands, and Gutierrez et al. (1974) and Gutierrez and Yaninek (1983) who used growth index models to capture the climatic limits of aphids in southeastern Australia. Elements of these initial studies (i.e., growth indices) are found in PBDMs (Gutierrez 1996).
Three CLIMEX-based assessments of T. absoluta potential range were made based on available thermal development data (Desneux et al. 2010;Santana et al. 2019;Han et al. 2019). At the onset of global invasion in 2006, data on the thermal biology of T. absoluta were available only in the range of favorable temperatures above 12°C, with data below 12°C reported starting in 2015, after invasion of central Europe occurred (Desneux et al. 2010) when its overwintering potential became evident (Van Damme et al. 2015). In 2010, the first CLIMEX analysis (Desneux et al. 2010) projected that only coastal southern Europe would be distribution; e average yearly cumulative mortality rates for low temperatures clipped to l Ã \7:9 C B 3.5 (values [ 3.5 in dark red); f average yearly cumulative daily mortality rates for high temperatures (l Ã [ 25 C ); and g average annual rainfall (mm) favorable ( Fig. 7a), this despite the pest having been recorded from central Europe (Desneux et al. 2010) ( Fig. 7b and SI, Fig. S3). Similarly, the second CLIMEX model analysis in 2019 (Han et al. 2019) was unable to predict most areas of known occurrence of T. absoluta in Europe (SI, Fig. S4). The same year, a third CLIMEX analysis with insights of its invaded range (Santana et al. 2019) used a lower thermal threshold of 7°C and a wider optimal range for development (14-25°C) based on data from Martins et al. (2016) not included in the previous CLIMEX studies ( Fig. 8a and b, SI, Fig. S4). This analysis predicted a Euro-Mediterranean and North American range for T. absoluta that compares favorably with the predictions of our PBDM ( Fig. 8c and d, SI Figs. S5, S6; see below). Note that the PBDM has a finer scale for the index of favorability, did not require occurrence data, and the full population dynamics for any location are available (e.g., Fig. 2).
Additional major limitations of the CLIMEX approach (Elith 2017) and by extension other correlative approaches used for assessing invasive species risk (Yates et al. 2018), include the lack of information about the phenology and vital rates that determine the population dynamics of species. These limitations preclude a detailed understanding of the observed distribution and relative abundance, and hence limit their transferability in geographic space and time to new areas and/or under climate change (Yates et al. 2018). Because correlative methods rely on occurrence records to assess invasive potential, suitable records capturing the full range of plasticity to environmental variables may not be available until after the potential invasive range of an invasive species has been occupied. Rossini et al. (2019) modeled the temperaturedependent age-structured dynamics of T. absoluta, but did not include leaf-mine microclimate or temperature-dependent mortality, and did not attempt assessment of its geographic distribution. Our mechanistic PBDM representation of the weather-driven biology and age-structured population dynamics (Gutierrez 1996) of T. absoluta captures the full range of its thermal biology, including the effects of relative humidity and solar radiation on microclimate temperatures experienced by the larvae and pupae inside leaf mines. Unlike correlative methods, PBDMs are not biased by the species' range or by spatial autocorrelation in climatic and distribution data (Journé et al. 2019;Briscoe et al. 2019). These attributes make our PBDM transferable to novel locations and conditions. Further, the PBDM for Tuta required no calibration, because the weather-driven biodemographic functions captured the moth's biology. When used in a GIS context as done here, PBDMs bridge the gap between bottom-up (e.g., field experiments) and top-down approaches (e.g., correlative ecological niche methods) to invasive species assessment (see SI, Discussion and Fig. S7).
Temperature is a dominant driver of T. absoluta dynamics, and recently Kahrer et al. (2019) and Campos et al. (2020) explored the biology in a range of low temperatures critical for assessing its overwintering survival in temperate regions, and hence its prospective geographic range and relative abundance. The moth's modest developmental threshold (7.9°C), and facultative diapause (Campos et al. 2020), combined with its high degree of cold hardiness (Kahrer et al. 2019) enable its northward range expansion that is ultimately limited by cold in northern areas roughly delimited at l Ã \7:9 C * 3.5. In southern areas, high temperatures affect reproduction and survival rates, particularly in hot-dry desert areas of North Africa, USA, and Mexico.
Had adequate data been available, and a PBDMbased pest risk analysis (Gutierrez 1996;Ponti et al. 2015a, b) been performed before T. absoluta invaded    (Riudavets et al. 2016) could have prevented its crossing of the ocean barrier to Spain. However, once established in Spain, quarantine measures proved ineffective, as range expansion in Europe, Africa, and Asia (Desneux et al. 2010) was aided by human spread (McNitt et al. 2019) and likely by its ability for long wind-assisted flight dispersal (Desneux et al. 2010;Biondi et al. 2018) as documented for another moth in the same Family Gelechiidae, the pink bollworm Pectinophora gossypiella (Stern and Sevacherian 1978). The invasive range of T. absoluta now includes many developing countries that are under-resourced to implement policies for preventing and managing biological invasions (Measey et al. 2019). In eight years, the species invaded the whole African continent (Biondi et al. 2018;Biondi and Desneux 2019) where it threatens the food security and livelihoods of subsistence farmers of solanaceous crops (Pratt et al. 2017;Aigbedion-Atalor et al. 2019). In East Africa alone, estimated annual losses to T. absoluta in mixed small farming systems are estimated to be 69.6-79.4 million USD (Pratt et al. 2017). Similar threats to food security and estimates of high economic losses occur in developing Asian countries (Biondi et al. 2018;Bahadur and Bikram 2019;Biondi and Desneux 2019;Han et al. 2019). The majority of farmers in Africa and Asia are smallholders (Hruska 2019), most of whom before recent major biological invasions such as T. absoluta and Spodoptera frugiperda, used no pesticides (Rwomushana et al. 2019;Hruska 2019).

Conclusions
Global society must do a better job of preventing and mitigating biological invasions that may only worsen under global change. The first tasks are to identify and characterize the risk posed; secondly, to institute measures for preventing their introduction; and thirdly, to have institutions and mechanisms to control them if introductions occur. Some species, however, may be particularly difficult to identify as risks, as they may expand their host range and increase their invasive potential only after invasion when released from native biotic regulation, as occurred with klamath weed (Hypericum perforatum), cottony cushion scale (Icerya purchasi), and cassava mealybug (Phenacoccus manihoti)-species that were quickly brought under biological control (Huffaker and Messenger 1976; Herren and Neuenschwander 1991). However, known pests of crops that are grown globally should be characterized and analyzed as potential invasive pests: T. absoluta is a bellwether case.
Given identification of risk, the PBDM approach offers a solid mechanistic pest risk assessment (Ponti et al. 2015a) of invasive species (Ponti et al. 2015b) lacking in most correlative approaches used to predict biological invasions (Novoa et al. 2020). Specifically, PBDMs have similar subcomponents (Gutierrez and Ponti 2013) that speed up identification of gaps in biological data and streamline data collection (Ponti et al. 2015a), and require relatively modest data and funding to parameterize and develop sound mechanistic descriptions of the biology of an invasive species (Gutierrez and Ponti 2013). In contrast to correlative ENMs (Yates et al. 2018), well parameterized PBDMs are transferable and can be used to gauge the invasiveness of species in new areas and/or under climate change using management-relevant metrics such as their phenology, relative abundance, and geographic distribution-before invasion occurs (e.g., North and Central America). PBDMs inform the biological mechanisms for the geographic distribution of species, including their prospective invasive range.
Compared to phenological models such as that developed by Suppo et al. (2020) for the invasive box tree moth Cydalima perspectalis, PBDMs enable assessing age-structured population dynamics that have distributed maturation times and are driven by temperature (including extremes), humidity, and other factors (e.g., nutrition) with time-varying effects on development, reproduction, and mortality-in addition to phenology. Because PBDMs generalize the biology of species across trophic levels using ecological analogies (Gutierrez 1996), they capture the biology of an invasive (or any) species in a detailed mechanistic way using a similar or lower number of parameters than used by phenological models or ENMs. Ecological complexity is tackled at the conceptual level and hence does not result in an increased number of parameters. This also enable explicitly accounting for interactions with other species in a mechanistic way (e.g., Gutierrez et al. 2005Gutierrez et al. , 2017. The PBDM approach addresses the ongoing lack of an effective ) and unified method (Dick et al. 2017) to assess and manage invasive species under climate change (Lehmann et al. 2020); a deficiency that has prevented invasion biology from becoming a more predictive science (Cuthbert et al. 2019). PBDM methods are applicable to potential invasive species of any taxa and can help unify the currently fragmented invasion science landscape (Dick et al. 2017) based on a physiologically-based paradigm of ecological analogies (Gutierrez 1996). The PBDM approach is not limited to terrestrial ecosystems (d'Oultremont and Gutierrez 2002) and can be extended to support ecological risk modeling in aquatic ecosystems (see e.g., Solovjova 2019). By design (Gutierrez and Ponti 2013), PBDM analyses would help address key open questions for devising management strategies for T. absoluta and other invasive species, which rate is increasing across most taxa (Seebens et al. 2017). We must and can do a better job of mitigating biological invasions under global change, but this requires commitment and readjustments of priorities in responsible agencies-it requires investigating the holistic biology of potential invasive species (see Campos et al. 2020). Climate Analytics Group (http://www.climateanalyticsgroup. org) for help in obtaining the NASA climate model weather data. The climate scenario for the USA and Mexico are from the NEX-GDDP dataset, prepared by the Climate Analytics Group and NASA Ames Research Center using the NASA Earth Exchange, and distributed by the NASA Center for Climate Simulation (NCCS). The climate scenario for the Euro-Mediterranean region was developed and provided by the Laboratorio Modellistica Climatica e Impatti at ENEA, Rome, Italy.
Authors' contribution LP and APG conceived and designed the work, performed the PBDM/GIS analysis, and wrote the initial draft of the manuscript. MRC and ND contributed biological data for T. absoluta. MN performed GIS analysis of CLIMEX studies, enabling comparisons of PBDM and CLIMEX maps. All authors contributed to the interpretation of data, substantively revised the manuscript, and have approved the submitted version.
Funding Open access funding provided by Ente per le Nuove Tecnologie, l'Energia e l'Ambiente within the CRUI-CARE Agreement. This research was supported by the Center for the Analysis of Sustainable Agricultural Systems Global (CASAS Global, http://www.casasglobal.org/), by Agenzia nazionale per le nuove tecnologie, l'energia e lo sviluppo economico sostenibile (ENEA), Rome, Italy, by the project MED-GOLD that received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 776467, and by the IPM Innovation Lab (USAID Cooperative Agreement No. AID-OAA-L-15-00001) for funding to N.D.
Data availability All data used in the analysis are publicly available or were sourced from the scientific literature (see Materials and Methods).
Code availability The source code for this PBDM is in Borland Pascal code that is embedded in a larger system base that includes PBDMs for 40 different species of plants, herbivores, parasitoids, predators, and pathogens that have been published as PBDM analyses based on the same methods and implemented in a GIS context (e.g., Gutierrez and Ponti 2014). These models are currently undergoing recoding in Python using an object-oriented programming paradigm for release as open source .The mechanistic PBDM for T. absoluta was developed based on data available in the literature, and is similar to the model developed for the light brown apple moth Epiphyas postvittana (Gutierrez et al. 2010). Like the other PBDMs developed in the last three decades, the Pascal code for the T. absoluta PBDM is not licensed nor it is deposited in a code repository, and is managed by the nonprofit scientific consortium Center for the Analysis of Sustainable Agricultural Systems Global (CASAS Global, http://www. casasglobal.org/). The PBDM algorithms as well as key innovative code for the distributed maturation time models with and without attrition, have been published (Gutierrez 1996, page 157).

Conflict of interest The authors have no conflicts of interest.
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/.