Implications of Projected Hydroclimatic Change for Tularemia Outbreaks in High-Risk Areas across Sweden

Abstract

Hydroclimatic change may affect the range of some infectious diseases, including tularemia. Previous studies have investigated associations between tularemia incidence and climate variables, with some also establishing quantitative statistical disease models based on historical data, but studies considering future climate projections are scarce. This study has used and combined hydro-climatic projection outputs from multiple global climate models (GCMs) in phase six of the Coupled Model Intercomparison Project (CMIP6), and site-specific, parameterized statistical tularemia models, which all imply some type of power-law scaling with preceding-year tularemia cases, to assess possible future trends in disease outbreaks for six counties across Sweden, known to include tularemia high-risk areas. Three radiative forcing (emissions) scenarios are considered for climate change projection until year 2100, incuding low (2.6 Wm⁻²), medium (4.5 Wm⁻²), and high (8.5 Wm⁻²) forcing. The results show highly divergent changes in future disease outbreaks among Swedish counties, depending primarily on site-specific type of the best-fit disease power-law scaling characteristics of (mostly positive, in one case negative) sub- or super-linearity. Results also show that scenarios of steeper future climate warming do not necessarily lead to steeper increase of future disease outbreaks. Along a latitudinal gradient, the likely most realistic medium climate forcing scenario indicates future disease decreases (intermittent or overall) for the relatively southern Swedish counties Örebro and Gävleborg (Ockelbo), respectively, and disease increases of considerable or high degree for the intermediate (Dalarna, Gävleborg (Ljusdal)) and more northern (Jämtland, Norrbotten; along with the more southern Värmland exception) counties, respectively.

1. Introduction

Climate change may cause shifts in geographical range, prevalence, and/or severity of some infectious diseases [1,2,3,4,5], including tularemia [6,7], a dangerous zoonotic disease caused by the intracellular bacterium Francisella tularensis [8] and widely prevalent in Europe, Asia, and America [8]. Transmission of Tularemia is usually caused by contact with infected rodents and hares, or by arthropod vectors [9]. In Europe, there is also a strong association between F. tularensis subsp. holarctica (Type B) and water conditions, with many humans reported to have contracted the disease around lakes and rivers [10]. Europe as a whole does not have a clear trend of tularemia outbreaks in recent decades, but rather a pattern of repeated local emergence and re-emergence throughout most countries [11]. In Sweden, however, the nationwide incidence of tularemia increased during the period 1984-2012 and the disease now occurs over a larger geographical area [12].

Previous studies have reported relationships between hydroclimatic factors and tularemia outbreaks [10,13,14,15,16,17], with some also evaluating future change scenarios [10,15,16]. However, perspectives and conclusions regarding future tularemia changes vary. For example, for tularemia in Sweden, Palo et al. [16] concluded that warming should not increase the frequency of tularemia outbreaks, whereas Rydén et al. [10] addressed a future scenario of an approximately 2 °C increasinge, concluding that increase in monthly summer temperature should be expected to increase the duration of tularemia outbreaks in Sweden, and Ma et al. [6] showed generally high tularemia sensitivity to hydroclimatic variability and change.

In view of the quite limited investigations so far, and their different perspectives and conclusions, this study aims to more comprehensively consider future climate change projections and assess their implications for tularemia incidence. This is done with focus on change trends in disease outbreaks along the steep climatic gradient spanned by different Swedish sites (counties) with relevant data and previously established statistical disease models, which are here combined with the latest outputs of a multi-model ensemble of global climate models (GCMs) from phase six of the Coupled Model Intercomparison Project (CMIP; [18]).

2. Methods

2.1. Six High-Risk Counties with Established Statistical Models of Tularemia

The cases considered for future disease trend projections are six Swedish counties (Norrbotten, Jämtland, Dalarna, Gävleborg, Värmland, and Örebro) distributed throughout Sweden (Figure 1). These counties fully or partly encompass seven previously identified tularemia high-risk regions that account for 56.4% of tularemia cases in Sweden, even though they contain only 9.3% of Sweden’s population, according to disease surveillance data for 1984–2012 [14]. Moreover, projection at county scale is consistent with the administration system in Sweden, where reported diseases are managed per county.

Previously established statistical models (fitted to data) have been reported and are of the same type but with different coefficients for each Swedish site [14]. Statistical models do not assume known functional relationships between infectious disease spreading and hydroclimatic and other environmental variables, but quantify the data-given statistical signals of such relationships so far [19]. These can further be used to quantify the implications of GCM-projected scenarios of future hydroclimatic changes for diseases like tularemia, for which adequate mechanistic knowledge is still lacking, e.g., on bacteria ecology and transmission routes to humans [10]. Figure 2 illustrates schematically the overall approach to such quantification developed and applied in this study.

The disease models used in this study were selected because they are, so far, the only quantitative models available that have been developed and tested for modeling of temporal variations of tularemia cases in direct relation to associated hydroclimatic conditions, based on actual data for targeted high-risk regions in Sweden from 1984 to 2012 [14] . The basic disease model type is power-law scaling, which relates the annual number of tularemia cases (Tul) to the preceding-year number of cases (Tul_lag) with exponent ${\beta }_1$, as Equation (1) [14,17]:

$$\mathrm{Tul}=\mathrm{EXP}\left({\beta }_0+{\beta }_2\mathrm{sRMA}+{\beta }_3{\mathrm{ST}}_{\mathrm{lag}}+{\beta }_4\mathrm{SP}\right)\ast {\mathrm{Tul}}_{\mathrm{lag}}{}^{{\beta }_1}=\mathrm{A}\ast {\mathrm{Tul}}_{\mathrm{lag}}{}^{{\beta }_1}$$

(1)

The scale factor A in Equation (1) is determined by disease-independent hydroclimatic variables of: summer temperature in the preceding year (ST_lag, °C); summer precipitation in the same year (SP, mm); along with standardized relative annual mosquito abundance (sRMA). The latter, derived from annual mosquito aboundance(RMA), in turn fully depends on hydroclimatic variables, as expressed by Equations (2)–(5) [14,17]:

$$\mathrm{sRMA}=\frac{{\log }_2\mathrm{RMA}-\mathrm{mean}\left({\log }_2\mathrm{RMA}\right)}{\mathrm{SD}\left({\log }_2\mathrm{RMA}\right)}$$

(2)

$$\mathrm{RMA}=\mathrm{Median}\left({\mathrm{RMA}}_{\left(\mathrm{t}\right)}\right)$$

(3)

$${\mathrm{RMA}}_{\left(\mathrm{t}\right)}=2^{{\mathrm{S}}_{\mathrm{N}\left(\mathrm{F}\right)}\mathrm{SM}}$$

(4)

$$\mathrm{SM}=-2.76+0.67{\mathrm{Q}}_1+0.62{\mathrm{Q}}_2+0.19\mathrm{T}$$

(5)

where mean and SD (standard deviation) denotes the sample mean and sample standard deviation of RMA; RMA is median value of daily mosquito aboundance ${\mathrm{RMA}}_{\left(\mathrm{t}\right)}$; $\mathrm{SM}$ is the standardized mosquito abundance; ${\mathrm{S}}_{\mathrm{N}\left(\mathrm{F}\right)}$ is the standard deviation of $\mathrm{S}\mathrm{M}$ observed in the mosquito modeling area, which is cancelled out when substituted into Equation (2); Q1 and Q2 (unitless) are maximum standardized river flows in two time periods preceding each RMA evaluation time t (36–42 days and 22–28 days, respectively, before time t), and mean temperature (T, °C) over 1–7 days before time t. Original model expressions also included winter days with low snow coverage (CW, days) in A, but the variation in Tul has been shown to be insensitive to CW [6]; as such, it has been omitted from this analysis, for simplicity and increased clarity. Mosquitoes are main disease vectors in the large boreal forest regions of Alaska, Sweden, Finland, and Russia [20,21,22], and the statistical disease models used in this study include mosquito abundance as a main variable; in a fitted statistical model, however, such a variable may also be a proxy for other vectors.

Table 1. Tularemia model coefficients assigned to the six Swedish study counties.

High-Risk County	$\mathbf{Intercept}\left({\mathit{\beta }}_0\right)$	${\mathbf{Tul}}_{\mathbf{lag}}\left({\mathit{\beta }}_1\right)$	$\mathbf{sRMA}\left({\mathit{\beta }}_2\right)$	${\mathbf{ST}}_{\mathbf{lag}}\left({\mathit{\beta }}_3\right)$	$\mathbf{SP}\left({\mathit{\beta }}_4\right)$
Norrbotten	−5.73	1.16	0.20	0.35	0.005
Jämtland	−11.47	0.93	0.75	0.84	0.003
Gävleborg (Ockelbo)	−2.86	−0.19	0.36	0.10	0.010
Gävleborg (Ljusdal)	−7.60	0.09	0.29	0.50	0.009
Dalarna	−10.74	0.37	0.43	0.67	0.008
Värmland	10.11	0.99	0.82	−0.39	−0.014
Örebro	−9.19	0.73	0.75	0.18	0.023

$\beta$: coefficients in the power-law scaling Equation (1); Tul_lag: the annual number of tularemia cases to the preceding-year number of cases; sRMA: standardized relative annual mosquito abundance; ST_lag: summer temperature in the preceding year; SP: summer precipitation in the same year.

lists the ${\beta }_0$–${\beta }_4$ coefficients in the power-law scaling Equation (1), as calculated and assigned to each Swedish county based on reported best fits to outbreak data for the associated high-risk sites [14,17]. For a high-risk site extending over more than one county, the associated site coefficient values were allocated to the county containing the largest proportion of the high-risk site. Gävleborg county contains two high-risk sites (Ockelbo, southern part of the county, and Ljusdal, northern part of the county) and was assigned two comparative sets of coefficients.

2.2. Climate Change Scenario Data, Human Tularemia Data, and Projection of Annual Cases

We used daily data for required hydroclimatic variables over 2015–2100 from GCMs that provide such outputs in phase six of the Coupled Model Intercomparison Project, CMIP6 (with raw output data downloaded from https://esgf-node.llnl.gov/projects/cmip6/). The required variables are air temperature T (corresponding to GCM variable: tas), precipitation P (pr), and runoff R (mrro; which in turn relates to river flow Q as R times the (constant) river catchment area). Because the finest available spatial resolution for daily runoff R is 100 km, we chose this spatial resolution also for the other hydroclimatic variables (T and P).

With regard to climate projection scenarios, we considered the low, medium, and high end scenarios of Shared Socioeconomic Pathways (SSPs) from the CMIP6 range of future radiation and emission pathways, i.e., SSP1-2.6, SSP2-4.5, and SSP5-8.5, respectively [23]. The SSP1-2.6 scenario represents low emissions driven by sustainable practices to produce radiative forcing of 2.6 Wm⁻² by 2100, leading to a multi-model mean warming projection of significantly less than 2 °C warming by 2100, while the SSP5-8.5 scenario represents sufficiently high emissions to produce radiative forcing of 8.5 Wm⁻² by 2100, leading to projected warming of 4.9 °C by 2100. The scenario SSP2-4.5 represents an intermediate radiative forcing level following continued historical patterns, leading to likely more realistic model projections than the other two more extreme scenarios. As such, we show in the following main results for the SSP2-4.5 climate scenario, and comparative results from the other two scenarios in Supplementary Material.

For each climate scenario, we only considered GCMs with daily outputs of all hydroclimatic variables required for full ${\beta }_0$–${\beta }_4$ quantification, as listed in

Table 2. List of global climate models (GCMs) considered in this study.

Model	SSP1-2.6	SSP2-4.5	SSP5-8.5
BCC-CSM2-MR	X	X	X
MRI-ESM2-0	X	X	X
NorESM2-MM	X	X	X
EC-Earth3	X	X	X
* INM-CM5-0	X	X	X
* INM-CM4-8	X	X	X
MPI-ESM1-2-HR	X	X	X
GFDL-ESM4	X	X	X
GFDL-CM4		X	X

* Not included in Gävleborg county because of data error. Remarks: BCC-CSM2-MR: The Beijing Climate Center Climate System Model, Version 2-MR; MRI-ESM2-0: The Meteorological Research Institute Earth System Model, Version 2.0; NorESM2-MM: The Norwegian Earth; System Model, Version 2-MM; EC-Earth3: The European Consortium Earth System Model, Version 3; INM; CM5-0: The Institute for Numerical Mathematics Climate Model, version 5.0; INM-CM4-8: The Institute for; Numerical Mathematics Climate Model, Version 4.8; MPI-ESM1-2-HR: Max Planck Institute Earth System Model, Version 1.2-HR; GFDL-ESM4: The Geophysical Fluid Dynamics Laboratory’s Earth System Model version 4; GFDL-CM4: The Geophysical Fluid Dynamics Laboratory’s climate model version 4. SSP1-2.6: Shared Socioeconomic Pathways (SSP) with emissions driven by sustainable practices to produce radiative forcing of 2.6 Wm⁻² by 2100; SSP2-4.5: Shared Socioeconomic Pathways with an intermediate radiative forcing level following continued historical patterns of 4.5 Wm⁻² by 2100; SSP5-8.5: Shared Socioeconomic Pathways with sufficiently high emissions to produce radiative forcing of 8.5 Wm⁻² by 2100.

.

According to the overall study approach illustrated in Figure 2, we selected for each target county the relevant daily outputs for GCM grid cells with at least 40% of their area located within the county (35% for Örebro county in the Geophysical Fluid Dynamics Laboratory’s (GFDL) model GFDL-ESM4 and GFDL-CM4). We further spatially averaged the values of each output variable in these grid cells to represent a corresponding county-average daily value, and further quantified the associated Tul result for each GCM from the disease model Equation (1). Finally, we averaged the Tul results across all GCMs to determine and illustrate GCM-ensemble mean annual Tul values, and associated inter-model standard deviations (as a GCM uncertaintly measure) around the ensemble mean, for each county and for each considered climate scenario. In addition, we also determined ensemble mean values of each GCM-projected hydroclimatic and related disease variable included in the basic tularemia model Equation (1), in order to also separately illustrate the projected temporal evolutions of these county-average variables in each projected climate scenario. Initial values of Tul_lag for year 2015 were determined from data on human tularemia outbreaks for each county obtained from the data repository of the Nordic project CLINF [24] as listed in Supplementary Material Table S1.

2.3. Typology of Model Behavior Under Mean Climate State in the Long-Term

For any projected combination of hydroclimatic variable values determining the scale factor A, the power-law scaling of $\mathrm{Tul}=\mathrm{A}\ast {\mathrm{Tul}}_{\mathrm{lag}}{}^{{\beta }_1}$ in Equation (1), with ${\beta }_1$ values given for the different Swedish counties in Table 1, follows the distinct types of behaviour illustrated in Figure 3a. Specifically, blue and green curves represent sublinear ($0<{\beta }_1<$1; most county cases, Table 1) and superlinear (${\beta }_1$ > 1; Norrbotten case, Table 1) conditions of positive ${\beta }_1$, respectively, and the orange curve represents sublinear negative ${\beta }_1$ conditions ($-1<{\beta }_1<$0; Gävleborg (Ockelbo) case, Table 1). All three (power-law scaling) types of curves intersect some point of the black 1:1 line, which in turn represents (linear ${\beta }_1\approx 1$ conditions with) unchanging number of tularemia outbreaks over time (Tul = Tul_lag).

For the most common sublinear positive case ($0<{\beta }_1<1)$, the power-law curve (blue in Figure 3a) intersects the 1:1 line at ${\mathrm{N}}_{\mathrm{A}}^{\ast }$, and implies convergence of Tul over time to this ${\mathrm{N}}_{\mathrm{A}}^{\ast }$ level, for any shift in long-term average hydroclimatic conditions (i.e., shift in A) and initial ${\mathrm{Tul}}_{\mathrm{lag}}$ level larger (N₀₁) or smaller (N₀₂) than ${\mathrm{N}}_{\mathrm{A}}^{\ast }$ (Figure 3b1); we refer to ${\mathrm{N}}_{\mathrm{A}}^{\ast }$ as the endemic convergence level, in consistency with the definition and use of this term also in Ma et al. (2019). In contrast, for the superlinear positive case (${\beta }_1>1)$, the power-law curve (green in Figure 3a) intersects at ${\mathrm{D}}_{\mathrm{L}}^{\ast }$, and implies divergence of Tul over time to increasingly greater or smaller values than ${\mathrm{D}}_{\mathrm{L}}^{\ast }$, for any hydroclimatic shift in A and initial ${\mathrm{Tul}}_{\mathrm{lag}}$ level larger (N₀₁) or smaller (N₀₂) than ${\mathrm{D}}_{\mathrm{L}}^{\ast }$, respectively (Figure 3b2); we refer to ${\mathrm{D}}_{\mathrm{L}}^{\ast }$ as a divergence level for future tularemia outbreaks. For the sublinear negative case ($-1<{\beta }_1<0$), the power-law curve (orange in Figure 3a) intersects at ${\mathrm{N}}_{\mathrm{O}}^{\ast }$, and implies convergent oscillation of Tul (i.e., with decreasing oscillation amplitude) around ${\mathrm{N}}_{\mathrm{O}}^{\ast }$ over time, for any hydroclimatic shift in A and initial ${\mathrm{Tul}}_{\mathrm{lag}}$ level larger (N₀₁) or smaller (N₀₂) than ${\mathrm{N}}_{\mathrm{O}}^{\ast }$ (Figure 3b3); we refer to ${\mathrm{N}}_{\mathrm{O}}^{\ast }$ as endemic oscillation level. There is no such case among the Swedish study counties but, for the sake of completeness, we note that a superlinear negative case (${\beta }_1<-1$) would imply divergent oscillation of Tul (i.e., with increasing oscillation amplitude) around a constant level analogous to ${\mathrm{N}}_{\mathrm{O}}^{\ast }$ for the sublinear negative case. Endemic convergence (for $0<{\beta }_1<1)$ and endemic oscillation (for $-1<{\beta }_1<0$) both imply an expected constant level toward/around which tularemia cases will tend to converge with time after a shift in long-term average hydroclimatic conditions (shift in A), so we combined these two terms into simply referring to the endemic level expected for any (positive or negative) sublinear power-law scaling of Tul with Tul_lag.

3. Results

3.1. Change Trends for Hydroclimatic and Related Mosquito Abundance Variables

Figure 4 shows Dalarna county as an example for the resulting GCM ensemble mean time series (solid lines) and corresponding inter-GCM ± 1 standard deviation (shaded area) of (a) yearly summer temperature (ST), (b) summer precipitation (SP), (c) standardized relative mosquito abundance (sRMA), (d) mean temperature (T) over 1–7 days before time t of the sRMA evaluation, along with maximum standardized river flows (e) Q₁ and (f) Q₂ over 36–42 days and 22–28 days before t, respectively, for climate scenario SSP2-4.5. Results are shown in terms of 20-year running averages in (a–c), and in terms of annual average values in (d–f). Black lines show average values over the time periods 2015–2057 and 2058–2100. Supplementary Figures S1 and S2 show results for all counties corresponding to panels (a–c) and (d–f), respectively. For comparison, corresponding results of T, SP, and sRMA for the other two climate scenarios are also shown in Supplementary Material Figure S3, with the 43-year averages for the two climate periods 2015–2057 and 2058–2100 also shown in histograms in Supplementary Material Figure S4.

For climate scenario SSP2-4.5, all counties exhibit similar weakly increasing trends for ST (Figure 4a, Figure S1a), with scenario SSP1-2.6 yielding even weaker and scenario SSP5-8.5 leading to steeper increase trends (Figure S3a). For SP, change trends are overall small for all counties and climate scenarios (Figure 4b, Figure S1b, Figure S3b). For sRMA in scenario SSP2-4.5, all counties have similar subtle increase trends (Figure 4c, Figure S1c), which also applies to the other two climate scenarios (Figure S3c), with the exception of Örebro county for SSP1-2.6. Overall, similar inter-GCM uncertainty levels are exhibited for these three hydroclimatically dependent variables among counties and climate scenarios. Figure 4d–f and Figure S2 further show projected ensemble mean results for scenario SSP2-4.5 of the three hydroclimatic variables underlying and determining sRMA from 2015 to 2100. Values of temperature T generally increase while values of standardised flow Q₁ and Q₂ exhibit mainly small-range fluctuations over time.

3.2. Tularemia Change Trends Under Projected Future Hydroclimatic Conditions

Figure 5 shows resulting projected ensemble means of annual tularemia cases for scenario SSP2-4.5 for each county, along with associated endemic levels (Insert figures; expected divergence level for Norrbotten is shown in Supplementary Material Figure S5) over two climate periods (2015–2057 and 2058–2100). Comparison with the other two climate scenarios is shown in Supplementary Material Figure S5. All values shown are normalized with the respective actual numbers of tularemia cases occurring in year 2015 in each county (listed explicitly in Supplementary Table S1). The normalization focuses on and facilitates direct comparability of overall change trends and associated uncertainty levels for and across different counties. Furthermore, it avoids giving a misleading impression of predicting absolute numbers of disease cases at specific future points in time. Even though change trends are relatively mild and largely similar among counties for the hydroclimatic variables determing the scale factor A in the disease Equation (1), the counties exhibit widely diverging resulting tularemia change patterns, with consistent trends in endemic levels.

Annual disease cases in Norrbotten are projected to surge towards much higher levels if no counter-measures are taken. Large and rapid disease increases also emerge for Värmland and Jämtland. In contrast, projected annual cases in Dalarna and Gävleborg (Ljusdal) show overall slower and smaller increase trends, while Gävleborg (Ockelbo) and Örebro exhibit overall and intermittent disease decreases (values <1) from the 2015 disease status.

Comparison with the other two climate scenarios (Figure S5) shows similar dramatic rise trends for Norrbotten, Värmland, and Jämtland, while Dalarna and Gävleborg (both sites) exhibit increase in endemic levels with higher scenario forcing level. Örebro exhibits mixed annual case occurrence among climate scenarios and lowest endemic level for the intermediate scenario SSP2-4.5, but overall projected endemic levels for Örebro are <1 (i.e., smaller than the 2015 number of cases) and near-zero across all scenarios. Inter-model uncertainties can be large, and are particularly dramatic for Norrbotten, Jämtland and Värmland, while they are relatively small for Gävleborg (Ockelbo).

4. Discussion and Conclusions

Results for the six Swedish counties show large tularemia sensitivity to relatively small hydroclimatic change trends, and large inter-GCM uncertainty levels for disease projections compared to those for the underlying hydroclimatic variables. High sensitivity to the power-law disease scaling characteristics is also evident in the widely different disease projection results under more or less similar hydroclimatic change trends among the counties. Among counties, the relatively southern counties of Örebro and Gävleborg (Ockelbo) exhibit periodic or overall tularemia declines, respectively, while the most northern counties Norrbotten and Jämtland, but also the southern Värmland, exhibit large increases, and the intermediate Dalarna and Gävleborg (Ljusdal) exhibit intermediate trends of mostly increases (depending on climate scenario) until 2100.

In general, projected long-term trends in endemic levels are closely related to the power-law exponent ${\beta }_1$ for scaling of Tul with Tul_lag. With Norrbotten having superlinear ${\beta }_1>1$ and Värmland and Jämtland having sub- but near-linear ${\beta }_1$ = 0.99 and 0.93, respectively, associated tularemia and endemic level impacts tend to be enhanced by projected shifts in hydroclimatic (scale factor A) conditions. The differences in projected tularemia cases among counties and some mixed results for the three climate scenarios challenge possible notions that climate change (and higher emission scenarios) will generally lead to (higher) increases of disease incidence. The different best-fit ${\beta }_1$ values for the different counties may then implicitly reflect geographic differences in, e.g.,: (1) demographics; (2) risk of pathogen exposure; (3) other local conditions/measures affecting vulnerabity/resilience to disease; (4) and interactions among hydroclimatic factors. This notion has also been discussed in other research using statistical disease modeling to project future disease burden under various climate scenarios. For example, Hales et al. (2002) [25] estimated that climate change would lead to 50–60% of the global population being at risk of dengue transmission by 2085, compared with 35% without the projected climate change. For falciparum malaria, however, Rogers and Randolph (2000) [26] projected that, by 2050, 23 million hosts would be gained in previously uninfected regions while 25 million human hosts would be lost in areas no longer suitable for transmission, which would lead to little net disease change in total.

With regard to the latter, the high-emissions scenario (SSP5-8.5) led to similar increases in summer temperature across the counties, but warming-related changes in precipitation and runoff affected tularemia results most. Other studies have found that human-related factors may play a more important disease role than climate, e.g., superior healthcare infrastructure might lead to net lowering of disease impacts even if climate change enhances pathogen ranges [19]. The history of widely studied diseases (e.g., malaria, yellow fever, dengue fever) also shows that human activities and their impacts on local ecology may affect disease spreading more significantly than climate change [27]. In addition to external drivers, internal complexity of climate-disease interactions also affects disease risk. For example, an increase in temperature may increase mosquito biting rates, parasite replication within mosquitoes, and mosquito development, but also increase mosquito mortality, making disease outcomes difficult to determine [28].

This study also has several limitations. First, in Sweden, high tularemia incidence usually appears in low- population areas [14]. So, although an area can have high outbreak rate projections (such as Norrbotten and Värmland in this study), local population levels may set lower upper limits for outbreaks. Second, the considered disease models do not explicitly take human behavioral factors into account, such as time spent on outdoor activities, which is an important factor for exposure of humans to the disease. Moreover, the statistical tularemia models have been rather inaccurate in estimating the magnitude of recent outbreaks, especially for the two most northern counties (Norrbotten and Jämtland) [14]. Thus, the present results for how hydroclimatic changes may impact future outbreaks should be used with caution, as comparative indications rather than in terms of absolute disease outbreak projections.

Uncertainties are inevitable in and among projected results of different climate models, and in all studies using climate model outputs in other types of models. The latter, propagated type of climate-related uncertainties are here shown to be amplified in disease models with high climate sensitivity, such as the tularemia models for Norrbotten, Jämtland, and Värmland in this study. Bring et al. (2019) [29] found relatively good model-data agreement for ensemble mean outputs of runoff and temperature in the Nordic-Arctic region, but worse agreement for precipitation outputs, and considerable doubt still remains about realism and accuracy of hydroclimatic results from individual GCMs. Contradictions inevitably emerge in disease projections for localized transmission routes [14] and future work needs to continue exploration of opportunities to improve projection realism and accuracy.

Well-archived data of infectious diseases in clinics and laboratories, along with adequately-recorded climate, hydrological, and other environmental as well as socio-economic data in recent years has made it feasible to develop statistical disease models, which can further be used in combination with related projected hydroclimatic and other types of data to quantify scenarios of possible future disease evolution (Figure 2). The available historical records along with forthcoming new scenario-projected model data and model-coupling methods will surely benefit more accurate large-scale assessments of future disease pressures and risks. In addition, advancements in mechanistic disease modeling are needed to bridge gaps and overcome weakness of statistical models, which may not be as relevant for other locations and new hydroclimatic and other environmental and societal conditions than the ones they were fitted to.

In conclusion, this study has quantified the implications of scenario-projected future hydroclimatic trends for possible future disease evolution, using site-specific, established, and parameterized statistical tularemia models. Results show highly divergent disease change trends and fluctuation levels around these for future climate change scenarios among Swedish counties, with scenarios of steeper future climate warming not necessarily leading to steeper disease increases. The directions of future tularemia change trends are robust in some counties, as seen from results across various future climate scenarios and their representations by different GCMs. Such robust change-trend projections are essential to identify and useful in pointing out needs for policy and management measures to avoid clear negative directions of future disease evolution, even though uncertainties about absolute future disease numbers may be large.