The Potential Global Climate Suitability of Kiwifruit Bacterial Canker Disease (Pseudomonas syringae pv. actinidiae (Psa)) Using Three Modelling Approaches: CLIMEX, Maxent and Multimodel Framework

In recent years, outbreaks of kiwifruit bacterial canker (Pseudomonas syringae pv. actinidiae, Psa) have caused huge economic losses to two major global kiwifruit producers, Italy and New Zealand. To evaluate the potential global risk areas of Psa, three modelling methods (MaxEnt, CLIMEX and a Multi-Model Framework, including Support Vector Machine or SVM) were used. Current global occurrence data for Psa were collected from different sources. The long-term climate data were sourced from WorldClim and CliMond websites. The model results were combined into a consensus model to identify the hotspots. The consensus model highlighted the areas where two or three models agreed on climate suitability for Psa. All three models agreed with respect to the climate suitability of areas where Psa is currently present and identified novel areas where Psa has not established yet. The SVM model predicted large areas in Central Asia, Australia, and Europe as more highly suitable compared to MaxEnt and CLIMEX. Annual mean temperature and annual precipitation contributed most to the MaxEnt prediction. Both MaxEnt and CLIMEX showed the probability of Psa establishment increased above 5 ◦C and decreased above 20 ◦C. The annual precipitation response curve showed that excessive rain (>1200 mm/y) constrains Psa establishment. Our modelling results will provide useful information for Psa management by highlighting the climatically susceptible areas where Psa has not established, such as the USA, Iran, Denmark, Belgium and especially South Africa, where kiwifruit has been planted commercially in recent years.


Introduction
The causal agent of bacterial canker of kiwifruit, Pseudomonas syringae pv. actinidiae (Psa) is a Gram-negative, obligate aerobic bacterium. Symptoms of the disease mainly appear at the beginning of spring when leaf spots become visible, appearing brown surrounded by a bright chlorotic halo [1]. Floral buds may turn brown and exude gum without opening. Translucent exudate on otherwise healthy tissues, such as canes, is also considered to be one of the main symptoms of Psa on kiwifruit (Figure 1a-d) [1][2][3]. Six biovars of Psa have been described worldwide, based on genomic fing and toxin production (biovars 1-6), though some studies have revealed genetic within Psa populations of up to 14 groups [21][22][23]. However, biovar 4 (Psa-Lv transferred to a new pathovar, Pseudomonas syringae pv. actinidifoliorum (Pfm) biovars 1 and 2 are reported to cause moderate damage, biovar 3 is the most vir and has inflicted significant damage in Italy and New Zealand [25]. For exam the early reports in 2013 stated that 77% of New Zealand's kiwifruit orchard fected with Psa [26], the most recent report indicates that this has increased up Psa can infect Actinidia chinensis var. deliciosa (green kiwifruit) and Actinid var. chinensis (gold kiwifruit), kiwiberry (Actinidia arguta) and wild Actinidia kol However, its survival on nonhost plants such as Alternanthera philoxeroid phyllales: Amaranthaceae), Setaria viridis (Cyperales: Poaceae) and Paulow (Scrophulariales: Scrophulariaceae) has also been documented [29,30]. In add neste et al. (18) reported that Psa can survive on Cryptomeria japonica (Japan Psa was first isolated from kiwifruit in Japan in 1984 and described as a new pathovar in 1989 [4]. In the late 1980s, Psa was reported from Korea [5] and from China in 1990 [6].
Six biovars of Psa have been described worldwide, based on genomic fingerprinting and toxin production (biovars 1-6), though some studies have revealed genetic variability within Psa populations of up to 14 groups [21][22][23]. However, biovar 4 (Psa-Lv) was later transferred to a new pathovar, Pseudomonas syringae pv. actinidifoliorum (Pfm) [24]. While biovars 1 and 2 are reported to cause moderate damage, biovar 3 is the most virulent form, and has inflicted significant damage in Italy and New Zealand [25]. For example, while the early reports in 2013 stated that 77% of New Zealand's kiwifruit orchards were infected with Psa [26], the most recent report indicates that this has increased up to 92% [27].
Psa can infect Actinidia chinensis var. deliciosa (green kiwifruit) and Actinidia chinensis var. chinensis (gold kiwifruit), kiwiberry (Actinidia arguta) and wild Actinidia kolomikta [28]. However, its survival on nonhost plants such as Alternanthera philoxeroides (Caryophyllales: Amaranthaceae), Setaria viridis (Cyperales: Poaceae) and Paulownia fortunei (Scrophulariales: Scrophulariaceae) has also been documented [29,30]. In addition, Vanneste et al. (18) reported that Psa can survive on Cryptomeria japonica (Japanese cedar). However, this was achieved by artificial inoculation and no evidence of multiplication was found (Psa populations on C. japonica decreased with time). Until the devastating Italian outbreak of 2008, in which Psa losses were estimated to be €2 million by 2009, it was assumed that the disease's economic importance was relatively low [31,32]. In New Zealand a 2012 publication estimated that the cost of Psa (as Psa-V) would be NZ$310-$410 million in the next five years, and NZ$740-$885 million over the next 15 years [33]. However, by 2014, only two years after these estimates, the cost of export losses alone mounted up to NZ$930 million [34]. Although, estimation of the economic impacts of Psa, or any other plant disease is beneficial, the evaluation of socio-economic and environmental impacts which can be devastating are often rare or ignored. For instance, in New Zealand an early estimation in 2012 claimed that Psa could result in up to 470 full-time job losses per year for the following four years. Such impacts can lead to anxiety and stress which need proper evaluations [33].
Environmental factors such as strong wind and rain are major vectors for introduction into new regions [9,35]. In addition to human-mediated spread though grafting materials, nursery materials and pollen, Psa has been shown to be carried by wind and it has also been shown that frequent spring hailstorms cause wounds that contribute to disease spread over large areas [1,25]. Moreover, the role of some sucking insects such as Metcalfa pruinosa (Hemiptera: Flatidae) as vectors of Psa in laboratory conditions has been investigated and confirmed [36].
Information about the epidemiology of the pathogen is limited. It has been reported that at temperatures between 10 • C and 20 • C, the pathogen is very active [37]. In other studies, the optimum temperature range for growth on new canes is estimated at 12 • C to 20 • C and 15 • C for leaf infection [6,38,39]. While studies in Japan and Korea reported very little symptom development at temperatures >20 • C during a 10-day period, hightemperature tolerance is generally estimated at 25 • C in France, Italy, and Portugal [6,25,40]. In previous studies, rain has been implicated as the most important factor in establishment and distribution of Psa but no detailed quantitative data are available as yet [41,42].
Despite the huge economic losses attributed to Psa in recent years, there are limited studies that evaluate the risk of its establishment in areas where kiwifruit is commercially grown and Psa has not been reported yet. Therefore, studies that could project the habitat suitability of Psa in noninfected areas can be very valuable. Reynaud et al. (2011, cited in [40]) used CLIMEX software package (version 3.0, CSIRO Publishing, Melbourne, Australia) to evaluate habitat suitability based on cold and dry stress parameters. In a preliminary study with limited data, two models (MaxEnt and CLIMEX) were used to evaluate the potential distribution of Psa, which both showed differences in projections in areas of interest [43]. In a study by Wang et al. [44], MaxEnt was used to model the current and future distribution of Psa in China based on different emission scenarios. They identified that environmental variables such as maximum April temperature, mean temperature of the coldest quarter, precipitation in May and minimum temperature in October contributed the most to the model projection. In a recent study focused on China, [45], an ensemble approach using three correlative models including generalized boosting models (GBM), random forests (RF) and classification tree analysis (CTA) was employed to predict the potential distribution of Psa in China under four general circulation models (GCMs). In all three models, precipitation in March contributed the most to model projection (28.2%) followed by the mean temperature of the warmest quarter (17.7%). Do et al. [46] developed a model to evaluate the accumulated potential damage of Psa on kiwifruit during the growing and overwintering seasons which can be set to run using hourly or daily mean air temperature. While their study was only focused on the green kiwifruit cultivar Hayward (A. chinensis var. deliciosa), they noticed that necrotic lesion length increased with temperature, with the longest observed at 35 • C. Beresford et al. [39] developed a weather-based mechanistic risk model to predict Psa development based on a multiplication and dispersal concept. Their model calculated the bacterial multiplication (M index) from a temperature function and the daily risk factor (R index), which is the three-day accumulation of the M index on days with total rainfall > 1 mm/day. They found that the optimum temperature for leaf infection was 21 • C and 15 • C in vitro and in planta, respectively. Kim and Koh [47] adapted the same risk model to predict Psa epidemics in Korea. Due to significant differences between the New Zealand and Korean climates, the temperature and rainfall functions of the model were modified. They found that the model was highly sensitive to rainfall, and disease incidence was high when rainfall was recorded in consecutive days. By testing different ranges of average temperatures, 24 • C was identified as the optimal threshold, with no change in disease incidence when daily maximum temperatures exceeded 27-30 • C. Kim and Koh [48] also developed an integrated modelling approach to investigate the potential for epidemics of Psa under climate change. The integrated modelling approach used three models. First, a climate suitability model was developed which predicted the suitable areas for growing kiwifruit. Then a chill-day model simulated the periods of kiwifruit flowering limited to the suitable areas predicted by the climate suitability model mentioned earlier.
Finally, these flowering dates were fed into Pss-KBB model [49] to simulate the potential risk of the epidemic.
Most studies predicting the distribution of Psa rely on a single model, or use the same type of models (e.g., only correlative models or only mechanistic models). Their predictions for some areas are sometimes different or contradictory which can add to the prediction uncertainty. It is therefore important to develop a consensus model which relies on the agreement of different modelling approaches to investigate the risk of Psa globally. In this way, effective readiness and response measures can be employed to minimize the potentially devastating impacts of Psa at locations of interest.
The objectives of the current study are to: (i) model the potential global distribution of kiwifruit bacterial canker (Psa) using ensemble models and model agreement, (ii) investigate agreement among these models' projections to identify the hotspots and estimate how they contribute to a better knowledge of Psa spread and risk; and (iii) characterize the potential environmental variables important in the distribution of Psa.
Three different types of species distribution models (SDMs) were used to address the objectives: The current distribution records of Psa and long-term climate data were used to project the potential distribution using a semimechanistic model (CLIMEX), a presence-background model (MaxEnt) and multimodel framework (MMF). These models have frequently been used in different studies to project the habitat suitability of different invasive species into new areas [50][51][52].

Psa Presence Data
At first, 324 presence points of Psa were collected from available studies and online sources [10,11,13,31,35,44]. In cases where coordinates were not available and only the name of localities of Psa presence were mentioned in literature, ArcGIS Pro v 2.7.3 and Google Earth Pro v. 7.3.4 were used to find the coordinates. Duplicate points were removed based on the resolution of the data (10 arc minute) to minimize the environmental bias and to avoid possible spatially autocorrelated presence points ( Figure 2 and Table S1). Finally, 201 points were used in the modelling process. Because there is no reliable information about the environmental needs for infection of kiwifruit vines by different strains or biovars, these strains were not differentiated among the presence records in this study. The records of biovar 4 were not included in this study as they are now proven to belong to a different pathovar [24].

Bioclimatic Variables
Climate data files (Met manager version 1.2) from CliMond website were downloaded for direct use in CLIMEX software (www.climond.org) [53]. Bioclimatic variables [54] that have been frequently used in species distribution models were sourced from the Worldclim website (Bioclim version 2.1) (www.worldclim.org, accessed on 1 December 2021). This dataset includes 19 variables that were derived from long-term temperature and rainfall data. For correlative models, a pairwise Pearson's correlation test was used to eliminate variables that were highly correlated. This led to selection of 10 variables with a correlation less than 0.75 to be used in the model training process. The variables used in MaxEnt and MMF included annual precipitation, annual mean temperature, minimum temperature of the coldest month, mean diurnal range, maximum temperature of the warmest month, temperature seasonality, precipitation of the coldest quarter, precipitation of the wettest quarter, and precipitation of the wettest month. In addition, variables were screened for ecological relevance [55]. To ensure uniformity among models, gridded climate data with 10′ resolution were used for all modelling approaches.

Species Distribution Modelling
CLIMEX as a semimechanistic model has been frequently used in plant disease distribution modelling efforts [56][57][58] and is considered a mechanistic model that enables the user to estimate the potential geographical distribution and seasonal abundance of a species in relation to climate [59][60][61]. Unlike statistical models, CLIMEX matches the patterns of climate and species' distribution by calibration rather than by a statistical fitting process. The parameter values used in [40] were initially used as starting points to calibrate the model. After adjusting the parameter values iteratively, the model was run repeatedly to achieve the highest ecoclimatic index (EI) values close to the known distribution of Psa and lowest EI values outside the range of Psa reported distribution. These criteria were used to obtain the closest match of suitable habitats projected by CLIMEX and the reported distribution patterns of Psa. The final parameter values used to develop the model are shown in Table 1. Finally, the model was assessed visually based on the match of known presence points and suitable habitat projected by the model.

Bioclimatic Variables
Climate data files (Met manager version 1.2) from CliMond website were downloaded for direct use in CLIMEX software (www.climond.org, accessed on 1 December 2021) [53]. Bioclimatic variables [54] that have been frequently used in species distribution models were sourced from the Worldclim website (Bioclim version 2.1) (www.worldclim.org, accessed on 1 December 2021). This dataset includes 19 variables that were derived from long-term temperature and rainfall data. For correlative models, a pairwise Pearson's correlation test was used to eliminate variables that were highly correlated. This led to selection of 10 variables with a correlation less than 0.75 to be used in the model training process. The variables used in MaxEnt and MMF included annual precipitation, annual mean temperature, minimum temperature of the coldest month, mean diurnal range, maximum temperature of the warmest month, temperature seasonality, precipitation of the coldest quarter, precipitation of the wettest quarter, and precipitation of the wettest month. In addition, variables were screened for ecological relevance [55]. To ensure uniformity among models, gridded climate data with 10 resolution were used for all modelling approaches.

Species Distribution Modelling
CLIMEX as a semimechanistic model has been frequently used in plant disease distribution modelling efforts [56][57][58] and is considered a mechanistic model that enables the user to estimate the potential geographical distribution and seasonal abundance of a species in relation to climate [59][60][61]. Unlike statistical models, CLIMEX matches the patterns of climate and species' distribution by calibration rather than by a statistical fitting process. The parameter values used in [40] were initially used as starting points to calibrate the model. After adjusting the parameter values iteratively, the model was run repeatedly to achieve the highest ecoclimatic index (EI) values close to the known distribution of Psa and lowest EI values outside the range of Psa reported distribution. These criteria were used to obtain the closest match of suitable habitats projected by CLIMEX and the reported distribution patterns of Psa. The final parameter values used to develop the model are shown in Table 1. Finally, the model was assessed visually based on the match of known presence points and suitable habitat projected by the model.

MaxEnt
MaxEnt is a presence-only or presence-background model which has been claimed to outperform most of existing modelling approaches [62]. A good explanation of MaxEnt, both with respect to its conceptual basis, its application and advantages can be found in Elith et al. [63]. In this study, suggestions by Merow et al. [64] and Syfert et al. [65] regarding feature selection and sampling bias were taken into account. Because enough presence data were available for Psa to allow model complexity, all available feature types i.e., linear, quadratic, product, threshold, and hinge were used as per the default settings of the MaxEnt model. In addition, to measure the variable importance a Jackknife test was performed. The response curves of the predictor variables were produced using available options in MaxEnt to investigate the dependency of predicted suitability on the selected variables. To validate the model output, MaxEnt was set to randomly use 75 percent of the presence data for training and 25 percent for testing the output.

Multimodel Framework
The multimodel framework was used to model species presence/absence data by seven different species distribution models: logistic regression (LOG), naive Bayes (NB), classification and regression trees (CART), conditional trees (CTREE), K-nearest neighborhood (KNN), support vector machines (SVM), and artificial neural networks (NNET). These models have been well-used in ecological studies and detailed descriptions about them can be found in the general literature on species distribution models [54]. All seven models were trained (or fitted) and tested using selected variables. The multimodel framework uses one-class support vector machines (OCSVMs) to select appropriate absence points out of large datasets. Instead of selecting a single best performing OCSVM model that can result in over-fitting, a set (an ensemble) of 100 models which had the lowest prediction errors were chosen to profile the background data. Pseudoabsence points were then selected from sites that were assigned 0 (probability of occurrence) on the profiled background data. As there were still many possible absence locations after this analysis, the generated pseudoabsence points were reduced by clustering locations that had similar environmental variables by a defined number of clusters to balance the number of presence locations [66].
Following the selection of the pseudoabsence data by the model, ten selected bioclimatic variables mentioned earlier were used in all seven models. Model validation was carried out to test the ability of the models to predict new data. This was done by bootstrap- ping (resampling data) and cross-validation (10-fold) to generate independent validation data. Model predictions were plotted globally. The predicted maps were exported in ASCII format for visualization ArcGIS Pro.

Model Consensus (Agreement)
The locations where all three models agreed on suitability of Psa establishment were identified by converting the three predicted maps into binary maps. To create the binary maps, the occurrence threshold of 0.5, the threshold of 10 percentile presence in the training data, and an Ecoclimatic Index (EI) threshold of >1 were used for, MMF, MaxEnt and CLIMEX, respectively [50,66,67]. These binary maps were then overlaid using "equal to Frequency" tool in ArcgGIs Pro to visualize the hotspots.

Model Performance
According to the current distribution of Psa, the CLIMEX model represented a good fit to the presence data. Table 1 shows the final parameter values used in the current study and those used in an earlier study [40]. The parameter values in [40] and [43] were used as starting points and the current values were obtained through iterative calibration until the map projection fitted the current distribution of Psa. The resultant parameter values were also checked for biological feasibility.
For MaxEnt, the area under the curve (AUC) was used to measure model performance. The AUC of both training and test data for the Psa model was 0.98, which indicates MaxEnt was able to discriminate the data very well ( Figure S1).
Multimodel, as an ensemble framework, produces a table of ten different performance criteria, which are separately ranked, and an overall score is calculated from these rankings to find the best-performing model. These performance criteria include accuracy, precision, sensitivity (recall), F.Score, Kappa index, specificity, true statistic skill (TSS) uncertainty, CV. error and AUC. A detailed description of these criteria can be found in [66]. Table S2 shows the variability of the model's performance according to ten performance criteria calculated by cross-validation and bootstrapping. Although variability among modelling methods is expected, there was a low overall variation in the performance among all nine models in MMF. For Psa, the model which had the highest rank in both validation methods was SVM, the results of which are interpreted and used for further analysis here.

Predictor (Environmental) Variables
According to CLIMEX model, the lower temperature threshold (DV0) was 5 • C, which is lower than the 8-10 • C that has been reported in similar studies [37,40]. In addition, the upper temperature threshold obtained was 28 • C which was higher than reported values in previous studies (25 • C) and those achieved by MaxEnt in this study [6]. These differences are expected as reported temperature thresholds are based on experiments conducted in controlled conditions [37].
By removing highly correlated variables through the Pearson correlation test (>75%), ten variables mentioned in the Methods section were used in MaxEnt ad MMF. The variable contribution analysis by MaxEnt showed that annual precipitation and annual mean temperature contributed the most to the model prediction followed by minimum temperature of the coldest month (Table 2). The Jackknife test implied that certain variables such as mean annual temperature (An_Mean_Temp), mean temperature of the coldest month (Mean_Temp_Coldest_Month), minimum temperature of the coldest month and annual precipitation were the most important variables and provided a reasonably good fit to the training data ( Figure S2). Furthermore, Figure S2 showed that no variable contained significantly useful information that was not previously included within the other variables. Compared to other variables, the gain in training data experienced a slight decrease when mean diurnal range was removed from the process of training data.
The response curves of some important variables such as mean annual temperature and annual precipitation showed the dependence of predicted suitability on both selected variables. For example, Figure 3A indicates a mean annual temperature of around 17 • C resulted in the highest probability of Psa presence which was in accordance with the CLIMEX results for the lower optimum temperature range (DV1-DV2) and that reported in the literature [37]. There was also a gradual drop in probability of Psa presence at temperatures above 17 • C and around 24 • C the probability of Psa presence became zero, which is in accordance with the report that Psa growth decreases at temperatures higher than 25 • C [37]. The response curve of annual precipitation (cbio12) indicated that excessive rain (>1200 mm/year) can decrease the probability of Psa presence ( Figure 3B).  The Jackknife test implied that certain variables such as mean annual temperature (An_Mean_Temp), mean temperature of the coldest month (Mean_Temp_Cold-est_Month), minimum temperature of the coldest month and annual precipitation were the most important variables and provided a reasonably good fit to the training data (Figure S2). Furthermore, Figure S2 showed that no variable contained significantly useful information that was not previously included within the other variables. Compared to other variables, the gain in training data experienced a slight decrease when mean diurnal range was removed from the process of training data.
The response curves of some important variables such as mean annual temperature and annual precipitation showed the dependence of predicted suitability on both selected variables. For example, Figure 3A indicates a mean annual temperature of around 17°C resulted in the highest probability of Psa presence which was in accordance with the CLIMEX results for the lower optimum temperature range (DV1-DV2) and that reported in the literature [37]. There was also a gradual drop in probability of Psa presence at temperatures above 17 °C and around 24 °C the probability of Psa presence became zero, which is in accordance with the report that Psa growth decreases at temperatures higher than 25 °C [37]. The response curve of annual precipitation (cbio12) indicated that excessive rain (>1200 mm/year) can decrease the probability of Psa presence ( Figure 3B).

Potential Distribton
Large areas in China, United States, Europe and South America were predicted as suitable by all three models (Figure 4a-c). Nevertheless, there were noticeable discrepancies among the models' projections in some areas such as southeastern parts of the USA,

Potential Distribton
Large areas in China, United States, Europe and South America were predicted as suitable by all three models (Figure 4a-c). Nevertheless, there were noticeable discrepancies Climate 2022, 10, 14 9 of 17 among the models' projections in some areas such as southeastern parts of the USA, Eastern Europe, Central Asia, and South America. In North America, while both MaxEnt and CLIMEX projected large areas in southeastern states of the United States as highly suitable, MMF did not project these areas as suitable and instead northern parts of United States were predicted as highly suitable. California was projected as highly suitable by all three models, in particular in the kiwifruit orchards centered in the Sacramento Valley. In South America, the most extensive areas that were projected as highly suitable for Psa were located around Buenos Aires and Mar del Planta in Argentina. In addition, MMF prediction for Argentina extended further south (around Comodoro Rivadavia) while these areas were not predicted as highly suitable by MaxEnt and CLIMEX. In addition to Argentina, the whole of Uruguay and southern Brazil were also projected as highly suitable. In Chile (where kiwifruit is grown), from Puerto Montt to Santiago were predicted highly suitable by all models.
Central Asia was not predicted as suitable by MaxEnt and CLIMEX and Eastern Europe was predicted marginally suitable, while MMF projected Central Asia and larger areas in Europe (extending toward Eastern Europe) as highly suitable. All models projected Japan and most parts of South Korea as suitable which agrees with current reports of the disease from these regions. All three models projected that northern parts of Iran where kiwifruit is grown commercially were highly suitable. Kiwifruit is also grown in a few orchards in some countries such as Laos, Philippines, Cambodia, Vietnam, India and Bangladesh. Among those countries only Laos was predicted suitable by all three models.
In Africa, some coastal areas in Morocco and Algeria were projected as highly suitable by all three models, but kiwifruit is not currently grown in these areas. Kiwifruit have been cultivated in South Africa only recently and all models project areas around Cape Town and Durban as highly suitable for long-term establishment of Psa.
In Australia, all three models predicted the southeast coastal areas such as Victoria are highly suitable, although so far only Psa-LV (biovar 4) has been reported from these areas [68]. Perth in Western Australia, where kiwifruit is grown on a small scale, was also projected as highly suitable for Psa by all models.
In New Zealand, the North Island was predicted highly suitable for Psa establishment by all three models with all areas having EI values greater than 25 according to the CLIMEX model. While SVM and CLIMEX projected large areas in both North Island and South Island of New Zealand as highly suitable, MaxEnt projection for the South Island was limited to northern parts around Nelson. The prediction agrees with the current distribution of Psa in New Zealand. The weekly growth index chart produced by CLIMEX ( Figure 5) shows that on most days of the year, Psa can establish in the Bay of Plenty, an important kiwifruit growing area, although the most suitable periods seem to be March-April and October-December. This higher growth index in winter was in accordance with peaks of Psa reported in New Zealand (www.kvh.org.nz, accessed 1 December 2021). In summary, all kiwifruit growing areas in New Zealand lie within the suitable category of CLIMEX confirming that New Zealand is highly suitable for long-term establishment of Psa.

Model Consensus
The consensus model highlighted the areas (hotspots) that all models agreed on suitability of Psa establishment ( Figure 6). From 201 presence points, 199 points were located in areas where all three models agreed on suitability. The remaining two points (in Sichuan (China) and Gijon (Spain)) were located in areas predicted suitable by two models. Areas that were predicted as suitable by all three models included New Zealand, Australia, large areas of southern and southeast China, Japan, South Korea, large areas of Europe, Turkey, coastal areas of South Africa, California, Chile, Southern Brazil, Uruguay, and Buenos Aires in Argentina. Climate 2021, 9, x FOR PEER REVIEW 10 of 17 In New Zealand, the North Island was predicted highly suitable for Psa establishment by all three models with all areas having EI values greater than 25 according to the CLIMEX model. While SVM and CLIMEX projected large areas in both North Island and South Island of New Zealand as highly suitable, MaxEnt projection for the South Island important kiwifruit growing area, although the most suitable periods seem to be March-April and October-December. This higher growth index in winter was in accordance with peaks of Psa reported in New Zealand (www.kvh.org.nz, accessed 01 December 2021). In summary, all kiwifruit growing areas in New Zealand lie within the suitable category of CLIMEX confirming that New Zealand is highly suitable for long-term establishment of Psa.

Model Consensus
The consensus model highlighted the areas (hotspots) that all models agreed on suitability of Psa establishment ( Figure 6). From 201 presence points, 199 points were located in areas where all three models agreed on suitability. The remaining two points (in Sichuan (China) and Gijon (Spain)) were located in areas predicted suitable by two models. Areas that were predicted as suitable by all three models included New Zealand, Australia, large areas of southern and southeast China, Japan, South Korea, large areas of Europe, Turkey, coastal areas of South Africa, California, Chile, Southern Brazil, Uruguay, and Buenos Aires in Argentina.

Model Consensus
The consensus model highlighted the areas (hotspots) that all models agreed on suitability of Psa establishment ( Figure 6). From 201 presence points, 199 points were located in areas where all three models agreed on suitability. The remaining two points (in Sichuan (China) and Gijon (Spain)) were located in areas predicted suitable by two models. Areas that were predicted as suitable by all three models included New Zealand, Australia, large areas of southern and southeast China, Japan, South Korea, large areas of Europe, Turkey, coastal areas of South Africa, California, Chile, Southern Brazil, Uruguay, and Buenos Aires in Argentina.

Discussion
Only four studies have been developed to model Psa distribution using SDMs [40,[43][44][45]. All these models focus on the potential distribution of Psa at a local scale and mostly rely on single models or the same type of models such as correlative models [43][44][45]. While some studies benefit from an ensemble approach [45], it is recommended that where possible, contributions to an ensemble model should be from different types of models and they should add information for building a better model [69]. The current study is the first to address the potential global distribution of an important plant pathogen using three different modelling approaches, a semi-mechanistic model (CLIMEX), a presencebackground correlative model (MaxEnt) and an ensemble of presence-pseudo absence models (MMF). To deal with uncertainty and benefit from different models and algorithm predictions, the models' outputs were combined into a consensus model which focuses on prediction agreements.
Each individual model showed a good fit to the current global distribution of Psa. While there was a high level of agreement among the three models on highly suitable areas in China, Australia, New Zealand, Europe, the USA and some parts of South America, compared to CLIMEX and MaxEnt, SVM projected larger areas in China, Central Asia and Europe as highly suitable.
The parameter values in the CLIMEX model were slightly different from the model developed in a preliminary study by [43]. Current CLIMEX model parameter values are more in agreement with the ranges mentioned in earlier studies [6,38]. The CLIMEX model indicated that the minimum temperature threshold (DV0) is 5 • C, which means the pathogen can stay active at this temperature. This is close to the finding by [44] where temperatures between 3.6-17.1 • C in October were found to be conducive to Psa establishment in China. However, in the CLIMEX model presented in the EPPO report [40] and [37], a minimum temperature threshold of 10 • C (DV0) was reported, which differs from the current findings and the information available on temperature needs of the pathogen. In the CLIMEX model presented in the EPPO report, the optimum temperature range is 20 • C to 25 • C, which is higher than the current findings (12 • C to 20 • C) and those of previous studies regarding the optimum temperature range for Psa growth [37,44]. The response curves produced by the MaxEnt model (Figure 3), also indicated that the probability of pathogen establishment is the highest between 10 and 20 • C with optimum being around 13 • C. This is in close agreement with the current CLIMEX findings on suitable temperature range. Both the CLIMEX and MaxEnt models in the current study indicated that above 20 • C, the probability of Psa establishment starts to decrease, which is similar to the findings in [44] that showed the probability of Psa establishment rapidly decreases above 21.2 • C.
The soil moisture parameters in the current study are slightly higher than [40], which indicates that based on the current CLIMEX model, the pathogen is favored by higher soil moisture content compared with the CLIMEX model presented in the EPPO report. In CLIMEX, soil moisture is used as a proxy for rainfall and evapotranspiration, meaning that a higher soil moisture index implies increased rainfall and humidity affecting species growth [70]. Thus a higher soil moisture parameter value reflects more optimal rainfall, as spread of the foliar and stem pathogen Psa is not directly affected by soil moisture, and infection through roots is practically absent [71]. While the ensemble model developed in [45] showed that annual precipitation plays a minor role in Psa spread, similar to the model in [44], our study shows that precipitation plays an important role in Psa establishment. The sensitivity of Psa to rainfall has also been highlighted in other studies [48,72] and is in agreement with the high soil moisture (SM) values indicated in the current CLIMEX model ( Table 1). The response of Psa to annual precipitation shows that while the pathogen is favored by moist conditions, excessive rain (above 1200 mm/year) will decrease the probability of pathogen establishment, probably because of the washout of inoculum [71].
The reasons for these differences in model outcomes might be due to coarser climate data and the fact that the previous CLMEX study was conducted when less information on the distribution of Psa was available. For example, Psa occurrences from Spain had not been used in the modelling effort mentioned in the EPPO report [40]. On the contrary, occurrence data for Perth in Western Australia were used in the EPPO study but these occurrence(s) proved to be erroneous later [72]. The number of Psa records used in model development can significantly affect results for species that are still spreading. In the current study we used 201 presence points while [43] used 82 presence points, explaining part of the differences found.
While the EI value estimated by the EPPO model [40] for New Zealand ranges from 15 to 35, indicating suitable to highly suitable areas for Psa establishment, the EI values estimated by the current model range from 50 to 90 suggesting an even higher suitability. Furthermore, our model projected the north of Iran (an area in which kiwifruit is grown commercially) as highly suitable but it was not projected as suitable or even marginally suitable by the CLIMEX model in the EPPO study [40]. To date there has been no report of Psa presence in this region. It should be noted, however, that one study mentions the presence of kiwifruit bacterial canker in Iran but refers to the causal agent as P. syringae pv. syringae [73]. In East Asia, where China is the major kiwifruit producer, all models perfectly fitted the occurrence data. At a national level, Yunnan and Guangxi province (where kiwifruit orchards cover 500 to 1000 ha respectively) were projected as suitable by all models. Although the number of predictor variables, the number of presence data and the resolution of the data used in [44] and 45 were different from the current study, overall, the model projection patterns for China are similar to those reported by [44] and 45 except that MMF projected larger areas in northern parts of China (Figure 4a-c). In Europe, all three models agreed on the suitability of some countries such as Italy, Spain, Portugal, France, Turkey, and Greece. The Trabazon Rize province in Turkey, where there are few reports of Psa presence, is predicted as suitable by all three models. The consensus model highlighted the areas where Psa establishment has been predicted by all three models. The areas of interest are California in the USA, northern parts of Iran, unaffected areas in Greece, Belgium, Denmark, and South Africa ( Figure 6).
While prediction models, especially SDMs, can be valuable tools in assessing the risk of damaging invasive species, the limitation, caveats, and inherent uncertainties should be considered while interpreting their results [74][75][76][77]. In general, for species that may not have changed their environmental range, the presence profiles are expected to be confined to a small area in multidimensional space. The widespread predicted presence of Psa may be an indication that this species is still spreading, and that we may underestimate the full range of its potential distribution [43,50]. This result questions the validity of applying presence-only models such as MaxEnt that rely only on current knowledge of species that are still spreading. Moreover, it is important to note that the presence of a species may not be recorded for several reasons such as, (1) difficulty in detecting the species (an important consideration for plant pathogens due to their microscopic size), (2) evolutionary change of the pathogen, becoming more aggressive, and (3) the sudden appearance of new pathways aiding species spread.
Despite the limitation of individual models, the consensus model based on model agreement presented in this study, composed of three different types of models as recommended by [78], minimized the uncertainties and resulted in predictions that reflected the observed occurrences and the predictions from previous models [40,43]. The consensus model had the highest sensitivity and lowest false negative rate, with 199 out of 201 presence points lying within areas predicted as suitable. In situations involving important biosecurity threats such as Psa, it is important to minimize false negatives and aim for models with lower emission errors [79,80]. While over-prediction can also be an issue due to implementing unnecessary increased biosecurity readiness and response programs and strategies, prior knowledge of species biology can be used to prioritize areas which are most at risk so that when limited resources are available, they could be efficiently allocated in the right place at the right times.
As there are no known cures for Psa, the control of the disease should primarily rely on preventive measures. The results from predictive modelling suggest that optimal moisture supply, avoiding overhead irrigation, may contribute significantly to reduce Psa spread. In addition, management measures can include preventive application of chemicals such as copper-based compounds and antibiotics (streptomycin) to contain the pathogen spread. However, effective management would not be achieved without proper management procedures in the field including orchard hygiene, controlled access to orchards, monitoring grafting and use of nutrients for a healthier orchard. In addition, due to environmental concerns about application of chemicals, alternative methods such as biocontrol strategies using specific bacteriophages to eliminate Psa are currently under study and some already showed promising results [81].
Apart from some discrepancies among the models, all three models highlighted suitable areas in which Psa may establish. The high level of agreement between the SVM model (from the multimodel framework) and CLIMEX and MaxEnt increased confidence with respect to model projections for novel areas. These results are particularly valuable indicators for countries where Psa is currently reported from limited localities and for the USA, Iran, South Africa, Belgium, and Denmark, where Psa has not been reported.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cli10020014/s1, Figure S1: Receiver operating characteristic (ROC) curve produced by MaxEnt for Psa; Table S1: The geographic coordinates of locations of Psa. Some presence data such as Australian points are excluded as they are now known to be a different pathovar; Table S2