The UK’s suitability for Aedes albopictus in current and future climates

The Asian tiger mosquito Aedes albopictus is able to transmit various pathogens to humans and animals and it has already caused minor outbreaks of dengue and chikungunya in southern Europe. Alarmingly, it is spreading northwards and its eggs have been found in the UK in 2016 and 2017. Climate-driven models can help to analyse whether this originally subtropical species could become established in northern Europe. But so far, these models have not considered the impact of the diurnal temperature range (DTR) experienced by mosquitoes in the field. Here, we describe a dynamical model for the life cycle of Ae. albopictus, taking into account the DTR, rainfall, photoperiod and human population density. We develop a new metric for habitat suitability and drive our model with different climate data sets to analyse the UK’s suitability for this species. For now, most of the UK seems to be rather unsuitable, except for some densely populated and high importation risk areas in southeast England. But this picture changes in the next 50 years: future scenarios suggest that Ae. albopictus could become established over almost all of England and Wales, indicating the need for continued mosquito surveillance.


Development and fecundity rates
According to Quinn [72], development rates do not continuously increase with increasing temperature but show an optimum after which they decrease again. We thus choose a quadratic polynomial function d(T ) = a · T 2 + b · T + c for development periods as shown in Figure S1, and took the inverse to get development rates, δ(T ) = 1/d(T ). Data for egg development did not show a clear relationship with temperature, so we chose the median of all data points as a constant. Egg laying rates were averaged over the whole lifetime of a mature female, excluding the first pre-bloodmeal period. A skewed Gaussian distribution was fitted to the data.

Survival rates
Fitting functions to daily survival rates derived from laboratory data were of the form s(T ) = a·exp(−( T −b c ) 6 ) (Fig. S2). We use a sextic rather than a quadratic exponent as data suggest an optimal plateau rather than a peak around an optimal temperature. The daily survival probabilities are subsequently translated into mortality rates suitable for a differential equations setting by taking the negative natural logarithm, µ(T ) = − ln(s(T )). Figure S2: Daily survival probabilities for summer eggs, juveniles, and adults and winter survival probability for diapausing eggs. Data points from Chang et al. [95] are only for L1 and L4 larvae.
For daily survival of adult mosquitoes, we only use the meta analysis of Brady et al. [98] as they investigate survival in the field which we are interested in. The field survival rates are calculated for daily or even monthly mean temperatures [99,100,101,102,103], so we decided to use daily mean temperatures to calculate adult mortality rates, rather than the full diurnal temperature cycle. We used a temperature-dependent factor of T 0.1 to account for the skewness in the daily survival probability, µ A (T ) = − ln(s A (T ) · T 0.1 ). Finding a temperature-dependent curve for diapausing egg survival was the most challenging task as there are no laboratory studies investigating the survival at constant temperatures for periods of up to five months. However, some studies analyse the effect of sub-zero treatments on the eggs' survival rates [104,105,106]. They suggest that hatching success does not change significantly after a 24h exposure to -2°C in the F1 generation [106] and another study found the length of the egg maturation period does not greatly effect survival rates [96]. We thus make the assumption that only the absolute minimum temperature experienced throughout the winter is shaping the survival curve. This survival probability is then applied only when eggs are hatching in spring. This is in line with the findings of [107] whose fitted daily winter mortality was too high to give meaningful results. In the end, we used a curve fitted to survival data of Italian eggs against minimum temperatures [105] (Fig. S2). There are some studies looking into the interaction of temperature and larval densities on mortality of Aedes mosquitoes but they do not show a clear pattern [108]. We thus keep temperature-related and density-related mortalities separated.

Diapause
The relationship between latitude, L, and the critical photo period CPP A after which half of the eggs laid go into a diapausing state shows a positive correlation [109]. We use the observations made in Japan and in the US after adaptation to derive a linear function for our model, see Figure S3. Figure S3: Relationship between latitude, L, and the critical photo period CPP A in Japan and the US after adaptation [109].
The hatching of diapausing eggs in spring is triggered by a critical temperature, CTT S = 11°C and photo period, CPP S = 11.25 hours [110]. Eggs then only hatch if mean temperatures over a period of 7 days reach this threshold [111].

SI.2 Model with 7 stages
Other studies use larvae and pupae instead of juveniles to model immature development, some studies also use host-seeking and gestating stages to model the gonotrophic cycle [112,113,111]. An advantage of additional stages is that the development rate from egg to adult can approximate a Gamma distribution, which is more realistic with delayed development (linear chain trick) [114]. We tested whether such an extended model structure with additional parameters (Fig. S5) improved the fit to observed container index (CI) data but it gave very similar results. Comparing model output with spatial observations failed to predict occurrences in parts of East Europe and north of 45°N. We thus chose the model with only five equations, as there is also more literature data available to parametrise the reduced model.

SI.3 Suitability index E 0
We define the suitability of a certain year and location as the number of diapausing eggs at the end of the year divided by the (small) number of diapausing eggs at the beginning of the year (E i = E d (x = 365)/E d (x = 1)). In our model, diapausing eggs do not experience a daily mortality (a survival probability is applied in spring instead, see SI.1) and thus accumulate over autumn. At the end of December, when almost all adult females have died off and no more eggs are produced, diapausing egg numbers E d have reached a stable plateau and E d (x = 365) can be used for E 0 calculation without taking an average over a certain period.
It would be also possible to use other stages for comparison, e.g. the number of adult females during the summer. It would be necessary though, to use an average over a certain time period for this variant, as numbers of juvenile or adults can fluctuate on a daily basis. As an example, we define as the averaged number of females over 30 days in July/August of one year divided by the averaged number of females over the same period in the previous year. Figure Table S1. Values for A 0 are slightly higher as relatively more females had been introduced than diapausing eggs. The suitability definition based on diapausing eggs has the advantage that simulations only have to be run with single calendar years (often easier with climate data sets) and that population numbers do not have to The definition of E 0 as the geometric mean of all E i analysed has another advantage; E n 0 = Π n i=1 E i gives an approximate value of population increase/decrease after n years; as an example let E 1 = 3, E 2 = 0.5 and E 3 = 1.5, then E 0 = 1.3 and thus slightly suitable, and E 3 0 = 2.3. Starting with one diapausing egg at the beginning of year 1, we can expect approximately 23 eggs at the end of year 3. This estimate does only hold in the first few years though, later on the mosquito population will approximate its carrying capacity and yearly egg numbers will more or less stagnate.
with DOY being day of year, ranging from 1 to 365, and L for latitude in degrees. We approximate the time of sunrise, t s , by ).
As the dynamical life cycle model for Ae. albopictus is computationally expensive to run, we took a subset of five GCMs out of 21, while keeping those models that represent the variability in the ensemble, i.e. GCMs that maximize uncertainty. For this, we first calculated land annual temperature indices for all GCMs and each emission scenario. We then kept the GCMs that were closest to the minimum, 25th percentile, median, 75th percentile and maximum temperature increase in the ensemble of the 21 GCMs (using the RMSE that minimized the distance to each categories, Tab. S2). Please see [116] for full details.  Figure S7 shows the suitability for 2016 only. The suitability over Western Europe (UK, France, Benelux, Germany) is higher than the 10-year mean, compare Figure 2 in the main text. Note that the used climate data set was incomplete for 2016, so that northern Italy and areas around the Black Sea appear grey.

SI.7 Comparison with temporal container index data
Italian container index (CI) data of Ae. albopictus eggs are taken for a model validation in time. The data for 2011 -2015 is available from http://www.zanzaratigreonline.it/ [117]. The region of Emilia-Romagna has nine districts and each district is monitored for Ae. albopictus eggs by ovitraps that are checked and numbers are counted on a biweekly basis. To exclude migration effects, we first checked if CI data was spatially correlated, which it was not, see Figure S8. Figure S8: Test for spatial correlation Pearson's correlation coefficient between the nine regions' container index data plotted against distance between regions. Though the correlation between sampling sites is generally high, the slight negative trend between distance and Pearson correlation is not significant.
We run our model with Emilia-Romagna climate data for the period 2010 -2015 to have one year of spin-up (cell codes 01421, 01867, 01699, 01138, 00369, 00774, 01983, 00977, 02231 for the according region). Figure  S9 shows the comparison between the empirical CI data and scaled simulated egg data over a period of five years in the according region.

SI.8 DTR vs. mean temperature runs
Running the model with and without DTR enables us to compare the suitability outputs. Figure S11 shows that the northern European regions (England, northern France, Germany, Czech Republic, Scandinavia) show an increased suitability in the model runs that consider the effect of DTR. Spain and the regions around Israel show decreased suitability, while the Balkan region is indifferent. The E-OBS rainfall data set is incomplete for northern Italy and north-east Africa in many years, so the high values in these regions might be not reliable. Figure S11: Suitability of the DTR run for 2006 -2016, relative to the run using only daily mean temperatures. The red contour line indicates a value of 1, above which the DTR runs gives more suitable results, and below which the mean temperature runs do.
Including the DTR in simulations affects various parameters involved in population growth. One of the key parameters is the temperature-dependent development rate of juveniles, δ J (T ). This development rate is non-linear, with a maximum at around 30°C, compare 1/δ J in Figure S1. Due to the non-linearity, being 8°C above T mean has a bigger positive effect on the mean development rate than being 8°C below T mean , if T mean is relatively low. If T mean is high, e.g. 30°C, then any divergence from this optimal development temperature only decreases the rate. As a result, including a DTR in simulations can have positive or negative effects on population growth, depending on T mean . Figure S12 shows the development of egg numbers over 365 days at constant and varying temperatures. For Figure 4 in the main text, we e.g. divide E 18°C,DT R=8 (x = 365) by E 18°C,DT R=0 (x = 365) to indicate relative population growth at T mean = 18°C and DTR=8°C. The fact that the population size has reached the carrying capacity at 28°C after 365 days but not quite at 18°C only has a small quantitative but not qualitative effect on the relative population growth.

SI.9 Equilibria & Stability
We analysed the system of differential equations for equilibria and their according stabilities. Assuming rainfall is not a limiting factor (P R > 0 and constant), we get a trivial equilibrium, (E * , J * , I * , A * , E * d ) = 0, and a non-trivial equilibrium: Checking the eigenvalues of the Jacobian shows that the trivial equilibrium is unstable between 13°C and 32°C, meaning that we have a positive equilibrium in this temperature range for any positive carrying capacity K. K is non-zero if rainfall or the human population density are non-zero. The optimum temperature that maximises population numbers is at 27°C, close to what Brady et al. [98] have found for a maximum adult life span.