The geodatabase described in Tubby et al. (2023), generated following the approaches used for examinations of Dothistroma species and Fusarium circinatum Nirenberg and O’Donnell range and impact (Drenkhan et al., 2016, 2020), was used to describe the current geographic distribution of L. acicola. A global consortium of specialist researchers contributed to data collection with the objective of collating records and locations of L. acicola worldwide. Certain data entry fields were mandatory, including pathogen presence/absence, identity of data holder, host species, date of record, and forest type (urban, natural, plantation). Voluntary fields allowed entry of more detailed data on disease severity, presence of other pathogens, local climate, soil type and management practices, among other variables (see Tubby et al., 2023; Supplementary Table 1 for details).

The geodatabase is projected into the WGS 1984 coordinate system. All data are publicly available on a web map application hosted by Mendel University in Brno, Czechia.¹ The geodatabase is a live system with periodic additions. In the current study the database version from the 16th of December 2022 was used for model development. The geodatabase contained 3,081 records from Asia, Europe, North America, South America, Oceania, and Africa. Only data from the wider environment were used (n = 2,970), nursery records were omitted. The database contained 896 confirmed reports of L. acicola and 2,074 negative records where L. acicola was marked as absent (Tubby et al., 2023).

The majority of the presence/absence data were confirmed by molecular methods (e.g., conventional PCR, qPCR, rtPCR), i.e., 79.8% of the absence data and 48.9% of the presence data (Tubby et al., 2023). In addition, 37.9% of the presence data and 6.0% of the absence data were identified by morphological examination with a compound microscope. A smaller proportion of data were identified by visual inspection only, 14.2 and 13.2% of the absence and presence data, respectively. Most presence/absence data came from official systematic surveys and are therefore considered highly reliable (Tubby et al., 2023).

The gymnosperm database was used to establish a global distribution of pines (Earle, 2021; accessed 24. July 2022). The database contained 6,109 records from the BRAHMS database at Conifers of the world (Farjon, 2021; accessed 11. November 2021), and four GBIF Occurrence download (GBIF, 2022; accessed 12. November 2021). The database considered 119 Pinus species. We assumed that all Pinus species were equally susceptible to L. acicola infection and that their distribution and susceptibility would not change under climate change.

The worldwide distribution of Pinus spp. was generalized to 1 × 1-degree cells taking into account the native range of taxa and extralimital occurrences, e.g., agroforestry, horticulture, or naturalized. In total 3,769 model cells with pines occurred (Figure 2).

Therefore, the spatial resolution of the model was 1 × 1-degree cells. The geographical distribution of L. acicola was then matched to the spatial resolution of the model. L. acicola was present in 222 and absent in 182 model cells (Figure 2). If a model cell contained positive and negative records of L. acicola, the model cell was marked as positive for L. acicola. The dataset contained nearly balanced presence/absence data. However, sampling bias was not investigated.

Initially, data for 44 independent variables were collected and converted into the model spatial resolution (Supplementary Table 1). Elevation and 19 bioclimatic variables were acquired from WorldClim version 2.1 climate data for 1970–2000 (Supplementary Table 1; Fick and Hijmans, 2017, 2020).

Four variables from the global high resolution soil water balance dataset were included (Trabucco and Zomer, 2010, 2019; average over 1950–2000 period): potential evapotranspiration (PET), actual evapotranspiration (AET), monthly fraction of soil water content available for evapotranspiration processes (as percentage of maximum soil water content; a measure of soil stress), and the Priestley-Taylor alpha coefficient (generalized as the ratio of annual AET over annual PET; a description of overall aridity stress on vegetation).

Furthermore, the global aridity index was obtained from the global aridity index and PET database (Zomer et al., 2022) and expressed as a ratio of mean annual precipitation over mean annual reference evapotranspiration for the averaged 1971–2000 time period.

Four topsoil properties were included in the database: organic carbon (% of weight), pH, reference bulk density (kg/dm³), and base saturation (%). Soil data was acquired from the harmonized world soil database (FAO et al., 2009).

Additionally, five variables were included from monthly global climate historical data computed in the framework of the sixth phase of the Coupled Model Intercomparison Project (CMIP6): near-surface relative humidity (%), near-surface wind speed (m/s), moisture in the upper portion of the soil column (MRSOS, vertical sum per unit area from the surface down to the bottom of the soil model of water in all phases contained in soil; kg/m²), surface downwelling longwave radiation (RLDS, radiative longwave flux of energy downward at the surface; W/m²), surface downwelling shortwave radiation (RSDS, radiative shortwave flux of energy downward at the surface; W/m²). All five variables were included in three variants (altogether 15 variables), i.e., average values in the period 1971–2000 for driest quarter, warmest quarter, and wettest quarter where quarter is a period of 3 months.

Bioclimate and soil variables were available with a spatial resolution of 30 s (approximately 1 km; Supplementary Table 2). However, relative humidity, wind speed, MRSOS, RLDS and RSDS were available at a much coarser scale, i.e., 0.5–2.8° (approximately 55.6–311.2 km). We decided to use a spatial resolution of 1-degree to match the spatial resolution of the Pinus dataset. Finally, the data with various spatial resolutions were resampled to a 1-degree model grid using bilinear interpolation.

Autocorrelation between variables was tested using the FindCorr function from DescTools package version 0.99.45 (Signorell et al., 2022) in R statistical software (R Core Team, 2022) using 0.90 pair-wise absolute correlation cut-off threshold. The FindCorr function automated a search for autocorrelated independent variables. FindCorr is based on the correlation matrix calculated using the “cor” function where Pearson’s correlation coefficient was chosen as the measure of linear correlation between two variables. Before using the “cor” function, it was checked that the assumptions of the Pearson’s correlation coefficient were met: (1) the variables are numerical and continuous, (2) there is a linear relationship between the variables, (3) the variables are normally distributed, and (4) the data are from a representative sample. All variables were numeric, continuous, had a visible linear relationship, and the data sample was representative. Variables that were not normally distributed were log transformed to achieve approximate normality (elevation, BIO3, all precipitation variables, aridity index, ETP, topsoil organic carbon, topsoil base saturation). Altogether ten highly autocorrelated variables were removed from further model development (Supplementary Table 1 and Supplementary Figure 1).

Further variable selection was done using random forests with the VSURF package version 1.1.0 (Genuer et al., 2010, 2019). VSURF incorporates three steps for variable selection based on random forests for supervised classification and regression problems. The first “thresholding step” concentrates on eliminating irrelevant variables from the dataset. The second “interpretation step” selects all variables related to the response for interpretation purposes. The third “prediction step” refines the selection by eliminating redundancy in the set of variables selected by the second step, for prediction purposes. A detailed description of the individual steps can be found in the VSURF package documentation (Genuer et al., 2010, 2019). One hundred random forests were grown in each step and each random forest had 300 classification trees.

Four modeling approaches were used: (1) generalized linear model (GLM) with binomial distribution with logit link function, (2) individual classification trees (random tree, J48, LMT), (3) bagging (with three different base learners: random tree, J48, LMT), and (4) random forest (RF).

Generalized linear models (GLMs) are frequently used, regression-based SDMs. GLMs can handle non-normal error distributions, additive terms and non-linear fitted functions (McCulloch et al., 2005). GLMs generalize “linear regression by allowing the linear model to be related to the response variable via a link function and by allowing the magnitude of the variance of each measurement to be a function of its predicted value” (Zhao, 2013). The GLM consists of three elements: (1) a probability distribution, (2) a linear predictor, and (3) a link function. A binomial probability distribution with logit link function is most commonly chosen for presence/absence data.

Classification trees are decision trees which are among the most popular tools for machine learning and data mining (Blockeel, 1998). “A decision tree is a hierarchical structure where each internal node uses a test on an attribute, each branch corresponds to an outcome of the test, and each leaf gives a prediction for the value of the dependent class variable” (Džeroski, 2001). The paths from root to leaf represent classification rules. A number of classification trees have been developed since the 1990s including J48, Random tree and LMT (logistic model tree). In our study we used three methods: J48, Random tree and LMT. The J48 method is used for generating a pruned or unpruned C4.5 decision tree (Quinlan, 1993). Random tree constructs a classification tree that considers a number of randomly chosen attributes at each node and performs no pruning (Breiman, 2001). LMT builds logistic model trees, which are classification trees with logistic regression functions at the leaves (Landwehr et al., 2005; Sumner et al., 2005).

”Bagging is an ensemble method that constructs the different classifiers by making bootstrap replicates of the training set and using each of these replicates to construct one classifier. Each bootstrap sample is obtained by randomly sampling training instances, with replacement, from the original training set, until an equal number of instances is obtained” (Breiman, 1996, 2001; Kocev et al., 2007). Breiman (1996) has shown that bagging can give substantial gains in predictive performance when applied to an unstable learner such as classification and regression tree learners. The classification tree learners J48, Random tree and LMT are also available in the bagging method.

A random forest is now widely used in SDMs (e.g., Mi et al., 2017; Watt et al., 2021; Koldasbayeva et al., 2022). “A random forest is an ensemble of trees, where diversity among the predictors is obtained by bagging, and additionally by changing the feature set during learning. More precisely, at each node in the decision trees, a random subset of the input attributes is taken, and the best feature is selected from this subset” (Breiman, 2001). “Random forests are a combination of tree predictors such that each tree depends on the values of a random vector sampled independently and with the same distribution for all trees in the forest” (Breiman, 2001; Kocev et al., 2007). It is an ensemble method which is better than a single classification tree because it reduces over-fitting by averaging the result (Jin et al., 2014; Carvajal et al., 2018; Nahar and Ara, 2018; Naseem et al., 2021). The final predictions in random forest are calculated from an ensemble of predictions from individual trees used as a median prediction. Different studies have shown that random forest outperforms other methods such as extreme gradient boosting, multiple regression (Watt et al., 2021), geographically weighted regression, artificial neural network, and support vector machine for regression (Chen et al., 2019).

Altogether, eight models were developed (Table 1). The best settings for each method were explored by changing its parameters repeatedly. The settings that gave the highest accuracy were kept as optimal. The models were validated using 10-fold cross-validation. The three best models were used as an ensemble for the predictions where the most frequent prediction per model cell was considered.

The models can be further improved if spatial autocorrelation is taken into account. Spatial autocorrelation is the correlation between spatially adjacent values of a single variable. In our study, we did not consider spatial autocorrelation. Therefore, the bias due to spatial autocorrelation was unknown (Diniz-Filho et al., 2003).

Attribute importance was calculated based on Gini impurity. For a decision tree, the impurity decrease from each feature can be averaged, and the features are ranked according to this measure. Gini impurity gives the probability of misclassifying an observation. Therefore, the lower the Gini the lower the likelihood of misclassification, i.e., the lower Gini the higher variable importance. The Gini impurity has a minimum (highest level of purity) of 0 and a maximum value of 0.5 indicating a random assignment of classes. Gini impurity for the ith node is given by Breiman et al. (1998):

where p_i,k is the proportion of samples that belong to class k in the ith node of the tree.

Model performance was evaluated with mean absolute error (MAE), root mean squared error (RMSE), true positive rate [TPR or sensitivity = TP/(TP + FN)], false positive rate [FPR or fall-out = FP/(FP + TN)], precision = TP/(TP + FP), ROC area (Receiver Operating Characteristic area), and Kappa = (observed accuracy—expected accuracy)/(1—expected accuracy). ROC area is an area under the curve (AUC), where the TPR is plotted on the Y-axis, and the FPR on the X-axis (Witten and Frank, 2005).

The potential distribution of L. acicola was computed for five global climate models (GCM, Table 2) and three combined pathways of Shared Socioeconomic Pathway (SSP) and Representative Concentration Pathway (RCP): SSP1-RCP2.6, SSP2-RCP4.5, and SSP5-RCP8.5 (van Vuuren et al., 2011; Riahi et al., 2017). In the remainder of the text, we refer to these simply as RCP scenarios and they are used to illustrate different pathways for future climate forcing. The predictions were made for eleven 30-year periods between 1971 and 2100. The results of the five GCMs were averaged on combined SSP-RCP (median) per 30-year period. Data on the presence or absence of L. acicola were for the period 1922–2021, but not all climate data were available for this period. Most climate data were available for the period after 1970. Therefore, a standard meteorological period of 30 years (1971–2000) was chosen as the reference period.

The model’s results were clipped by the global distribution of pines (Earle, 2021) which resulted in maps of potential distribution of L. acicola. Potential distribution in the reference period was compared with the known distribution of L. acicola (Tubby et al., 2023) to get the realized distribution (Figure 1).

Uncertainty of the GCMs predictions were classified as (1) certain when all five GCMs predicted the same result, (2) low when one GCM had a different result, and (3) medium when two GCMs had different results.

The database was prepared and managed in Microsoft SQL Server 2016 (version 13.0.6419.1). Variable selection was done using R statistical software (R Core Team, 2022). The models using machine learning methods were developed with Weka 3.8.6 (Hall et al., 2009; Frank et al., 2016), GLM was developed with R. Maps were drawn with ESRI ArcMap 10.6.1. Charts were drawn with R and Microsoft Excel version 2304.

Eight of 44 independent variables were selected for the model development (Table 3).

There were 3,769 model cells covering the native and extralimital range of pines (Table 4). Europe and North America shared almost the same number of model cells containing pines, i.e., 33.1 and 32.5%, respectively. Asia followed with 25.6%. All other continents had a non-native distribution of pines covering under 5% of model cells.

The worst performing model was GLM by almost all evaluation indicators. Individual classification trees performed better than GLM. Among individual classification trees, the best was LMT and the worst, random tree. The bagging approach performed better than individual classification trees and GLM. Among bagging the best classifier was random tree and the worst J48. Random forest performed nearly the same, but slightly better than bagging with LMT. The best model was bagging with random tree, with a very high AUC (0.938) and high precision (0.878). The best three models were selected for the climate change predictions: bagging—random tree, bagging—LMT, and random forest (Table 5).

Variable importance, based on average Gini impurity, ranged from 0.33 to 0.40. The most important variables were RLDS, DIURNG, and PRECWARMQ. Elevation and TMPWETTQ had nearly the same Gini impurity. PRECSEASON was also identified as important for L. acicola distribution. The least important variables were MRSOS and RSDS. However, differences between Gini impurities of the variables were small (Figure 3).

Data on the occurrence of L. acicola were mainly from Europe and North America (Figure 4). BSNB was also detected in a small number of cases in Asia and South America, but not in Oceania and Africa. Most of the absence data were recorded in Europe. Absence data were available for all continents in the geodatabase.

The actual distribution of L. acicola in the reference period covered 5.9% of Pinus spp. area globally (Tubby et al., 2023), but the model ensemble indicated suitable environmental conditions for the pathogen across 58.2% of Pinus spp. cover in this period (Figure 4).

The most suitable conditions for L. acicola in 1971–2000 were in Asia, North America and Africa, where the potential distribution was computed to cover 75, 69, and 62% of continental Pinus spp. area, respectively, taking into account the median of five GCMs and the RCP4.5 climate change scenario. In contrast, only 35% of Pinus spp. area in Europe fitted the model conditions of L. acicola in 1971–2000 (Figures 4, 7).

The comparison between the actual distribution and the potential distribution of L. acicola, taking into account median predictions based on historical data of five GCMs in the period 1971–2000, showed that 81.6% of model cells were classified correctly, i.e., 40.9% cells predicted to be absent for L. acicola and 40.7% present. In some cases, the model predicted presence of L. acicola but actual distribution data showed an absence of the disease (4.0% of model cells, false positive predictions). This type of potential expansion of actual range was mostly detected in Europe (13 model cells) and in Africa (three model cells). In contrast, in some locations the model ensemble predicted absence of the disease, but actual distribution showed presence (14.4% of model cells). These false negative median predictions were again, mostly in Europe (50 model cells), with a few in North America (six model cells) and Asia (one model cell) (Figure 5).

Globally, actual distribution of L. acicola covered 10.1% of potential distribution in the reference period (Table 6). The highest realized distribution of L. acicola was determined in South America with disease known to be present in 22.7% of potentially suitable areas. In second place was Europe, with 21.1% of realized distribution. Asia had only 2.3% realized distribution. In Africa and Oceania, the disease was absent.

The different climate change scenarios (five GCMs, three RCPs) showed a positive trend in possible range expansion of L. acicola for the period 2001–2100. The average model predictions toward the end of the century (2071–2100) show the potential distribution of L. acicola rising to 62.2, 61.9, 60.3% of Pinus spp. area for RCP2.6, RCP4.5, RCP8.5, respectively. Predictions for different RCPs were almost consistent up to 2031, smaller differences are expected after 2031. However, on a global scale the relative change in median prediction from 2001 to 2100 was only 4.0, 3.7, 2.1% of Pinus spp. area for RCP2.6, RCP4.5, RCP8.5, respectively. In contrast, the 95% confidence interval was 35.7–82.3% of Pinus spp. area in the period 1971–2000 and 33.6–85.8% in the period 2071–2100 (Figure 6).

The linear relationship between diseased pine area and time was statistically significant (p < 0.001). Linear correlation for all three RCPs was positive and very high to medium, i.e., Pearson correlation coefficient measured 0.97, 0.84, 0.64 for RCP2.6, RCP4.5, RCP8.5, respectively. The average increase of potential global distribution of L. acicola measured 0.05, 0.04, and 0.01% per year for RCP2.6, RCP4.5, RCP8.5, respectively.

On all continents potential distribution of L. acicola will be higher at the end of 21st century than in the reference period 1971–2000. The highest relative change in the potential distribution of L. acicola, compared to the potential distribution in 2000, considering the median of five GCMs and the RCP4.5 climate change scenario until the end of 21st century, is expected in Oceania, Australia, and Africa with 175, 26.8, and 14.1% increases, respectively, compared to the potential distribution in 1971–2000. By 2071–2100 the model predicts an increase in the potential distribution of L. acicola of 12.2% in Europe, 10.5% in South America and 10.0% in Asia. In North America a negligible change of 3.1% is expected. Intriguingly, the potential distribution of L. acicola in Europe generally decreases between 2040 and 2070, although thereafter rapidly rises from 2071 to 2100 (Figure 7).

A spatial comparison between the potential distribution of L. acicola in the reference period 1971–2000 and the median prediction of five GCMs considering RCP4.5 in the period 2071–2100 showed that most of the model cells (87.7%) will not change status (i.e., presence/absence of L. acicola will not change). The model predicted that L. acicola could expand its potential distribution to an additional 8.0% of Pinus spp. area, to cover a total of 61.9% of Pinus spp. area. Most of this expansion was registered in Europe (3.2%), Asia (2.1%), and North America (1.8%). Expansion was detected in a northward direction in the Northern Hemisphere, and southward in the Southern Hemisphere. A reduction in the potential distribution of the disease was predicted in 4.3% of model cells, mostly in Europe (2.6%), North America (1.3%) and Asia (0.3%) (Figure 8).

Each model’s cell uncertainty of GCMs predictions was classified as low and medium as described in the methods. Low and medium uncertainty of GCMs predictions were identified in 46.7–51.8% of model cells for three RCPs in the period 1971–2100 (Figure 9). The lowest uncertainty was for RCP2.6 and the highest for RCP8.5, although differences between RCP uncertainties were small. Uncertainty increased from 1971 to 2100.

The highest uncertainty of the predictions was for Oceania and Europe and the lowest for Australia and South America (Figure 10). Australia, South America, Africa and Asia had similar uncertainties of the model ensemble predictions, e.g., uncertainty for 2071–2100 for those continents ranged from 22.6 to 34.6%. Uncertainty of the predictions for North America was medium. Uncertainty in model predictions increased over the period 1971–2100 for all continents except South America, where levels of uncertainties decreased.

There was high confidence (low uncertainty) in the model predictions for RCP4.5 in the period 2071–2100 for most of the south Europe, east Asia, most of North America except for the central states, Africa and Australia (Figure 11). Most of New Zealand had medium and low uncertainties. In Australia medium and low uncertainty was concentrated to Victoria, New South Wales and South Australia. In South Africa most of the model cells with low and medium uncertainties were found on the coastline between Cape Town and Gqeberha. In North America, uncertain model predictions were mostly found in Rocky Mountains and Great Plains. Uncertainty of GCMs predictions for central, east and west Europe were low to medium.