Climate tolerances of Philaenus spumarius should be considered in risk assessment of disease outbreaks related to Xylella fastidiosa

The bacterium Xylella fastidiosa (Xf) is an invasive insect-borne pathogen, which causes lethal diseases to important crops including olives, citrus, almonds and grapes as well as numerous forest, ornamental, and uncultivated plants. Outbreaks of Xf-related plant diseases are currently occurring in the Mediterranean region, causing substantial losses to various agricultural sectors. Several models have recently been published to identify which regions are at highest risk in Europe; however, such models did not consider the insect vectors, which constitute the key driver of short-range Xf spread. We fitted bioclimatic species distribution models to depict the macroclimatic preferences of the meadow spittlebug Philaenus spumarius L. (1978) (Hemiptera: Aphrophoridae), the major epidemiologically relevant vector currently responsible for Xf spread in the Europe. Many regions of Western Europe and Mediterranean basin are predicted by models as highly climatically suitable for this vector, including all regions where severe Xf have occurred so far. Conversely, the driest and warmest areas of the Mediterranean basin are predicted as little suitable for P. spumarius. Models forecast that agricultural-important parts of the southern Mediterranean area might experience a substantial decrease in climatic suitability for P. spumarius by the period 2040–2060. Areas predicted as highly suitable just for the bacterium but not optimal for this vector are apparently still free of severe Xf outbreaks, suggesting that climate tolerances of P. spumarius might partly explain the current spatial pattern of Xf outbreaks in Europe and should always be considered in further risk assessments.


Introduction
The bacterium Xylella fastidiosa (Xf) is an insect-borne plant pathogen that induces severe diseases in a wide array of ornamental and agricultural crops including Pierce's disease of grapevines (PD), olive quick decline syndrome (OQDS), almond leaf scorch (ALS) and citrus variegated chlorosis (CVC) among others (Delbianco et al. 2019). Xylella fastidiosa colonizes the xylem of plants and when the bacterial load is sufficiently large reduces xylem sap circulation. This disruption of water transport along xylem vessels causes decline in yield and fruit quality and ultimately may lead to the death of plants (Hopkins 1989).
Native to Americas, Xf has invaded many regions including Asia, the Middle East and Europe (EFSA PLH Panel., 2019). In 2013, an OQDS outbreak reported in the Apulia region of South-Eastern Italy represented the first confirmed introduction and establishment of Xf in Europe (Saponari et al. 2013). Thereafter, the bacterium has been detected on a diverse range of crops and ornamental plants in France, Spain, Italy, Portugal, Germany and Israel (https:// gd. eppo. int/ taxon/ XYLEFA/ distr ibuti on) (Fig. 1). Although Xf has only recently been officially detected in Europe, epidemiological data indicate that the Xf introduction is likely an older event in various Mediterranean areas (Soubeyrand et al. 2018;Moralejo et al. 2020), suggesting that the pathogen could be widespread in southern Europe and that outbreaks might occur in new areas if favorable ecological conditions are met. Since no efficient treatment exists to cure infected trees, control strategies implemented in the European Union mainly rely on pathogen early detection, removal of susceptible hosts in localities experiencing outbreaks, control of vectors and application of restrictive measures on planting material (EFSA PLH Panel., 2019).
The worldwide economic impact of Xf is massive, e.g. annual losses caused by PD were estimated to reach ca. 104 million $US in California (Tumber et al. 2014) while the costs of removing trees affected by CVC were estimated at  (Bove and Ayres 2007). Recent OQDS outbreaks caused the death of hundreds of thousands of olive trees and induced an unprecedented economic crisis in the Apulian agriculture (Saponari et al. 2019). In addition, the impact of Xf went beyond purely economic losses, as olive cultivation is a centuries-old pillar of Mediterranean landscape, gastronomy and identity. Recently published bioeconomic models forecast that the oliveinfecting Xf subspecies pauca might cause economic losses amounting to several billion euros over the next 50 years in Spain, Greece and Italy (Schneider et al. 2020); however, such models overlooked the pivotal factor mediating either Xf dispersal, or long-term establishment in a new area, i.e. insect vectors presence and abundance.
Growth of Xf in the plant xylem strongly depends on temperature conditions, with cold temperatures substantially reducing bacterial load in infected grapevines ("cold curing" phenomena; Feil & Purcell, 2001). Previous bioclimatic models aimed at predicting climate suitability for Xf in Europe concluded that the Mediterranean countries were the most suitable (EFSA PLH Panel et al. 2019;Godefroid et al. 2019;Schneider et al. 2020). However, these predictions focused solely on the climatic tolerances of the bacterium and incompletely represented Xf epidemiology by omitting the vector ecology and distributions. Accounting only for the environmental tolerances of the bacteria itself is likely not sufficient to assess risk of Xf-related agricultural losses as outbreaks incidence also depends on vector field abundances, host preferences, dispersal patterns and transmission efficiencies (Purcell 1981;Daugherty and Almeida 2009). Transmission of Xf is positively correlated to vector abundance and the time they spend on the host, so severe Xf outbreaks may result from vectors with abundant or dense populations that reside for long periods on a Xfsusceptible plant, despite having a relatively poor propensity (i.e. defined by Irwin & Ruesink (1986) as the probability of a vector transmitting a pathogen under field conditions) (Purcell, 1980).
The meadow spittlebug is a univoltine species overwintering as eggs; eggs hatch in late winter or early spring and nymphal instars develop producing characteristic foams that presumably act as a thermoregulation and defensive strategy (Weaver and King 1954;Tonelli et al. 2018). Adults usually appear in spring and live until late summer or fall, even though few adults of P. spumarius were occasionally sampled in late winter in France and Italy (Bodino et al. 2019;Cunty et al. 2020). In Mediterranean regions characterized by high precipitation seasonality, adults usually spend spring feeding on herbaceous vegetation and migrate to woody hosts-including a large range of trees, shrubs and agricultural crops-in late spring or early-summer concomitantly to rainfall decreases when herbaceous vegetation dries out Cruaud et al. 2018;Morente et al. 2018;Antonatos et al. 2019;Bodino et al. 2020).
The meadow spittlebug is a highly polyphagous species (Ossiannilsson 1981;Cornara et al. 2018) with a wide geographic range, including Europe, Asia, Middle East, North Africa, North America, and New Zealand . All available experimental and distributional data suggest that the meadow spittlebug is sensitive to humidity and temperature, cool and moist conditions being optimal for its local long-term multiplication. For example, an increase in summer temperatures and drought episodes were suggested to have induced a local decline in P. spumarius populations in coastal California in the last decades (Karban and Strauss 2004;Karban and Huntzinger 2018).
The present study aims at depicting the climate tolerances of P. spumarius, the main epidemiologically relevant vector involved in Xf spread in Europe. The ultimate objective of this study is predicting the current and future climate suitability of Europe for this vector. To achieve this goal, we used bioclimatic species distribution models (SDMs) that establish a statistical relationship discriminating presence and absence of a species as a function of environmental conditions (Peterson et al. 2011). These models may serve to predict the potential distribution of a species in new space and time and are thus widely used in invasion risk analyses and assessments of global change effects on biodiversity 1 3 (Peterson et al. 2011). To validate our predictions and highlight the crucial importance of accounting for vectors in further risk assessments, we then explored the link between the SDM-derived climate preferences of P. spumarius and the spatial pattern of Xf outbreaks currently occurring in Europe.

Bioclimatic predictors
We downloaded bioclimatic descriptors from the CHELSA database (Karger et al. 2017) at a resolution of 30 arc-seconds (ca.1 km), which represent average historical climatic conditions during the 1979-2013 period (hereafter referred to as current climate conditions). We downscaled all rasters to a 2.5 arc-minutes resolution (i.e. averaging values among neighbor pixels), to gain computational time. We used three temperature descriptors in further ecological niche modeling, i.e. the mean temperature of the warmest quarter of the year (bio10), the mean temperature of the coldest quarter of the year (bio11) and the average of March, April and May mean maximal temperatures (Tmaxspring). To reflect moisture conditions, we calculated a climatic moisture index (CMI) proposed by Willmott and Feddema (1992) for two periods of the year, i.e. the warmest quarter of the year (cmi_summer) and the coldest 8-month period of the year (cmi_coldest_8months). We built this CMI using the envirem R package (Title and Bemmels 2018). The CMI ranges from − 1 (very arid conditions) to + 1 (very humid conditions). Calculation of CMI was based on CHELSA monthly maps of mean, minimum and maximum temperature and precipitation as well as maps of extraterrestrial solar radiation generated using the palinsol R package (Crucifix 2016). We selected these covariates for further ecological niche modeling because they are presumably key climatic conditions that are crucial for survival of the different life stages of the meadow spittlebug (Weaver and King 1954;Whittaker and Tribe 1996;Karban and Strauss 2004).
We predicted the climate suitability of the Mediterranean area for P. spumarius by the period 2040-2060 using future climate predictions reported in the fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC). To account for uncertainty in available forecasts of global change trends, we used future climates simulations from three global circulation models (GCMs), i.e. the Model for Interdisciplinary Research on Climate version 5 MIROC5 (Watanabe et al. 2011), the Canadian Earth System Model version 2 (CanESM2; Arora et al. 2011) and the HadGEM2-CC (Bellouin et al. 2011). These simulations of future climates were extracted from the CHELSA database, at a resolution of 30 arc-seconds and downscaled to a 2.5 arc-minutes resolution (i.e. averaging values among neighbor pixels). We considered two Representative Concentration Pathways (RCPs), i.e. RCP45 and RCP85 that represent moderate and high future greenhouse gas emissions, respectively (Van Vuuren et al. 2011).

Distribution data
We collected presence records of P. spumarius in scientific literature and biodiversity databases such as the United States Federal Resource for Biological Occurrence Data (BISON-https:// bison. usgs. gov/# home) and the public database Global Biodiversity Information Facility (GBIFhttps:// www. gbif. org/). The reliability of records was assessed and erroneous, ambiguous, or non-reliable records were discarded from this dataset. As we aim at depicting continental-scale macroclimate suitable conditions for establishment of P. spumarius, we also discarded presence records that could presumably be related to the availability of moist microclimatic conditions in particular habitats (e.g. riparian areas, borders of lakes and rivers, irrigated fields or wetlands). As an example, we discarded some museum occurrences related to specimens collected in arid areas of Iran, Turkey, Kyrgyzstan as well as great plains and Rocky Mountains in North America, which likely reflect microclimatic moist conditions available in the vicinity of lakes and rivers. Such occurrence filtering is particularly challenging when microclimatic conditions profoundly affect the distribution of a species as this is the case for P. spumarius. We performed this filtering according to expert-based knowledge on biology and distribution of this species. Besides, we also assembled presence data collected by different Spanish plant protection agencies and research institutions (i.e. "Instituto de Ciencias Agrarias" at CSIC, Madrid, Spain; "Servicio de Sanidad Vegetal de la Junta de Andalucía" based in Sevilla and Jaén, Andalusia, Spain; Sanidad Agrícola Econex S.L. based in Murcia, Spain). As further ecological models were fitted using climatic descriptors relying on average monthly values of precipitation and temperature (see below), we discarded the few presence records from southern Hemisphere distributed in New Zealand. A total of 17,443 presence records were collected (Appendix S1). To avoid using redundant information in models, we only kept one occurrence by cell in the climatic rasters. After discarding these duplicated records, we obtained a final presence dataset amounting to 5,510 records.

Absences and background data for P. spumarius
The nature of absence data is one of the most crucial aspects in SDM (Lobo et al. 2010). While it is widely accepted that true absences (i.e. localities unsuitable for long-term establishment of the species under study) are preferable (Peterson et al. 2011), SDMs are usually fitted using "pseudo-absences" or "background points" (i.e. localities considered as available for the species under study and whose the status of presence/ absence is unknown) when reliable true absences are lacking. In the present study, we randomly generated a total of 20,000 background points equally divided into both North America and in Europe, excluding pixels associated to at least one presence of P. spumarius in a raster at 10 arcminute resolution (Appendix S2).
We recognize that the relative scarcity of P. spumarius in some areas of North America and Europe might be partly explained by biotic interactions (e.g. competition, predation, host distribution) or landscape conditions. Also, as P. spumarius is invasive in the Americas, we recognize that its invaded distribution might not be at climatic equilibrium. We suggest, however, that both these hypotheses are unlikely since P. spumarius is a highly polyphagous species, which is therefore little influenced by host availability, and whose distribution in North America has remained relatively stable for many decades (Weaver and King 1954). We assumed that the current distribution of the meadow spittlebug in both Europe and the Americas likely reflects its macroclimate tolerances. Besides, we cannot discard the hypothesis that a genetic bottleneck occurred during the invasion process and has made invading populations unable to fill the fundamental niche of this species (Allendorf and Lundquist 2003).
We also compiled a set of field-collected absences of P. spumarius recorded by Spanish plant security agencies in southern regions of Spain. The set of absences was retrieved from four years (2015-2018) of thorough and opportunistic field sampling performed in south-western Spain (i.e. in the Andalusian provinces of Sevilla, Cadiz, Malaga, and Huelva). This sampling consisted of visiting a total of 516 localities in various habitats (i.e. mostly natural habitats, roadsides or olive groves with natural cover vegetation presumably suitable for P. spumarius establishment). The sampling was focused on P. spumarius distribution and aimed at visiting different kinds of habitats and climatic regions. A sampling point consisted in searching nymphs and/or adults in vegetation during ca. 15 min in a diameter of ca. 100 m. Sampling of adults was performed using a sweeping net. Sampling was usually performed by two or three samplers with high experience and expert knowledge on spittlebugs ecology. Samplers followed main roads by car and performed sampling in areas with environmental conditions the most suitable as possible (suitable cover vegetation, good moisture conditions, suitable landscape, etc.). More details on this sampling survey are provided by Durán (2018). We assembled this absence dataset with absence records collected during the 2016-2018 period by Morente et al. (2018) in Andalusia autonomous community, Spain. To minimize the risk of using false absences, the field-collected absences situated in the same pixel as any available presence record on a bioclimatic raster at 10 arc-min resolution were discarded from this dataset. We purposely used a coarser resolution raster to filter absences to, as much as possible, avoid including false absences in the dataset. This filtered set of absences encompassed 103 records. We combined these absences with the presences collected during the same field surveys to build a spatially and temporally independent evaluation dataset that will be used for evaluation of models (referred as to "Pres-Abs-Spain1"; Appendix S3).
We warrant that these absences might also encompass climatically suitable localities and could partly reflect unsuitable landscape conditions as olive cropping in southern Spain is particularly intensive (i.e. frequent soil tillage, weeds removal and application of chemical inputs, which likely negatively impacts populations of spittlebugs). However, given the high numbers of sampled localities, the experience of samplers and available knowledge on the physiology of P. spumarius, we assume that these absence data likely reflect a climate-driven distribution pattern in southern Spain. This assumption was strongly supported by preliminary analyses showing that most of bioclimatic models fitted using only distribution data available from North America predicted low climate suitability in these absence localities (see results and discussion parts).

Abundance data
An abundance dataset was collected in the Murcia autonomous community, Spain, during spring 2019. Spittlebugs adults were collected using a sweeping net in 20 agricultural fields (i.e. 4 olive groves, 13 almond groves and 3 vineyards), twice a month from April to August 2019 (i.e. a total of 8 sampling dates for each field; Appendix S4). As natural vegetation was seldom in the core of most of fields because of soil tillage, sampling was also performed in their close surroundings. At each sampling date, a hundred sweeps were performed in each of three components of landscape, i.e. cover herbaceous vegetation, the canopies of cultivated plants and the canopies of surrounding trees and shrubs. We summed the total spittlebug adults collected in each field and its surroundings (referred to as "Murcia-abundance" dataset). This abundance dataset is another available spatially and temporally independent evaluation dataset that will be used for evaluation of models.

Modeling approaches
Many SDM algorithms exist and may yield variable forecasts, warranting the usefulness of inspecting outputs of different techniques (Araújo and New 2007). We modeled the climatic niche of P. spumarius using four different modeling approaches, i.e. boosted regression trees (BRT; Elith et al. 2008), generalized additive models (GAM; Guisan et al. 2002), multiple adaptive regression splines (MARS; Leathwick et al. 2005) and Maxent (Phillips et al. 2006). When starting this study, we also fit random forest models and generalized linear models. However, these models were discarded as they repeatedly yielded highly non-ecologically realistic outputs (data not shown).
Models were fitted using the "biomod2" (Thuiller et al. 2016) package in R (R core Team 2013). The BRT models were fitted using default settings implemented in biomod2 package. For the GAM and MARS models, we did not account for interactions between considered variables to ensure simplicity and interpretability of models. In MARS, the "polynomial" model type was fitted (Thuiller et al. 2016). Maxent can theoretically use six features classes (FCs), namely linear (L), quadratic (Q), product (P), threshold (T), categorical (C) and hinge (H) (Phillips et al. 2006). To avoid as much as possible model overparameterization and the fitting of too complex and nonecologically realistic responses to climate, we evaluated best-predicting parameters among different combinations of regularization multipliers (RM; values ranging from 0.5 to 4 with an increment of 0.5) and FC combinations ("L", "LQ", "H", "QH"; Radosavljevic and Anderson 2014). We selected the combination of RM and FC associated to the highest mean AUC.
In model calibration, we equaled the weighted sum of presences and the weighted sum of absences (Barbet-Massin et al. 2012). The predictive power of models was evaluated using the area under the curve (AUC) of the receiving operator curve (ROC; Fielding & Bell, 1997) and the True Skill Statistic (TSS; Allouche et al., 2006). We fitted models using 80 percent of data whereas the remaining 20 percent of records were kept for model evaluation. We also evaluated models using the spatially and temporally independent dataset "Pres-Abs-Spain1".
To account for uncertainty induced by cross-validation methods, we performed five replicates of every model. For each repetition, we processed a consensus prediction of the models, which is the average predicted climatic suitability across the different modeling methods. To reduce uncertainty in projections, we combined predictions from models with AUC > 0.7 when processing the consensus. We weighted individual model contributions to the consensus predictions by their AUC scores. We evaluated the predictive power of each consensus prediction using the same part of the data reserved for model evaluation.

Bioclimatic predictors selection
The selection of the environmental dataset used in SDM has a crucial impact on the accuracy of potential geographic range predictions (Pliscoff et al. 2014). It is widely recognized that using ecologically relevant climate descriptors that have a direct impact on the study species, enhances model transferability to novel space and time (Peterson et al. 2011). Also, as ecological data are often autocorrelated, a traditional random data partitioning into calibration/evaluation datasets may artificially inflate accuracy metrics. It is thus usually preferable to estimate model transferability using an evaluation dataset that is spatially or temporally independent from the calibration domain (Randin et al. 2006). Before fitting definitive full models, we thus applied a preliminary intercontinental transferability evaluation approach to select the best-predictive environmental dataset used in further definitive full models. This approach consists of fitting models using data from a region and evaluating their predictive power using evaluation data available for another geographic area (Randin et al. 2006). This approach therefore identifies proximal descriptors that have a direct causality to the biology of the study species and avoids potential issues concerning differences in predictor correlation structures among different geographic areas or time periods (Jiménez-Valverde et al. 2009). This framework is particularly suited for a widely distributed polyphagous invasive species such as P. spumarius whose geographic distribution is wide, well-documented, presumably at equilibrium (i.e. the species has likely reached most of areas climatically suitable in the biogeographic regions where it occurs) and minimally constrained by biotic factors (e.g. plant host distribution).
We first visually compared the distribution of the values for the five bioclimatic variables extracted at presence records among North America and Eurasia. Our assumption is that ecologically relevant descriptors should yield similar limits and distributions in both geographic areas. We then fitted models using only occurrences and background points of P. spumarius situated in the Americas and calculated the average predictive accuracy of consensus predictions (AUC and TSS) on the remaining presences available for the rest of the world and background data generated in Europe ("eval1"). We also calculated the average predictive accuracy of the consensus predictions (AUC and TSS) on the "Pres-Abs-Spain1" dataset ("eval2"). To select the climate dataset used in further definitive full models, we checked for ecological reliability of the modeled response curves to climate and accuracy metrics of models on these two evaluations datasets.
Applying this approach, we tested nine climate datasets representing different combinations of few bioclimatic descriptors (Appendix S5). At purpose, we fitted models using few bioclimatic variables to avoid model overfitting as much as possible (Randin et al. 2006). We adopted a such conservative approach according to the purpose of the study, i.e. minimizing omission error and transferring model to new time and space. Since correlation between environmental descriptors may have a negative impact on model performance, we ensured that the climate datasets did not encompass highly correlated variables (Pearson's correlation index < 0.8 considering climatic values extracted at all presence and background localities).

Full model calibration and evaluation
Sampling bias is a common challenge for ecological modelers. Sampling bias occurs when some regions are oversampled when compared to others, which may lead to models reflecting sampling effort rather than the true physiological preferences of the considered species. Even though accounting for sampling bias has been suggested to lead to better predictions in some situations (Varela et al. 2014), this cannot be done blindly as high densities of presence records in a particular region of the climatic space may simply reflect higher climatic suitability and therefore higher detectability of the species under study. When fitting preliminary full models with all occurrences/background data, we observed that some models yielded bimodal response curves to climates slightly unrealistic according to niche theory (Austin 2007). Our hypothesis was that these slightly unrealistic bimodal modeled response curves could come from a too high spatial clustering of occurrence in some regions of Europe or North America. To attempt at improving ecological reliability of outputs, we slightly reduced spatial aggregation of presence records by keeping only one occurrence in a radius of 10 km. While we are aware that some information can be lost by reducing the number of records used in models, we obtained more biologically realistic unimodal responses curves to climate descriptors after this slight thinning of the presence dataset.
We calibrated full models using the previously selected best set of bioclimatic predictors (i.e. tmaxspring and cmi_coldest_8months; see results below) and applying the previously described four modeling methods. As above-described, we evaluated models using 20% of data kept for evaluation as well as the spatially and temporally independent "Pres-Abs-Spain1" evaluation dataset. We also fit a generalized linear model (quasi-Poisson family) to address the impact of the log-transformed average predicted climatic suitability on the abundance of P. spumarius in south-Eastern Spain ("Murcia-abundance" dataset).

Selection of bioclimatic descriptors
After analyzing results of our preliminary analyses, we concluded that the most ecologically realistic estimates of the potential distribution of P. spumarius were provided by models fitted using both covariates Tmaxspring and cmi_coldest_8months. The selection of this climate dataset was motivated by three main observations: (1) Density curves of climate conditions extracted at presence records indicate that Eurasian and North American geographic ranges of P. spumarius have similar limits and distributions for the covariates Tmaxspring, cmi_summer and cmi_coldest 8 months (Appendix S6). Conversely, Eurasian and North American ranges of P. spumarius yield much different limits and distributions for both bio10 and bio11 (Appendix S6). This is not surprising as both covariates Tmaxspring and cmi_coldest 8 months reflect temperature and moisture conditions during the development of eggs and nymphs, which are very sensitive to desiccation and high temperatures (Weaver and King 1954;Halkka et al. 2006). (2) Both covariates were ranked by models among the best performing climate covariates to predict the European range of P. spumarius when fitting models using only data available from America (Appendix S5). The average consensus provided by models fitted with only American occurrences and both covariates Tmaxspring and cmi_coldest_8months yield good evaluation metrics (mean AUC ± standard deviation = 0.91 ± 0.01; mean TSS ± standard deviation = 0.72 ± 0.03) and accurately predicted the distribution of P. spumarius in Europe (mean AUC ± standard deviation = 0.78 ± 0.003; mean TSS ± standard deviation = 0.38 ± 0.04) as well as the presence/absence dataset collected by Spanish plant protections agencies (mean AUC ± standard deviation = 0.79 ± 0.002; mean TSS ± standard deviation = 0.44 ± 0.05) (Appendix S5). Moreover, the positive relationship between the logtransformed average consensus climate suitability and abundance of P. spumarius adults in Murcia autonomous community, was highly significant according to a GLM (P-value < 0.05). This suggests that this climate dataset constitutes a good proxy to estimate the range of P. spumarius.

Predictions of definitive full models
We observed slight differences in outputs among the different algorithms. We suggest that selection of the modeling approach must be partly selected according to biological reliability of modeled response curves to climate descriptors rather than only a combination of evaluation metrics, which can be highly sensitive to bias in calibration and evaluation data. After having meticulously analyzed the modeled response curves to climate and evaluation metrics, we suggest that the algorithm GAM was the approach yielding the most realistic outputs. On the one hand, the other models were also promising but showed some slight weaknesses, e.g. BRT and Maxent infer non-null climate suitability in regions with very high Tmaxspring, superior to the maximum value of Tmaxspring experienced by P. spumarius (Appendix S8). Similarly, we suggest that the decline in climatic suitability inferred by MARS at high values of Tmaxspring might be too abrupt to be extremely realistic (Appendix S8). On the other hand, GAMs yielded continuous responses to moisture and relatively regular bell-shaped responses to Tmaxspring, which is more realistic from a biological perspective (Appendix S8). In the following manuscript, we therefore only will present predictions and discuss outputs related to the GAM models. The GAM models yielded high predictive accuracy (mean AUC ± standard deviation = 0.85 ± 0.01; mean TSS ± standard deviation = 0.59 ± 0.01) indicating excellent model discrimination between presences and background points. They accurately predicted the presence/absence dataset collected by Spanish plant protections agencies (mean AUC ± standard deviation = 0.77 ± 0.01; mean TSS ± standard deviation = 0.11 ± 0.01). Similarly, the positive relationship between the log-transformed average GAM-derived climate suitability and abundance of P. spumarius adults in Murcia autonomous community, was highly significant according to a GLM (P-value < 0.05).
According to GAM models, the Mediterranean areas where severe Xf outbreaks currently occur (i.e. Balearics islands, Alicante region in Valencia autonomous community, north of Portugal, lowlands in Corsica, southeastern France, Apulia and Tuscany regions in Italy) yield high climatic suitability for P. spumarius (Fig. 2). Models predict that large parts of colder and humid regions of Europe (e.g. Atlantic and Mediterranean coastal regions of France, northern and coastal regions of Portugal, northern and high-altitude parts of Spain, Northwestern Italy, southern parts of Apulia, Western Greece and Northwestern coasts of Turkey) are currently highly climatically suitable for P. spumarius (Fig. 2). Lower climate suitability was, however, predicted in arid and warm regions from the Iberian Peninsula, eastern Greece, northern parts of Apulia region, and central arid mountainous regions of Turkey (Fig. 2).
Bioclimatic models predict that most of the Mediterranean region will experience a decrease in climatic suitability for P. spumarius by the period 2040-2060 (Fig. 3). Predictions are overall similar among the different GCMs used in   in the entire Mediterranean area. This prediction reflects the mean predicted climate suitability across five generalized additive models replicates the study (Fig. 3). Geographic areas that will experience the most important decrease in climatic suitability for P. spumarius by the period 2040-2060 encompass central and southern Iberian Peninsula (particularly Portugal), coastal regions of Turkey, central and southern France and most of lowlands in Italy (Fig. 3).

Climate tolerances of the meadow spittlebug should be considered when assessing Xf outbreak risk
The availability of accurate forecasts of the current and future distribution of pests is a crucial prerequisite for the implementation of monitoring and control strategies in a rapidly changing world (Venette et al. 2010). Many harmful pests are currently experiencing rapid distribution shifts due to ongoing climate change or human activities (Parmesan 2006) and the plant pathogen Xf is likely not an exception to this rule. Several studies have recently addressed the environmental tolerances of the bacterium Xf with the aim to predict which regions are at high risk in Europe (Godefroid et al., 2019;EFSA PLH Panel, 2019;Schneider et al., 2020). Moreover, some bioeconomic models have predicted that the ongoing climate change will likely promote the emergence of Xf in crucial Mediterranean agricultural areas with dramatic socio-economic consequences (Schneider et al. 2020). However, all these studies overlooked the pivotal factor mediating Xf spread, i.e. the distribution and abundance of insect vectors (Purcell 1981;Daugherty and Almeida 2009), and consequently tend to overestimate the risk in regions where efficient vectors are rare or absent. Indeed, according to Purcell (1981), the frequency and intensity of Xf outbreaks are positively correlated to abundance of efficient vectors in the field, the time they spend on hosts and their transmission efficiency.
In the present study, we modeled the macroclimatic tolerances of P. spumarius-the only epidemiologically relevant Xf vector described in Europe (Cornara et al. 2017(Cornara et al. , 2020Cruaud et al. 2018;Cavalieri et al. 2019). We focused on this species because the relative contribution of other putative vector species in Xf epidemiology in Europe is still unclear. The models presented here predict a marked spatial variation in climatic suitability throughout the European continent for the meadow spittlebug. Overall, bioclimatic models predict high climatic suitability in most of Western Europe while the warmest and driest Mediterranean regions-mainly situated in Iberian Peninsula, eastern Greece, and central Turkey-are predicted as moderately or little suitable (Fig. 2). This finding is crucial since some of the geographic regions predicted as little suitable for P. spumarius encompass economically major agricultural areas predicted at high risk for long-term establishment of the bacterium itself (e.g. Southern Spain, which is the world's largest producer of olive; Godefroid et al., 2019;EFSA PLH Panel, 2019;Schneider et al., 2020). In addition, the present study suggests that large parts of southern Europe including areas where severe Xf outbreaks currently occur (e.g. Apulia province in Italy and Alicante region in Spain)-could experience a decrease in climatic suitability for P. spumarius by the period 2040-2060 (Fig. 3). These forecasts of current and future climate suitability will have crucial practical applications when assessing both the economic impact of Xf and the implementation of control strategies. For example, bioeconomic models estimated massive putative economic losses due to spread of Xf in the European olive-producing regions over 50 years (Schneider et al. 2020); however, the present study mitigates these conclusions and calls for a complete reassessment of potential Xf impact to European agriculture under current and future climate conditions. Noticeably, we found that Xf outbreaks have occurred so far only in Mediterranean areas predicted as highly suitable for P. spumarius by the models presented here (i.e. Balearic Islands, Alicante region, north of Portugal, lowlands in Corsica, southeastern France, Tuscany and Apulia, Italy; Fig. 2 -Appendix S8). Conversely, no severe Xf outbreak has ever occurred in Mediterranean areas predicted as little suitable for P. spumarius, even though some of these regions (e.g. most of southern Iberian Peninsula) were predicted as suitable for the establishment of different Xf bacterium subspecies (Godefroid et al. 2019;Schneider et al. 2020). One hypothesis to explain this pattern is that climate conditions preclude the establishment of local abundant populations of efficient vectors in these regions and consequently, drastically reduce probability of Xf transmission to plants and subsequent spread of the pathogen. To properly verify this hypothesis, we should, however, ideally assess that Xf was introduced in these regions in the past and did not spread because of the lack of abundant vector populations. Although it is likely that the bacterium Xf was long ago introduced in most of agricultural regions of Europe (Soubeyrand et al. 2018;Moralejo et al. 2020), further research should be needed to support this hypothesis. This result emphasizes the importance of accounting for vectors' ecological characteristics when assessing Xf outbreak risk and highlights the risk of gathering unreliable predictions by excluding vectors when modeling the spread of a vectorborne plant pathogen.

Accuracy of models and predictions
Models were purposely calibrated using few ecologically relevant climate descriptors, which is a conservative approach to avoid overfitting, reduce omission risk and enhance model transferability. We meticulously inspected the modeled species' response curves to bioclimatic covariates and tested the transcontinental transferability of models to ensure ecological reliability and transferability of climate suitability forecasts. This modeling effort highlighted three main P. spumarius-climate relationships: (1) High winter temperatures were not identified as a main climatic constraint acting on the distribution of P. spumarius in the Mediterranean region and North America. When starting this work, one hypothesis was that the absence or scarcity of P. spumarius in South-Eastern USA, inlands of California and some regions of southern Europe might partly be explained by failure in eggs development induced by too high winter temperatures (Weaver and King 1954). Distributional data and the transcontinental transferability approach suggest, however, that this hypothesis is unlikely. Indeed, high densities of P. spumarius were encountered in European regions characterized by high winter temperatures like those observed in South-Eastern USA where this species is absent. Also, preliminary models calibrated with both winter temperatures and spring temperatures as covariates (data not shown), did not identify high winter temperatures encountered in the region of study as a major constraint for P. spumarius presence.
(2) Some presence records of P. spumarius were associated with Mediterranean localities characterized by extremely hot summer temperatures like those observed in the regions of South-Eastern USA where this species is absent. Conversely, we found a good match between the superior limit of average maximal spring temperature experienced by P. spumarius in both Eurasia and America. Therefore, spring temperatures were identified as a better proxy than summer temperatures to predict the distribution of P. spumarius in the Mediterranean region. Summer humidity was also identified as a poor proxy to predict the distribution of P. spumarius in the Mediterranean region. This is not surprising as high densities of P. spumarius can be found in many European regions characterized by various summer rainfall regimes ranging from moist to extremely dry. Adults of P. spumarius usually switch to trees and shrubs in early summer in the Mediterranean region (Morente et al. 2018;Bodino et al. 2020), which might act as an efficient strategy to cope with warm and dry summer conditions (3) Models identified spring temperature and the CMI during the coldest 8-month period of the year as good predictors to forecast the geographic range of P. spumarius in the Mediterranean region. This is consistent with existing knowledge on the physiology of the meadow spittlebug (Weaver and King 1954;Cornara et al. 2018) as both climatic covariates reflect temperature and moisture conditions occurring during the development of eggs and nymphs, which are life stages particularly sensitive to variation in humidity and temperature (Weaver and King 1954). Weaver and King (1954) found higher egg hatching rates at higher levels of humidity under controlled conditions. Halkka et al. (2006) found that nymphal survival rate in P. spumarius is partly explained by temperature and humidity levels of spring. The negative effect of low moisture during the coldest period of the year on P. spumarius might be either direct, by causing desiccation of eggs or nymphs, or indirect by precluding establishment of an abundant cover vegetation, which provides moist microhabitats for eggs and nymphs.
Overall, bioclimatic models yielded excellent accuracy metrics and the predictions ideally fit available data on the distribution pattern of P. spumarius in Europe. For example, models predict low climatic suitability in warm and dry areas of the Iberian Peninsula, such as the Guadalquivir basin, the lowlands of Murcia autonomous community, the olive-producing regions in Andalucía autonomous community and southern parts of Portugal (Fig. 2), which is congruent with field observations (Drosopoulos and Quartau 2002;Durán 2018;Morente et al. 2018;Popova 2020;Torres 2020). Conversely, models predict high climate suitability in humid and cooler places of Spain (e.g. Balearic islands, Alicante region, the "Sierra of Aracena", Huelva), which are well-known hotspots of P. spumarius (data not shown). In addition, a highly significant statistical relationship was found between the number of P. spumarius adults collected in agricultural fields of South-Eastern Spain (Murcia autonomous community) and model-predicted climatic suitability. We are aware that the scarcity of P. spumarius in dry and warm lowlands of southern Spain might be partly explained by additional factors such as landscape structure, intensive agricultural practices and human disturbances. The present study suggests, however, that climate conditions likely contribute to shaping this pattern since bioclimatic models fitted using only distributional data available from North America predicted a similar climatic suitability pattern in Spain (Appendix S5 & S8). Besides, model predictions are also congruent with field observations made in Italy, e.g. densities of P. spumarius are positively correlated to altitude in central Italy (Santoiemma et al. 2019). Finally, preliminary field data suggest that densities of P. spumarius are usually higher in cool and moist areas of Corsica island, France (Chartois et al. 2019), which is congruent with the climate suitability maps presented here.
We urge, however, caution when interpreting SDM outputs since such correlative approach only depicts the realized climatic niche of species, i.e. a subset of their full environmental tolerances constrained by biotic interactions and dispersal constraints (Hutchinson 1957). One of the main assumption of models is that P. spumarius distribution is at climatic equilibrium in the biogeographic regions where background points were generated (i.e. North America and Europe). We suggest this is a reasonable assumption as the meadow spittlebug is a widespread and highly polyphagous organism whose absence or scarcity in some areas is congruent with available knowledge on its physiology (Weaver and King 1954;Karban and Strauss 2004). We suggest that calibration of mechanistic models based on physiological data is a key future research path to provide better estimates of the potential climate suitability of Europe for this vector. Finally, we remind that the models fitted here are only based on bioclimatic predictors and thus just provide estimates of continental-scale climatic suitability. The distribution of spittlebugs is shaped by many additional factors that were not accounted for in the present study, e.g. biotic interactions, dispersal range, plant-insect interactions, microclimatic conditions, insects' feeding behavior, cropping practices and landscape structure.

Conclusions
The present study provides estimates of the potential suitability of the Mediterranean basin under current and future climate conditions for the main epidemiologically relevant vector of Xf. Bioclimatic models highlight that some of the most economically important agricultural regions situated in southern parts of the Mediterranean area are not climatically optimal for the establishment of the main vector of Xf. This finding might explain why Xf outbreaks have never occurred in these regions despite these regions being predicted as suitable for multiplication of the bacterium in plants. Importantly, models predict that a substantial part of the Mediterranean regions might become less suitable for the main Xf vector by the period 2040-2060 due to ongoing climate change. The maps provided here are easily interpretable tools that provide valuable information to design monitoring strategies and assess Xf-related risk in Europe.

Author contributions
MG performed statistical analyses with the help of TS; all authors but TS participated in data collection: MG led the manuscript writing; all authors contributed to manuscript writing.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10340-021-01413-z. Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work has been financially supported by the Spanish Ministerio de Ciencia, Innovación y Universidades, Grant/Award Number: AGL2017-89604-R and the European Union Horizon 2020 research and innovation program under grant agreements no. 635646 POnTE (Pest Organisms Threatening Europe), and no. 727987 XF-ACTORS (Xylella Fastidiosa Active Containment Through a multidisciplinary-Oriented Research Strategy). The first author of this study was funded by the fellowship "Ayudas destinadas a la atracción de talento investigador de la Comunidad de Madrid". Daniele Cornara participation in this work was supported by a research grant in the frame of European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 835732 XYL-SPIT.

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