Analysis of meteorological effects on the incidence of tick- borne encephalitis in the Czech Republic over a thirty-year period

The goal is to analyse the effect of climatic factors on the incidence of tick-borne encephalitis (TBE) in three adjacent, high-incidence regions of the Czech Republic (CZ) that differed considerably in the baseline incidence and in the timing of high incidence period. The basis for the analysis were the meteorological data from a 30-year database (1982-2011) of daily air temperature and precipitation measurements and time-matched reports of TBE cases defined by the onset of symptoms. A semi-parametric Poisson regression model revealed a statistically significant relationship between the average daily ambient air temperature and incidence of TBE. The shapes of the dependence were similar in all three regions under study. The analysis of the available data did not show a statistically significant relationship between the precipitation amount and incidence of TBE. Other statistically significant predictors were the trend over calendar years, monthly seasonality and trendseasonality interaction (the shape of the monthly seasonality varies over years), and day of week. The climate changes had differentiated effects on TBE increase depending on the local environmental conditions as can be illustrated by the variations between the three regions not only in the TBE incidence level but also in the shape of the trend over years. Correspondence to: Marek Malý, National Institute of Public Health, Šrobárova 48, 100 42 Prague 10, Czech Republic, Tel: +42


Introduction
In the early 1990s, tick-borne encephalitis (TBE) was on the rise, to different extents, across the whole European area where it is known to spread. Similar trends were also seen in other tick-borne infections. At the same time, considerable climate changes were observed. In this context, several working hypotheses emerged about the causal links between these two phenomena [1]. A number of papers appeared in an attempt to elucidate the situation, ranging from in-depth analyses of the regional conditions to ungrounded speculations about the trigger effect of the political changes following the collapse of communism. To address the problem in its complexity, multiple factors and the extent of their involvement need to be considered [2][3][4][5][6]. These factors can be divided into the following areas: 1) changes in the ecology and spread of Ixodes ricinus tick, as the main vector of the TBE virus (TBEV) in Europe; 2) changes in land-use resulting in changed habitats and biocenoses, including ticks and their animal hosts maintaining the TBEV in the wildlife zoonotic cycle; 3) socio-economic changes influencing human behaviour and leading to more frequent and closer contact with the nature; 4) changes impacting directly the TBEV which is found in the nature as a heterogeneous population of more and less virulent variants whose proportions can vary under selection pressure. The factors listed above either can be influenced by the climatic factors or, conversely, can modify, to a certain extent, the climatic factors.
The onset of the infectious process is always associated with an infected tick bite of human (which is also true of the food-borne infections where the source of infection is an infested lactating animal). Research into the ecology of I. ricinus has brought about a body of evidence (both laboratory and field) for the crucial effect of the meteorological factors such as the air temperature and humidity on the occurrence and host-seeking activity of this vector, as described by many authors [2,[6][7][8]. From the perspective of the importance of the vital role of I. ricinus as the vector, significantly influenced by the meteorological conditions, in the epidemiology of TBE, the meteorological factors in general and ambient air temperature in particular can be considered as the driving force behind the reported changes in the incidence of TBE [9]. That is why the present study focused on the role of the meteorological factors in the epidemiology of TBE.
The objectives were as follows: 1. To analyze the effect of the meteorological factors (ambient air temperature and precipitation rate), characterized by the daily averages for each study area, on the incidence of TBE; 2. To use as the basis for the analysis a 30-year database of meteorological measurements generally accepted as adequate to determine the climatological characteristics of the study area. To use standard meteorological measurement data directly from the high incidence areas; 3. To analyze the trend over years and seasonality of TBE cases (after adjustment for the effect of the air temperature); 4. To compare the effects of the meteorological factors in three adjacent high incidence regions that differed considerably in the baseline incidence and in the timing of the high incidence period.

Meteorological data
The whole-year data on the average daily air temperature and daily precipitation amount were derived from the database of the Czech Hydrometeorological Institute (CHMI) in Prague. They originated from 26 CHMI stations provided with standard equipment and distributed evenly across the high TBE incidence areas of the regions under study. From the point measurements done at sites of the CHMI stations, the averages were calculated for the respective areas using the standard procedures taking into account the local geographical characteristics. There were selected nine CHMI stations in South Bohemian Region (average altitude above sea level of 474 m), eight stations in Central Bohemian Region (average altitude above sea level of 336 m), and nine stations in Highlands Region (average altitude above sea level of 538 m). These CHMI stations are located between geographical coordinates 49°03´40´´ -50°22´02´´ N and 13°55´27´´-16°07´14´´ E.

Statistical analysis
Daily TBE incidence data (given for each three regions separately) were modelled via semiparametric Poisson regression [10]. The Poisson model was stratified on region and took into account several biologically important explanatory variables (year, month, ambient temperature, precipitation rate, day of week). The model can be viewed as a generalization of the generalized linear model (GLM) class [11], namely as an instance of penalized-spline version of the generalized additive model (GAM) class [12,13]. The model estimation was based on the penalized likelihood and done in R (R Development Core Team, 2014) using the mgcv package [13]. The flexible semiparametric model was used to study effect of different variables upon TBE inhomogeneous Poisson intensity [14] in a detailed way not influenced by heavy a priori assumptions about the shape of the (nonlinear) relationships. The intensity was expressed as the incidence per 100000 permanent inhabitants of a region. To this end, we worked with Poisson model with canonical log link and offset. All tests were interpreted at 5% significance level.
During the data analysis, we formulated three alternative models: (1), (2), (3), labelled in the direction of decreasing complexity. Final model was selected with the using of the Akaike Information Criterion (AIC) [15] in order to balance quality of the fit and model complexity.
where N ymd is the number of TBE cases in day d of month m and year y (this is the response of interest).
µ ymd is the Poisson intensity of daily TBE incidence P y is the population of the region in year y (note that in the models (1), (2), (3), it acts as an offset) I(.) is an indicator (assuming 1 if its argument is true and 0 otherwise) β 1 ,..., β 7 are parameters estimated from data (they correspond to the effects of different levels of the day of week factor, analogous to the "textbook" parametrization in ANOVA models) T ymd,I is average ambient temperature for the given region, l days before the day ymd of the TBE incidence. In our modelling, for the models with offset we use the lag l=21 days. We experimented with other (both shorter and longer) lags. They gave generally worse AIC and worse interpretation of the estimates than l=21, however. Nevertheless, it was not our primary goal to estimate the lag precisely or to formally compare it to estimates obtained from other sources. In fact, the data do not contain too much information about precise lag and hence we pick 21 lag just as a reasonably well fitting (not necessarily globally optimal) value.
R ymd,I is average precipitation intensity for the given region, l days before the day ymd. In our modelling, we use the lag l=21 days. We experimented with other (both shorter and longer) lags. They gave generally worse AIC and worse interpretation of the estimates than l=21, however. are smooth functions estimated from data (as thin-plate splines). Roughness penalty is estimated from data via cross validation.
are smooth marginal effects of the trend over years, seasonality (in monthly resolution), (lagged) temperature, (lagged) precipitation intensity. are smooth interaction terms: trend over years -seasonality interaction (i.e., deformation of the seasonality from year to year), temperature -precipitation interaction. In both cases, the smooth interaction is modelled as a tensor-product spline.
Model (1) is the most complicated, while models (2) and (3) result from progressive simplification as special cases. All models include control for the day of week, for the trend over years and seasonal component. The models differ in the complexity of the effect of the meteorological variables. The day of week was originally considered to be a nuisance variable which is not of direct interest, but which is used to control for administrative inaccuracies when entering TBE reports, but as noted by one of the anonymous referees, it might be of some interest as it could reflect the periodicity in the recreational activities. The most complicated model (1) covers the effect of both the air temperature and precipitation, with interaction. Model (2) also covers the effect of both the air temperature and precipitation, but without the interaction (additive on the log link scale and multiplicative on the incidence scale). Model (3) only covers the effect of the air temperature. When comparing models (1) and (2), the temperature-precipitation interaction is tested and when comparing models (2) and (3), the precipitation effect is tested after excluding the interaction.
It is important to note that the model is stated so that in its systematic part, it captures both temperature (or temperature and precipitation) as the terms of primary interest and the residual part which describes inter annual changes not related to temperature (or temperature and precipitation). The residual part contains effects of other variables which might be potentially important in influencing the TBE incidence (like animal vector density, intensity of human visits to forested areas etc.) but whose data we do not have at our disposal for the period covered by the TBE incidence inventory. Term

Characteristics of the meteorological factors by region over 1982-2011
The input data from the analyses conducted to characterize the climatic background (air temperature and precipitation) are presented in Figures 1 and 2. The ambient air temperatures, despite great variability, show a slow upward trend. From the perspective of the precipitation, several-year periods of higher precipitation rates alternate with those of lower precipitation rates, but without any obvious trend. A sharp deviation in 2002 reflects the exceptional precipitation events that caused extensive floods over a great part of the Czech Republic.

Effect of the meteorological factors (average daily temperature and average daily precipitation amount) on the incidence of TBE by region
The parameters for models (1), (2), and (3) were estimated separately for each region and the models were compared. Neither the temperature-precipitation interaction in model (1) nor the effect of the precipitation in model (2) was statistically significant in any of the regions. Based on the AIC, the more complex models did not provide a clearly better fit. For example, the AIC values for South Bohemian Region were 12019.36 for model (1), 12022.13 for model (2), and 12024.40 for model (3). Therefore, concordantly in all three regions, the simplest model (3) which only takes into account, apart from the effect of seasonality, trend over years, and day of week, the effect of the air temperature was selected for the detailed analysis. It needs to be underlined that the statistical non-significance of the effect of the precipitation in models (1) and (2) definitely does not imply that the precipitation in general has no effect -we only were not able to prove its existence on the currently analysed data, unlike for the temperature effect. Table 1 shows the p-values from model (3) by region. All effects in the model are statistically significant for any region at a 0.05 significance level.  The relationships between the air temperature and TBE incidence estimated from model (3)  The effect of the average daily air temperature increases from zero to 12-14°C while higher temperatures result in saturation, with the curve plateauing or, at temperatures close to 20°C, slightly moving downwards and upwards (but the small plateau level changes are well within the estimation error). Although the input data (particularly TBE incidence data) considerably differ between regions and the models for the regions are estimated independently of each other, there is a clear concordance between the regions in the curve shapes. Differences are only seen in TBE incidence levels at higher air temperatures and mirror those in TBE incidence between the regions. The ascendant part of the curve for Central Bohemian Region has a less sharp slope in comparison with South Bohemian Region and Highlands Region. The temperature matching the peak incidence and the pattern and shape of the peak area of the curve are concordant in all three regions.
At this point, it is interesting to note significance of the day of week effect as it might reflect to some extent the weekly varying risk of the TBE (in connection with the weakly varying recreational activity), but in this data, it might be convoluted with other effects (including minor administrative inaccuracies), so that we do not discuss its details further.

Seasonality in the pattern of the daily TBE incidence within the year
As can be seen from the results in Table 1 for the effects of the day of week, monthly periodicity, and shape change of the monthly periodicity between years, the seasonal component of the daily incidence of TBE is clearly evident. Figure 4 illustrates the (certainly non-negligible) variation of the seasonal component between years in the Highlands Region after the correction for temperature. Low before 1993, the incidence of TBE slowly increased over the following years and, since 2006, a considerable and persistent rise was occurring. The curves have some characteristics in common over years, but differ in the time of peaking and pattern, ranging from unimodal in some years to bimodal in others, with a side peak appearing early in the autumn.

The trend over years in the pattern of the annual incidence of TBE
The model used also makes it possible to extract and analyze the residual trends over years in the incidence of TBE after correcting for temperature. The trends are highly statistically significant in all three regions (Table 1). Figure 5 shows clearly similar inter-annual trends (exp(S Y (y)) components) in the Central Bohemian Region and Highlands Region until 2004 with a subsequent sharp rise in the Highlands Region. The South Bohemian Region with a considerably higher incidence of TBE differs from the Central Bohemian Region and Highlands Region, on the one hand, in the incidence of TBE and, on the other hand, in the trend pattern. Unlike the Central Bohemian Region and Highlands Region, the South Bohemian Region shows a higher incidence of TBE over the whole study period, with a clear rise at the turn of the 1980s and 1990s and a certain drop in the last decade of the study period.
The above-mentioned modelled components of the monthly seasonality and trend over years plus their interaction corresponding to seasonality shape deformation from year to year (i.e., the term ) can be integrated together to yield a more comprehensive picture. The resultant estimated patterns of TBE incidence per 100,000 population by region are illustrated in Figures  6-8 using different scales on the y axis for reasons of legibility. Clear differences can be seen between the regions. The Central Bohemian Region and South Bohemian Region showed increase in TBE incidence, although to different levels, in the first half of the 1990s. In the Highlands Region, the first peak is evident in the second half of the 1990s and another more pronounced one appears since 2007. The South Bohemian Region exhibited the highest absolute number of TBE cases over the whole study period of 30 years as is clearly reflected in the pattern of the figure. In comparison with the two other regions, the seasonal component plays a greater role as reflected in the mostly bimodal patterns of the annual incidence. The hints of the sporadic incidence of TBE in the winter periods between 1995 and 2000 are also more pronounced. The patterns of the curve in the Highlands Region differ considerably between three periods which nearly match the three decades under study as already commented for Figure 4. The annual slump of the curves to zero is typical of the winter period. It needs to be pointed out that in the highest incidence years the curve shows some sporadic cases even in the winter. It means that, under suitable conditions, TBE may also emerge in this season of the year.

Discussion
The thirty-year study period was selected to span three consecutive decades with different TBE incidence: a low-incidence decade, a sharp  rise decade, and a persistent high-incidence decade (with some yearon-year variations). The end year of the study is 2011. The following years are not considered because of the possible impact of the increasing vaccination coverage of the population against TBEV as a result of the ongoing education campaigns that turned out to be effective mainly in the South Bohemian Region [16].
Great attention was paid to the climatological characterization of the study area. The input data were the measurements of the professional CHMI stations operated in a standard mode in the areas where TBE cases occurred. The criterion was the availability of a 30year continuous data series which is generally accepted as adequate to determine the climatological characteristics of the area under study. The selected CHMI stations are representative of the high-incidence areas and not of the whole administrative regions. E.g. in the South Bohemian Region, the Šumava mountains' higher altitudes where TBE cases emerge sporadically were not included in the study.
Detailed field research has demonstrated the predominant role of ambient temperature on both the activity of ticks [17] and its long-term effects on their distribution to areas previously not inhabited by them in the CZ, primarily to higher altitudes [18].
Studying the influence of abiotic factors and their changes on the epidemiology of the TBE in CZ there have been analysed in parallel also other factors that may affect the incidence of TBE human cases and was estimated the degree of such influence. Kříž, et al. [4] documented that incidence of TBE is not linked with the socio-economic condition of     the human population in the CZ. Further on, there were studied: foodborne TBE cases caused by consuming unpasteurized milk products of grazing economic animals [19], potential role of game animals in TBE epidemiology [20] and different types of recreational out-door activities influencing the TBE incidence [9]. It was not found an important role in any respects under studies. The main information is summarized by Kříž, et al. [21] Apart from temperature -other factors are considered. The very important biotic factor is human behaviour, predominantly the preventive vaccination against TBE virus, highly reducing the TBE human cases prevalence. However, TBE vaccination campaign started in the CZ only in the nineties. Recently (2013) an average of the vaccination coverage in the CZ reached only 23% [22]. The highest coverage was reported from South Bohemian Region (33.2%). There is obvious also in our results ( Figure 5).
An important step was to determine the time lag between the onset of symptoms and the beginning of the infectious process, i.e., being bitten by an infected tick under certain meteorological conditions which were crucial for the analysis. This relationship is expressed by the time lag which corresponds to the incubation period of the disease and was set to 21 days, based on statistical appraisal and expert opinion as an approximation of the real incubation length of which we have not enough information.
The study area extends over three high incidence adjacent regions exceeding the whole-country average. An important factor in this selection was also the outcome of the analysis of the local patient medical history data from 2001-2010 [16], most patients were shown to have acquired the infection in their respective regions of residence. The endemic to imported TBE cases ratios in the above-mentioned period were 1456:7 in the South Bohemian Region, 596:7 in the Highlands Region, and 685:50 in the Central Bohemian Region.
The high incidence rates of TBE in the three regions were achieved from different baseline levels, by different ways, and in different years. While the Central Bohemian Region and South Bohemian Region reported higher TBE incidence than the whole-country average as early as in 1982-1991, only sporadic TBE cases emerged in the Highlands Region at that time. Unlike other high-incidence areas in the CZ, the Highlands Region only exhibited increase in TBE cases since 1997 when the local TBE occurrence exceeded for the first time the average whole-country incidence per 100,000 population. In 2006, the incidence of TBE in the Highlands Region was even twice as high as the whole-country average, as pointed out by Danielová, et al. [23]. In that study the authors only had the meteorological data available for the analysis from a single CHMI station located at the northern boundary of the Highlands Region, but they already suggested the air temperature increase as the most probable explanation for the rise in the incidence of TBE. The results of this study based on 30-year data along with the observations from nine CHMI stations (distributed across the whole study area) confirm the above-mentioned assumption as well as the effect of the air temperature, bringing about the changes in the occurrence and host seeking activity of I. ricinus, on the incidence of TBE in the 1990s (and later).
The type of the landscape is another clue to a better understanding of the situation in the Highlands Region. Orographically, the area belongs to the Bohemian-Moravian Highlands and TBE cases emerged at altitudes close to the previously reported upper borderline of the range of I. ricinus ticks. When scarce, the local tick populations only caused sporadic cases of TBE. In the early 1990s, the local I. ricinus populations became more abundant, similarly to other areas at higher altitudes in the Czech Republic [18,24]. The consequently increased circulation of TBEV, at first among wild animal hosts, eventually led to a higher incidence of TBE in humans. Due to the long life cycle of ticks (two years on average in the CZ), there was a several-year time lag in the high incidence of TBE between the Highlands Region and the other regions. No changes were revealed in the land use and demographic or socio-economic characteristics in the Highlands Region over the decade 1997-2006 which could provide an alternative explanation for the sharp rise in the incidence of TBE [23].
In the 1990s when climate changes were observed globally, TBE cases were on the rise in the areas with active foci of TBEV, sustained by local populations of ticks that became abundant enough to maintain the TBEV in the wildlife zoonotic cycle to an extent that may pose human health risk. This was the case of the Central Bohemian Region and South Bohemian Region. In contrast, in the Highlands Region, there were latent foci, i.e., areas where the TBEV was present, but in an extent, that does not pose human health risk, as its vector I. ricinus was so scarce that it hardly sustained the circulation of the TBEV in the wildlife. (During the latency of foci, TBEV could long survive in ticks only as a result of vertical transmission without the involvement of warm-blooded hosts which may have led to partial attenuation of the virus.) First of all, the local populations of I. ricinus ticks must have become more abundant due to the climate changes in the early 1990s. The consequent activation of TBE foci was responsible for a sharp rise in TBE cases reported after a time lag.
Models used to investigate the interaction between the ambient air temperature and precipitation rate turned out to be helpful in elucidating the possible role of the climate change in the rise of TBE cases in the 1990s. A separate analysis of correlation between the precipitation rate and TBE incidence by year and [17] showed the effect of precipitation in the second half of the growing season, i.e., late summer to early autumn as stated previously [25]. Our results indicate that this effect was not associated with the air temperature. It means that the direct effect of the precipitation rate on the seasonal variations in TBE incidence in certain years is evident; nevertheless, from the long-term perspective of the three decades, the air temperature had the most important effect on the increase in the overall incidence of TBE. The effect of the precipitation or incidence -precipitation interaction in the models used was statistically non-significant in all study regions but this does not imply that the precipitation has no effect on TBE incidence at all. High resolution relative air humidity data which would probably have a stronger predictive power for the incidence of TBE than the precipitation amount were not available for the analysis.
Even after filtering out the direct effect of temperature increasing year-to-year trend in TBE incidence is obvious which may be related to the global climate changes.
It is worth noticing that the similarity of the curves in Figure 3 representing the dependence of TBE on the ambient air temperature in different regions is by no means obvious a priori and corroborates similar general influence of the temperature. This is an interesting finding opening a way to consideration of different temperature scenarios from climatic models as one of the driving forces of the TBE incidence. Potential effect of global warming can be relatively easily simulated by shifting the values along the vertical axis, using climate forecasting model outputs [26,27].