### Anopheles occurrence data

Occurrence data of the dominant *Anopheles* vectors in Europe come from the distribution maps of Sinka et al. [5]. These authors used the Malaria Atlas Project library, a literature research and expert opinion maps to produce maps of *Anopheles* occurrences. Occurrence data refer basically to the period 1985–2009. The presence-absence information of the six *Anopheles* vectors (*An. atroparvus*, *An*. *labranchiae*, *An*. *messeae*, *An*. *sacharovi*, *An*. *sergentii* and *An*. *superpictus*), which is provided as Shapefiles in the electronic supplementary material of Sinka et al. [5], was rasterized to a 0.25° resolution to match the grid resolution of the observation-based climate data.

### Climate data

#### Observational data

Mean, minimum and maximum temperature as well as precipitation data were taken from the daily 0.25° E-OBS dataset version 17 provided by the European Climate Assessment & Dataset (ECA&D [16]). A European-Mediterranean domain was selected, covering 75.375–25.375°N and 19.875°W49.875°E. Data in the time period 1950–2009 were selected. The time period 1985–2009, which is mainly used for subsequent analysis, was filtered for missing values for each month separately. A particular month was considered complete if there were less than three missing days per month, and the time series was considered complete if there were less than four missing months in the period 1985–2009. Grid boxes which did not meet these conditions were removed. Monthly temperature means (in Kelvin K) as well as monthly precipitation totals were calculated from the daily data.

#### Model data

RCM simulations carried out in the framework of EURO-CORDEX (European branch of the Coordinated Regional Climate Downscaling Experiment [17]) are employed. In the present study data with the grid resolution of about 50 km (0.44° on a rotated grid) are considered and conservatively regridded to 0.25° in order to match the E-OBS grid resolution. KNMI-RACMO22E driven by the GCM EC-EARTH as well as CLMcom-CCLM4-8-17 driven by MPI-ESM-LR were chosen as RCM simulations. The two RCMs were selected according to their good performance over Europe in the observational period [18]. Historical runs for the period 1950–2005 (KNMI-RACMO22E) and 1960–2005 (CLMcom-CCLM4-8-17) as well as scenario runs for the period 2006–2100 under RCP4.5 and RCP8.5 scenario assumptions [19] were available for subsequent analyses.

#### Land cover data

Land cover information was used in the form of 22 categories of land cover from the GlobCover project [20]. GlobCover3 2009 V.2.3 land cover map is derived by an automatic and regionally-tuned classification of a time series of global ERIS (Medium Resolution Imaging Spectrometer Instrument) Fine Resolution mosaics for the year 2009. The land cover classes, defined with the United Nations Land Cover Classification System, have a 300 × 300 m resolution. To match the E-OBS grid resolution, land cover information was calculated as percent coverage of each class in a particular 0.25° grid box.

### Correction of model data

Transfer functions were defined to match the RCM output of a variable *P*_{m} (temperature, precipitation) in the historical period with the statistical properties of a variable *P*_{o} from the E-OBS observations. The non-parametric empirical quantile method suggested in [21] and implemented in the *qmap* package in R was used. The transfer functions were subsequently used to correct the RCM output of the historical and the future periods. According to Gudmundsson et al. [21], the transformation is defined as:

$$ {P}_o={F}_o^{-1}\left({F}_m\left({P}_m\right)\right) $$

(1)

where *F*_{m} is the cumulative distribution function (CDF) of *P*_{m} and \( {F}_o^{-1} \) is the inverse CDF (quantile function) corresponding to *P*_{o}.

The empirical CDFs were approximated using empirical percentiles at a fixed interval of 0.01. Values in between the percentiles were approximated using a linear interpolation. A threshold for the correction of the number of wet days was estimated from the empirical probability of non-zero values in *P*_{o}. The correction of the daily values at each grid box was done for each month separately to account for seasonality.

The performance of QM was assessed by using a split-sampling validation approach. For KNMI-RACMO22E the historical model output comprises the years 1950–2005, which was split into the 30-year calibration periods 1950–1979 and 1976–2005. Bias correction was done for each of the two calibration periods and the performance was validated in the two, from the correction independent, 26-year periods 1980–2005 and 1950–1975, respectively. For CLMcom-CCLM4-8-17 historical model output was available for the period 1960–2005 and the 25-year calibration periods 1960–1984 and 1981–2005 were used for setting up the transfer functions between observed and modelled values. The statistical bias correction was subsequently validated in the independent 21-year periods 1985–2005 and 1960–1980, respectively. Root mean square error (RMSE) between observed and modelled values was used as performance measure. RMSE is calculated for the arithmetic means across all grid cells, assuming an equal weight for each cell. Thus, it serves as a cumulative measure of the bias over the considered domain.

### Distribution modelling and projection

Climate data (mean, minimum and maximum temperature and precipitation) of each month separately in the time period 1985–2009 as well as the land cover data served as predictors for vector occurrences. In order to quantify the relationships between vector occurrences with climate and land cover variables and to map and project the occurrences under present and future climate conditions Boosted Regression Trees (BRTs) were used. Detailed descriptions of BRT are provided by Elith et al. [14] and Hastie et al. [22]. BRT combines regression trees and boosting. BRT attempts to minimize a loss function, which involves jointly optimizing the number of trees, learning rate, and tree complexity. The learning rate is used to shrink the contribution of each tree as it is added to the model. Slowing the learning rate increases the number of trees required. In general, a smaller learning rate (and consequently a larger number of trees) is preferable. Tree complexity (number of nodes in a tree) relates to the interaction order in the predictand. With increasing tree complexity, learning rate must be decreased if sufficient trees are to be fitted to minimize predictive error. The tree complexity should reflect the correct interaction order in the response variable. However, since an adequate tree complexity is usually unknown, it is best evaluated using independent data. As in Elith et al. [14] the optimal number of trees, learning rate and tree complexity were estimated with a cross-validation approach, using deviance reduction as performance measure. The *dismo* and *gbm* packages in R were used to assess the optimal number of boosting trees using 10-fold cross-validation. In the present study models were developed with 50% of the data, and were validated with the remaining data. Tree complexity of 2 up to 8, and learning rates of 0.005, 0.1 and 0.5 were evaluated.

The modeling of the vector distributions using BRT requires both presence and absence data. The lack of confirmed absences in the occurrence data was addressed by the production of artificial absence data, called pseudo-absences. The pseudo-absences are all grid boxes outside of the suitable area, which is estimated by a rectilinear surface range envelope [23]. Following the recommendation of Barbet-Massin et al. [24] the same number of pseudo-absences as presences was tested (ratio 1:1, with 1000 presences and 1000 pseudo-absences randomly selected from the available data). Additionally, ratios of 5:1 (5000 pseudo-absences and 1000 presences) and 10:1 (5000 pseudo-absences and 500 presences) were also tested, since Sinka et al. [5], although using different predictor data and BRT setup, found the best overall performance for the European and Middle Eastern *Anopheles* species with a ratio of 10:1 pseudo-absences to presence. Model validation was subsequently done using the remaining independent data not used for model building. The BRT model was used to predict vector occurrences to the independent data and the result was taken for evaluation of the model. As statistics on predictive performance deviance, correlation, discrimination and Kappa were estimated and results were also evaluated visually. Details on cross-validation and performance measures can be found for instance in [25,26,27].

Subsequently, the best performing BRT configuration was used to project vector occurrences under future climate change. For this purpose, the bias-corrected RCM data were taken as new predictor data in the BRTs. The projected occurrences were evaluated for the historical period 1985–2005, and the two scenario periods 2040–2060 and 2080–2100.

### Potential malaria transmission stability

The vector stability index (VSI) of Kiszewski et al. [15] is used to generate maps of the future potential malaria transmission stability under climate change:

$$ \boldsymbol{VSI}=\sum \limits_{\boldsymbol{m}=\mathbf{1}}^{\mathbf{12}}{\boldsymbol{a}}_{\boldsymbol{i},\boldsymbol{m}}^{\mathbf{2}}{\boldsymbol{p}}_{\boldsymbol{i},\boldsymbol{m}}^{\boldsymbol{E}}/-\mathbf{\ln}\left({\boldsymbol{p}}_{\boldsymbol{i},\boldsymbol{m}}\right) $$

(2)

where m is month; i is vector; a is the human-biting proportion (0−1); p is the daily survival rate (0−1); and E is the length of extrinsic incubation period in days (for *P. vivax* E = 105/T-14.5). VSI was calculated for each vector i. The parameters a and p for each *Anopheles* species were taken from the publication of Kiszewski et al. [15]. Within the calculation of the length of extrinsic incubation period E, mean temperature (T) of the historical time slice 1985–2005 as well as of the future time slices 2040–2060 and 2080–2100, taken from the bias-corrected RCM data, were used. The vector-specific VSI results were integrated into an overall information by multiplying the VSI value with the modelled occurrence probability for each vector at each grid box. The VSI of that vector which has the highest value of combined VSI and occurrence probability at a particular grid box was subsequently mapped.