 Research article
 Open Access
 Published:
Tickborne encephalitis (TBE) cases are not random: explaining trend, low and highfrequency oscillations based on the Austrian TBE time series
BMC Infectious Diseases volume 20, Article number: 448 (2020)
Abstract
Background
Why human tickborne encephalitis (TBE) cases differ from year to year, in some years more 100%, has not been clarified, yet. The cause of the increasing or decreasing trends is also controversial. Austria is the only country in Europe where a 40year TBE time series and an official vaccine coverage time series are available to investigate these open questions.
Methods
A series of generalized linear models (GLMs) has been developed to identify demographic and environmental factors associated with the trend and the oscillations of the TBE time series. Both the observed and the predicted TBE time series were subjected to spectral analysis. The resulting power spectra indicate which predictors are responsible for the trend, the highfrequency and the lowfrequency oscillations, and with which explained variance they contribute to the TBE oscillations.
Results
The increasing trend can be associated with the demography of the increasing human population. The responsible GLM explains 12% of the variance of the TBE time series. The lowfrequency oscillations (10 years) are associated with the decadal changes of the largescale climate in Central Europe. These are well described by the socalled Scandinavian index. This 10year oscillation cycle is reinforced by the socioeconomic predictor net migration. Considering the net migration and the Scandinavian index increases the explained variance of the GLM to 44%. The highfrequency oscillations (2–3 years) are associated with fluctuations of the natural TBE transmission cycle between small mammals and ticks, which are driven by beech fructification. Considering also fructification 2 years prior explains 64% of the variance of the TBE time series. Additionally, annual sunshine duration as predictor for the human outdoor activity increases the explained variance to 70%.
Conclusions
The GLMs presented here provide the basis for annual TBE forecasts, which were mainly determined by beech fructification. A total of 3 of the 5 years with full fructification, resulting in high TBE case numbers 2 years later, occurred after 2010. The effects of climate change are therefore not visible through a direct correlation of the TBE cases with rising temperatures, but indirectly via the increased frequency of mast seeding.
Background
The tickborne encephalitis (TBE) virus is a flavivirus persisting in a natural transmission cycle between small mammals and ticks. Humans can be infected, but they are ecologically deadend hosts [1]. TBE vectors in Central Europe are predominantly ticks of the genus Ixodes, especially Ixodes ricinus, the castor bean tick [2]. Since TBE can be a serious disease in humans [3], it is notifiable in almost all endemic TBE areas. Despite the availability of efficient vaccines [4, 5], TBE cases in Central Europe has risen sharply in recent decades [6]. In 2018, historical maximum values of 584 cases in Germany [7] and 377 cases in Switzerland [8] were registered. In Austria, 154 cases were the highest reported since 1994, although more than 80% of the population is vaccinated [9, 10]. Without vaccination probably more than 800 TBE cases per year would occur in Austria. Looking at the long time series of TBE cases, some of which date back to the 1950s [11], the question arises of how the temporal variations of these TBE cases can be explained. It can be taken into account that climate and environmental variables, averaged over large areas such as Central Europe, explain biological relationships much better than those with high local accuracy, as discussed in the fundamental papers on patterns and scales in ecology from Levin [12] and Hallett et al. [13]. Additionally, it can be taken into account that different mechanisms act on different time scales.
For example, longterm TBE trends, that have been observed over many decades, have been linked to factors such as demographic trends, changes in land use and associated wildlife density, or changes in human recreational behavior and related exposure [14]. Not least, climate change has been discussed as a possible driver [15, 16]. While in Sweden TBE incidence was significantly related to milder winters and higher spring and autumn temperatures [17], for the Baltic countries it was stated that climate change cannot explain the increase in TBE cases [18]. Here, it is assumed that climate change plays a only minor role in explaining the trend of Austrian TBE cases. Instead, the demographic development of the population is assumed to be the most probable cause for the rising TBE trend. This longterm trend in the Austrian TBE time series is superimposed by cyclical fluctuations. The duration period of these cyclical fluctuations was determined by Zeman [19] for 6 time series of TBE cases in Austria, the Czech Republic, the German federal states Bavaria and BadenWuerttemberg, Slovenia, and Switzerland. Calculating the power spectra from the detrended time series of Austrian TBE cases results in 2 dominant periods of the oscillations. The first has a period of 10 years (lowfrequency oscillations), the second has a period of 2–3 years (highfrequency oscillations) [19].
It is wellknown that the large atmospheric circulation variability is responsible for population and disease fluctuations [20]. Atmospheric circulation variability is also referred to as climate variability and is often described by socalled teleconnection indices. The best known of these climate variability or anomaly indices is the El Niño Southern Oscillation (ENSO), which occurs in areas around the tropical Pacific, especially in the southern hemisphere. ENSO triggered Malaria, Dengue, Rift Valley fever and other vectorborne disease outbreaks [21]. The ENSO impact on outbreaks reaches as far as the south of the USA, where a rodentborne hantavirus outbreak was associated with the 1997–1998 El Niño [21]. The most studied climate variability of the northern hemisphere is the North Atlantic Oscillation (NAO). It has been linked to a variety of disease outbreaks in the USA and Western Europe [22].
For example, Hubálek [23] studied 14 viral, bacterial and protozoan notifiable human diseases in the Czech Republic and their association with NAO indices, but no correlation was found for the tickborne diseases TBE and Lyme borreliosis. Palo [24] also found no correlation between NAO and the number of Swedish TBE cases. Another teleconnection index describing the largescale atmospheric circulation variability is the Scandinavian index (SCAND). It is less known than ENSO and NAO, and there is currently only one study that correlates human disease data, the UK asthma mortality, with SCAND fluctuations [25]. Since SCAND describes the largescale atmospheric circulation variability from Central Europe to Central Asia, it is hypothesized that it is suitable for describing the 10year oscillations in TBE cases in Austria.
The 2–3 year oscillations might be caused by the variations in beech fructification, which is responsible for the population dynamics of small mammals [26]. This also describes the oscillations in the population of I. ricinus whose larvae and nymphs feed mainly on yellownecked mice Apodemus flavicollis and bank voles Myodes glareolus [27] and thus contribute to the natural TBE virus transmission cycle. Brugger et al. [28] demonstrated that with the beech fructification index 2 years prior, the annual average temperature of the previous year and the past winter temperature, the I. ricinus nymphal density can be described with great accuracy. However, some peaks in TBE time series cannot be explained by tick density. An example are the extraordinary high TBE numbers in 2006, which were observed in some European countries (not in Austria). They were explained by recreational behavior of humans, i.e. more outdoor activities in the extremely warm year 2006 [29]. Because there are no longterm studies, a simple hypothesis is pursued, according to which more sunshine hours should lead to more outdoor activities and thus to a higher exposure of the human population.
Using the Austrian time series, we aimed at identifying the demographic and environmental factors associated with the trend and the oscillations by using generalized linear models (GLMs). Both the observed and the predicted TBE time series are to be subjected to a spectral analysis. The resulting power spectra should then indicate which predictors are responsible for the trend, the highfrequency and the lowfrequency oscillations, and with which explained variance they contribute to the TBE oscillations. The model development presented here differs fundamentally from the usual approaches, according to which significant predictors are selected from a large number of possible predictors by stepwise modeling [30], whereby their contribution to the frequency spectrum is not taken into account. So far, only two GLMs have been developed to predict TBE time series. The first is a GLM to predict the numbers of Swedish TBE cases by using December precipitation and red fox (Vulpes vulpes) or mink (Mustela vison) abundance as predictors [31]. The second GLM confirmed the high correlation between red fox density and TBE cases [24], which underpin the causal link between beech fructification, small mammals, and their predators, red foxes and mink.
Methods
Two ways of representing annual TBE time series are in use. On the one hand the absolute number of annual TBE cases is indicated, on the other hand the TBE incidence, i.e. the number of annual TBE cases per 100,000 inhabitants. Here, absolute numbers of annual TBE cases are used to allow the demographic parameters of the human population (Fig. 1) to be used as predictors. Thus, for example, the variance of the TBE time series explained by the population growth can be determined.
Demographics of the human population
Since the Austrian population has risen sharply in recent years, the demographic development must be taken into account. It is described by the birth rate, the mortality rate, and the net migration rate. Figure 1 shows the official demographic data [32]. According to this, the human population increased by more than 1.2 million in the period 1979–2018. This is mainly due to the net migration rate, i.e. the difference between annual immigration and emigration. Four major net migration (immigration) events occurred within the 40year period 1979–2018. The 2 most outstanding immigration events were caused by the Yugoslavian Civil War in 1991 and the Syrian Civil War in 2015. Net migration peaks were also observed 1981 after the suppression of the anticommunist social movement Solidarnosc in Poland and during 2001–2005 after the labor market has been opened further [33]. The difference between the birth rate and the mortality rate, the reproduction rate, is on average just 3,000 people per year. This is one order of magnitude lower than the mean net migration rate of about 30,000 people per year. Here, the total population N_{TOT} and the net migration N_{MIG} were used as predictors and listed in Table 1.
Human TBE time series
The official human TBE time series in Austria for the period 1979–2018 is analyzed. This 40year time series was documented by the Department of Virology, Medical University of Vienna, acting as the national reference laboratory for TBE virus infections. Only hospitalized patients with a serologically confirmed recent infection with TBE virus were counted as cases and published together with the vaccination coverage of the Austrian population [34]. In Austria, TBE is a notifiable disease and thus accuracy of the records is very high. Using the official vaccination coverage a hypothetical time series without TBE vaccination was estimated.
Here, N_{TBE} are the annual TBE cases documented by the national reference laboratory for TBE virus infections, VC is the official vaccination coverage within the interval [0, 1], and N is the hypothetical TBE cases without vaccination. Values of N_{TBE}, VC and N are listed in Table 1. In the following, only the hypothetical TBE cases N are used to investigate the natural trend and the oscillations in the Austrian TBE time series.
Climate teleconnection
To describe the decadal changes of the largescale climate in Central Europe several teleconnection indices are available. Here, the Scandinavian index (SCAND) developed by Barnstone and Livezey [35] was used, which the authors originally called the Eurasia1 pattern. With the help of SCAND an atmospheric circulation pattern, i.e. the spatial arrangement of northern hemisphere high and lowpressure systems, is characterized by a single index value. Like ENSO and NAO, the SCAND is therefore well suited for investigating correlations of largescale atmospheric circulation patterns with the cases of vectorborne diseases. A time series of the monthly SCAND is provided for the period 1950 to present on the Climate Prediction Center (CPC) website of the National Oceanic and Atmospheric Administration [36].
The SCAND describes an atmospheric circulation center over Scandinavia, with weaker centers of opposite sign over western Europe and eastern Russia/western Mongolia. Positive values are associated with belowaverage temperatures across western Europe and central Russia. It is also associated with aboveaverage precipitation across central and southern Europe and belowaverage precipitation across Scandinavia (Fig. 2). For the TBE endemic areas in Austria and southern Germany, this means that high SCAND values represent cooler and rainier periods. In turn, low SCAND values describe aboveaverage warm and dry periods.
To determine the optimal correlation between TBE and SCAND, socalled crosscorrelation maps (CCMs) were used. With CCMs optimal time lags and accumulation periods of predictors can be determined [37]. As known from vector biology, the best correlation between arthropod vectors or disease cases caused by pathogens they transmit and environmental temperature is obtained when temperature data were averaged over the period of the life cycle of the vector. For example, the life cycle of Culex pipiens, the vector of West Nile virus, is about 2–3 weeks during the mosquito activity period. With CCMs, 18 days were determined as the optimal averaging period [38]. If one wants to describe the dynamics of the Bluetonge virus vector Culicoides obsoletus by temperature, the somewhat longer averaging period of 37 days was estimated according the typical length of the life cycle of biting midges [39]. For TBE and its main vector I. ricnus the averaging period of the climate predictor SCAND should therefore be 2–6 years [40], 3–6 years [41], or 4–6 years [42]. In fact, the optimal averaging period determined by the application of CCMs is 4 years. Ideally, predictors should be normally distributed, which frequently can be achieved with a logtransformation. Therefore, the Scandinavian index SI used here is derived from the logtransformed monthly values of SCAND [36], which were averaged over 4 years. The SI values are listed in Table 1.
Beech fructification
The natural transmission cycle of the TBE virus depends on the availability of suitable hosts for the main vector I. ricinus. Preferred hosts of I. ricinus larvae are among others small rodents [43]. As no observations of small rodents are available for the long time series investigated in this study, the fructification index of the European beech (Fagus sylvatica) was applied for indicating the rodent density. Beechnuts are a basic food source for small rodents resulting in population peaks one year after mast seeding [44, 45]. Higher host densities cause higher densities of larvae of the TBE vectors I. ricinus. Two year after mast seeding, higher densities of I. ricinus nymphs are observed [28], which may be responsible for peaks in human TBE time series. Since the mast seeding is continentalscale synchronized [46], only a single time series, the beech fructification index published by Konnert et al. [47], is used here. This index is defined as the annual seed production and is divided in the following four classes: (0) absent, i.e. no fructification, (1) scarce, i.e. sporadic occurrence of fructification, but not noticeable at first sight, (2) common, i.e. clearly visible fructification, and (3) abundant, i.e. full fructification, also known as mast seeding. The values of the beech fructification indices 2 years prior F_{year−2} are listed in Table 1.
Annual sunshine duration
It is hypothesized that human recreational behavior affects the number of TBE cases. For example, higher outdoor activities increase exposure and they should therefore increase the TBE cases. Since there is no established predictor for human outdoor activities, the annual sunshine duration in hours is used as such here. With largescale considerations in focus, this should be representative for Central Europe. Of all meteorological services in Central Europe, only the German Weather Service offers such open data [48]. This is averaged over the entire region of Germany and should also be representative of the smaller neighboring countries such as Austria. Thus, in addition to the averaged logtransformed Scandinavian index SI, and the beech fructification index 2 years prior F_{year−2}, the annual sunshine duration SD is the third largescale predictor used for the analysis of TBE time series (Table 1).
Statistical modelling
All statistical analysis and modeling was done with the Language and Environment for Statistical Computing R [49]. GLMs were used to describe relationships. In the course of this, an overdispersion was observed since the dispersion parameter was generally greater than 1. This overdispersion was taken into account by using negative binomial models implemented with the R package mass [50].
To assess the necessary conditions for the application of GLMs, the predictors used (total population, net migration, Scandinavian index, beech fructification index, and annual sunshine duration) were tested for collinearity. This test is commonly used to select from a large number of predictor variables those that are most strongly correlated with the target variable, here the TBE cases. In addition, the so chosen predictor variables should be only weakly correlated with each other. Otherwise, their number can be further reduced. Here, however, a different approach is pursued: a small number of biologically well interpretable predictors are given. The check for collinearities (Fig. S1) therefore has only a control function, all correlations between the individual predictors are well below the threshold of R=0.7 [51]. With the R package psych [52] additional model diagnoses were created. These include examining the model errors for randomness (residual vs. fitted plot, scale location plot) and normal distribution (normal Q–Q plot), both of which are prerequisites for the applicability of GLMs [53]. Cook’s distance (residual vs. leverage plot) was used to test which TBE observations have the greatest influence on the regression. Outliers can be defined and eliminated if necessary [53], but this was not applied here.
While the statistical methods described above are intended to ensure the reliability of the selected model, it is particularly interesting how well the annual TBE fluctuations are described by the chosen predictor variables. Therefore, the models were verified by the rootmeansquare error (RMSE) and the explained variance (R^{2}). The advantage of these verification measures is that with the RMSE the error is specified in units of the target variable, i.e. the TBE cases, and R^{2} is well known.
If different GLMs are developed, the best model can be chosen with the help of the Akaike information criterion (AIC). The AIC estimates the quality of each model, relative to each of the other models. For better interpretability, however, the adjusted R^{2} (R\(^{2}_{{adj}}\)) is given here. In general, the use of additional predictors in GLMs leads to higher R^{2} values, even if they do not make a significant contribution to the model. With R\(^{2}_{{adj}}\) this is considered, which also determines model performance, relative to each of the other models [53].
A key objective of this study is to describe the causes of the trend as well as the lowfrequency and the highfrequency oscillations of the TBE cases. For this purpose, the power spectra of both the observed and the predicted TBE cases are calculated as described e.g. by [54]. The predictors for the GLMs should be selected so that the power spectra of the predictions match those of the observations as closely as possible.
Results
In 4 steps GLMs (negative binomial regression models) were developed, which demonstrate the influence of the selected predictors on the model performance as well as on the power spectrum of the predicted TBE time series. Thus, a final model was stepwise developed, which explains more than twothirds of the variance in the observations.
The first GLM uses only one predictive variable, the total population N_{TOT}. Figure 3(GLM1) shows the TBE cases without vaccination N (grey bars) with the predicted TBE cases N_{GLM1} (red line) representing a good approximation of the linear trend calculated from the observations N (black line). A rankorder correlation coefficient after Spearman of R=0.29 between the total population and the observed TBE cases was estimated (Fig. S1). The corresponding power spectrum for the observations shows 2 maxima. The first one is located at a period of 3 years (highfrequency oscillations), the second one is located at a period of 10 years (lowfrequency oscillations). The power spectrum of the model shows no maximum but only red noise, as expected for the trend.
The second GLM was extended to explain the lowfrequency oscillations in addition to the trend. Additional predictors used were the transformed Scandinavian index SI and the annual net migration N_{MIG}. Both contribute to the 10year TBE oscillation. The rankorder correlation between the TBE cases and the SI was R=0.52 (Fig. S1). In Austria, periods of high SI are related to relatively cool summers with aboveaverage precipitation (Fig. 2). Another highly significant contribution to explain the TBE cases is provided by the Austrian net migration rate, which can be considered as an socioeconomic predictor. The net migration N_{MIG} is negatively correlated with the numbers of the TBE cases (R=0.14). This suggests that new arrivals are less exposed to TBE virus infections, although they are responsible for the longterm population growth and thus also for the longterm increase in TBE cases. Since there is no study on this topic so far, it is hypothesized that the majority of immigrants from abroad initially settle in big cities where they are less exposed to TBE foci. The model therefore reduces the overestimated TBE cases during random net migration events. Figure 3(GLM2) shows the extended GLM verified with an error of RMSE=122 TBE cases and an explained variance of R^{2}=0.44 (R\(^{2}_{{adj}}\)=0.48). The power spectrum clearly shows that the 10year oscillation of the observed TBE time series is well associated with the 2 predictors SI and N_{MIG}.
The third GLM considers the beech fructification index 2 years prior F_{year−2} as an additional predictor for the highfrequency oscillations. The fact that the mast seeding 2 years prior has a high influence on the density of the main TBE vectors I. ricinus has already been shown by Brugger et al. [28]. Here it is demonstrated that this also applies to the TBE cases (R=0.38). Figure 3(GLM3) shows the model extended by the predictor beech fructification index F_{year−2}. The power spectrum clearly shows that the fructification index contributes to the explanation of the highfrequency oscillations, although the power is slightly too low compared to the spectrum of the observed TBE cases. The period of the lowfrequency oscillations, on the other hand, fits those of the observed TBE time series very well. The explained variance increases to R^{2}=0.64 (R\(^{2}_{{adj}}\)=0.66), resulting in a further reduction of the error of RMSE=98 cases. Considering the uncertainties in the observed TBE cases and the fact that time series of disease cases are generally difficult to explain, this result can be classified as very good. The parameter estimates and the significance levels of the predictors of this GLM are summarized in Table 2, where the discrete values of the fructification index are modeled as a factor. The factor F_{year−2}=0 is set as default and the factors F_{year−2}=1, 2 or 3 are considered by different parameter estimates. Thus, the GLM requires only the four predictors N_{TOT},N_{MIG}, SI and F_{year−2}, all of which contribute to the model very significantly (p <0.001).
The fourth GLM was extended by a predictor for the outdoor activity of humans. So far, no influence of human recreational behavior on the annual TBE time series has been considered, which should lead to a further improvement of the model. A climatic parameter that should have a plausible influence on an increased outdoor activity of humans is annual sunshine duration SD in hours. It is directly (without a lag time) correlated with the TBE cases (R=0.27). Thus, the correlation between SD and N is similar to the correlation between F_{year−2} und N. Since there is no appreciable collinearity between SD and F_{year−2} (Fig. S1), the consideration of SD results in an increased model performance. Figure 3(GLM4) shows the results of the GLM with the additional predictor SD. It explains the variations of the TBE cases even better, namely with R^{2}=0.70 (R\(^{2}_{{adj}}\)=0.70), which leads to a further reduction of the error of RMSE=89 cases. The parameter estimates and the significance levels of the predictors of the final model are summarized in Table 3. Again, all predictors contribute significantly to model performance, with most pvalues being very significant. Of course, the relative contribution of the fructification index decreases at the expense of sunshine duration, as both predictor variables are responsible for highfrequency oscillations. Statistical features for the final GLM as described in “Statistical modelling” section are provided in Fig. S2, the AIC values for the stepwise developed models in Table S1.
Discussion
The climate of the TBE endemic areas in Austria and neighboring countries is still characterized by the boreal coniferous climate of the northern hemisphere in the 1960s. According to the wellknown KöppenGeiger climate classification it is known as Dfb climate, a boreal climate with rain at all seasons and warm summers [55]. Until today, this boreal coniferous climate has almost completely retreated from the Alps and has been replaced by a warm temperate Cfb climate [56]. It is called, according to the prevailing tree species in natural forests, beech climate. In order to take this into account, spruce monocultures threatened by climate change are being gradually replaced by deciduous or mixed forests across the region of the greater Alps. These more speciesrich forests also provide better living conditions for the most important TBE virus vector I. ricinus, resulting in higher I. ricinus densities [57, 58]. Beech fructification is a predictor of the intensity of the natural TBE virus transmission cycle between small mammals and ticks, with a high beech fructification index increasing the population density of small mammals and of I. ricinus larvae one year thereafter. One more year later, significantly higher densities of questing I. ricinus nymphs are responsible for the more frequent transmission of the TBE virus to humans [28]. However, effects on TBE cases due to these changes in forestry were only visible in the last decade, where higher frequencies of years with full fructification of beech (mast seeding) are responsible for several peaks in the TBE time series. A total of 3 of the 5 years with full fructification occurred after 2010 (Table 1). The effects of climate change are therefore not visible through a direct correlation of the TBE cases with rising temperatures, but indirectly via the increased frequency of mast seeding (F_{year−2}). In combination with the rapidly increasing human population (N _{TOT}, N _{MIG}) and a slight decline in vaccination coverage (VC), this explains the major effects of rising numbers of Austrian TBE cases observed after 1995 under real conditions with vaccination (see N_{TBE} in Table 1). Additionally, 10year oscillations are associated with the largescale distribution of atmospheric high and low pressure systems (SI) resulting in the third model version GLM3, which explains 64% of the variation of Austrian TBE cases. In Austria, periods of high SI are related to relatively cool summers with aboveaverage precipitation (Fig. 2). This may influence the density of the TBE vector I. ricinus, which does not like hot summers and extreme drought. Remarkably, oscillations with periods of 10 and 3–4 years were also observed for the TBE vector Ixodes persulcatus in Russia. Years with high tick densities follow frequently those with the peak population density of small mammals [41]. However, no direct correlation between human TBE cases and the TBE vectors I. ricinus and I. persulcatus has yet been published. Similar ecological connections are also known from the oak forests of eastern North America [44].
With the additional predictor sunshine duration SD in GLM4 the explained variance continues to increase. The hypothesis was that the exposure of the population increases with an increasing number of hours of sunshine. This hypothesis seems confirmed, as the explained variance in GLM4 increased to 70%. But that must not hide the fact that further studies on the contribution of human behavior to the cases of TBE are needed. It should be noted that GLM4 cannot be used for TBE forecasts because neither social behavior itself (e.g. during the unexpected SARSCoV2 pandemic 2020) nor SD can be estimated for the next 1–2 years. There is also a fraction of unexplained variance of 30%, which needs further research. In particular, rare extreme events are difficult to detect by statistics, because low case numbers result in low significance. For example, Dautel et al. [59] have shown for Germany that extreme low temperatures in January and February 2012, in combination with the lack of a protective snow cover, led to decreasing numbers of I. ricinus nymphs as well as very low numbers of human TBE cases in the same year (also recognizable in the Austrian TBE time series). The inclusion of this and similar field studies could help to improve future predictions.
Another aspect that has not been considered so far concerns the gender and age distribution of TBE cases within the population. The age distribution of the TBE cases shows a maximum at 55 years [3], with generally more men being infected with the TBE virus [60]. It has not yet been investigated whether the increasing number of TBE cases are related to the aging society.
It should also be noted that alternative predictors for the explanation of TBE cases are mentioned in the literature. In Slovenia, for example, a correlation was found between TBE cases and roe deer density 3 years ago [61]. A high roe deer (Capreolus capreolus) density is interpreted as a high host density for the TBE vector I. ricinus and consequently responsible for high TBE cases. For Austria, the results of Knap and Avs̆ic̆Z̆upanc [61] could be confirmed, but the explained variance of the GLM4 decreased with the use of roe deer density instead of SI from 70% to 58%. This is because the correlation between the roe deer density and the TBE cases is largely due to the concurrent trends. The estimation of unknown wildlife densities from hunting index generally leads to some uncertainties. In addition, Austrian hunting data [62] from the statistical database STATcube provided by Statistics Austria [63] are not available in near realtime, as the hunting year differs from the calendar year. It covers the period from April 1 to March 31 of the following year. Wildlife data are therefore generally less wellsuited as predictors than climate data, especially with regard to a possible forecast of the next year’s TBE cases.
Conclusions
The TBE models presented here confirm the work of Hallett et al. [13] that largescale indices can predict ecological processes very well, probably better than local weather and climate parameters. As the beech fructification indices of the years 2017 and 2018 are responsible for the TBE cases in 2019 and 2020, GLM3 can also be applied to forecast the TBE cases of the next 2 years. This is possible because for the forecast mainly the highfrequency oscillations caused by F_{year−2} are interesting. The predictors relevant to the trend and the lowfrequency oscillations, on the other hand, can be extrapolated by simple methods such as persistence or linear interpolation. GLM3 is also applicable to the neighboring countries such as Germany or Switzerland using the same predictors, since the Scandinavian index is representative for all of Central Europe and also the beech fructification is largescale synchronized. This has be examined in a followup study that was carried out during the review process of the paper presented here [64]. The verification with independent TBE cases from 2019 has demonstrated the good performance of the forecasts.
Finally, it should be noted that the findings presented here can subsequently be used to create process models of the type susceptibleinfectedrecovered (SIR). These models represent the highest stage of development in epidemiological modelling as, unlike statistical models, they map the dynamics of population health based on the underlying processes of disease transmission. A first SIR model on the dynamics of the Austrian human TBE cases was presented by Rubel [65].
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Abbreviations
 TBE:

Tickborne encephalitis
 GLM:

Generalized linear model
 ENSO:

El Niño Southern Oscillation
 NAO:

North Atlantic Oscillation
 SCAND:

Scandinavian index
 SI:

Transformed Scandinavian index
 VC:

Vaccination coverage
 SD:

Sunshine duration
 CCM:

Cross correlation map
 RMSE:

Rootmeansquare error
 AIC:

Akaike information criterion
 SIR model:

Susceptibleinfectedrecovered model
References
 1
Kahl O, Pogodina VV, Poponnikova T, Süss J, Zlobin VI. A short history of TBE In: Dobler G, Erber W, Bröker M, Schmitt HJ, editors. The TBE Book, 2. edn. Chap. 1. Singapore: Global Health Press: 2019. p. 11–18.
 2
Lindquist L, Vapalahti O. Tickborne encephalitis. Lancet. 2008; 371:1861–71.
 3
Kaiser R. The clinical and epidemiological profile of tickborne encephalitis in southern Germany 1994–98. Brain. 1999; 122:2067–78.
 4
Heinz FX, Holzmann H, Essl A, Kundi M. Field effectiveness of vaccination against tickborne encephalitis. Vaccine. 2007; 25:7559–67.
 5
Heinz FX, Stiasny K, Holzmann H, GrgicVitek M, Kriz B, Essl A, Kundi M. Vaccination and tickborne encephalitis, Central Europe. Emerg Infect Dis. 2013; 19:69–76.
 6
Süss J. Tickborne encephalitis 2010: Epidemiology, risk areas, and virus strains in Europe and Asia  An overview. Ticks Tick Borne Dis. 2011; 2:2–15.
 7
Robert KochInstitut. SurvStat@RKI 2.0. 2019. https://survstat.rki.de (visited on February 7, 2019).
 8
Bundesamt für Gesundheit. FrühsommerMeningoenzephalitis (FSME): Ausweitung der Risikogebiete (in German). BAGBulletin. 2019; 6:12–4.
 9
Kunz C. TBE vaccination and the Austrian experience. Vaccine. 2003; 21(S1):50–5.
 10
Kunze U, Böhm G. FrühsommerMeningoEnzephalitis (FSME) und FSMESchutzimpfung in Österreich: Update 2014. Wien Med Wochenschr. 2015; 165:290–5.
 11
Dobler G, Erber W, Bröker M, Schmitt HJ. The TBE Book, 2nd edn.Singapore: Global Health Press; 2019.
 12
Levin SA. The problem of pattern and scale in ecology. Ecology. 1992; 73:1943–67.
 13
Hallett TB, Coulson T, Pilkington JG, CluttonBrock TH, Pemberton JM, Grenfell BT. Why largescale climate indices seem to predict ecological processes better than local weather. Nature. 2004; 430:71–5.
 14
Zeman P, Benes C. Spatial distribution of a population at risk: An important factor for understanding the recent rise in tickborne diseases (Lyme borreliosis and tickborne encephalitis in the Czech Republic). Ticks Tick Borne Dis. 2013; 4:522–30.
 15
Daniel M, Danielová V, Fialová A, Malý M, Kŕiz̆ B, Nuttall PA. Increased relative risk of tickborne encephalitis in warmer weather. Front Cell Infect Microbiol. 2007; 8:90.
 16
Gray JS, Dautel H, EstradaPeña A, Kahl O, Lindgren E. Effects of climate change on ticks and tickborne diseases in Europe. Interdiscip Perspect Inf Dis. 2009; id=593232.
 17
Lindgren E, Gustafson R. Tickborne encephalitis in Sweden and climate change. The Lancet. 2001; 358:16–8.
 18
Sumilo D, Asokliene L, Bormane A, Vasilenko V, Golovljova I, Randolph SE. Climate change cannot explain the upsurge of tickborne encephalitis in the Baltics. PLoS ONE. 2007; 2:500.
 19
Zeman P. Cyclic patterns in the central European tickborne encephalitis incidence series. Epidemiol Infect. 2017; 145:358–67.
 20
Stenseth NC, Mysterud A, Ottersen G, Hurrell JW, Chan KS, Lima M. Ecological effects of climate fluctuations. Science. 2002; 297:92–6.
 21
Kovats RS, Bouma MJ, Hajat S, Worrall E, Haines A. El Niño and health. The Lancet. 2007; 362:1481–9.
 22
Morand S, Owers KA, WaretSzkuta A, McIntyre KM, Baylis M. Climate variability and outbreaks of infectious diseases in Europe. Sci Rep. 2013; 3:1774.
 23
Hubálek Z. North Atlantic weather oscillation and human infectious diseases in the Czech Republic, 1951–2003. Europ J Epidemiol. 2005; 20:263–70.
 24
Palo RT. Tickborne encephalitis transmission risk: Its dependence on host population dynamics and climate effects. Vectorborne Zoon Dis. 2014; 14:346–52.
 25
Majeed H, Moore GWK. Influence of the Scandinavian climate pattern on the UK asthma mortality: a time series and geospatial study. BMJ Open. 2018; 8:020822.
 26
Clement J, Maes P, van Ypersele de Strihou C, van der Groen G, Barrios JM, Verstraeten WW, van Ranst M. Beechnuts and outbreaks of Nephropathia epidemica (NE): of mast, mice and men. Nephrol Dial Transplant. 2010; 25:1740–6.
 27
Galfsky D, Król N, Pfeffer M, Obiegala A. Longterm trends of tickborne pathogens in regard to small mammal and tick populations from Saxony, Germany. Parasit Vectors. 2019; 12(1). https://doi.org/10.1186/s1307101933822.
 28
Brugger K, Walter M, ChitimiaDobler L, Dobler G, Rubel F. Forecasting next season’s Ixodes ricinus nymphal density: the example of southern Germany 2018. Exp Appl Acarol. 2018; 75:281–8.
 29
Randolph SE, Asokliene L, AvsicZupanc T, Bormane A, Burri C, Gern L, Golovljova I, Hubalek Z, Knap N, Kondrusik M, Kupca A, Pejcoch M, Vasilenko V, Z̆ygutiene M. Variable spikes in tickborne encephalitis incidence in 2006 independent of variable tick abundance but related to weather. Parasit Vectors. 2008; 1:44.
 30
Whittingham MJ, Stephens PA, Bradbury RB, Freckleton RP. Why do we still use stepwise modelling in ecology and behaviour?J Anim Ecol. 2006; 75:1182–9.
 31
Haemig PD, De Luna SS, Grafström A, Lithner S, Lundkvist A, Waldenström J, Kindberg J, Stedt J, Olsén B. Forecasting risk of tickborne encephalitis (TBE): Using data from wildlife and climate to predict next year’s number of human victims. Scand J Infect Dis. 2011; 43:366–72.
 32
Statistics Austria. Bevölkerungsstand und veränderung (in German). Vienna; 2018. https://www.statistik.at (visited on November 23, 2018).
 33
Bischof G, Rupnow D. Migration in Austria. Innsbruck: Innsbruck University Press; 2017.
 34
Heinz FX, Stiasny K, Holzmann H, Kundi M, Sixl W, Wenk W, Kainz W, Essl A, Kunz C. Emergence of tickborne encephalitis in new endemic areas in Austria: 42 years of surveillance. Euro Surveill. 2015; 20(13):21077.
 35
Barnstone AG, Livezey RE. Classification, seasonality and persistence of lowfrequency atmospheric circulation patterns. Mon Wea Rev. 1987; 115:1083–126.
 36
NOAA. Northern Hemisphere Teleconnection Patterns, Scandinavia (SCAND). Maryland, USA: Climate Prediction Center of the National Oceanic and Atmospheric Administration (NOAA); 2019. http://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml.
 37
Brugger K, Walter M, ChitimiaDobler L, Dobler G, Rubel F. Seasonal cycles of the TBE and Lyme borreliosis vector Ixodes ricinus modelled with timelagged and intervalaveraged predictors. Exp Appl Acarol. 2017; 73:439–50.
 38
Lebl K, Brugger K, Rubel F. Predicting Culex pipiens/restuans population dynamics by interval lagged weather data. Paraisit Vectors. 2013; 6:129.
 39
Brugger K, Rubel F. Bluetongue disease risk assessment based on observed and projected Culicoides obsoletus spp. vector densities. PLoS ONE. 2013; 8(4):60330.
 40
Gray JS. The development and seasonal activity of the tick Ixodes ricinus: a vector of Lyme borreliosis. Rev Med Vet Entomol. 1991; 79:323–33.
 41
Balashov YS. Demography and population models of ticks of the genus Ixodes with longterm life cycles. Entomol Rev. 2012; 92:323–33.
 42
Kahl O, Petney TN. Biologie und Ökologie des wichtigsten FSMEVirusÜberträgers in Mitteleuropa, der Zecke Ixodes ricinus (in German) In: Rubel F, SchiffnerRohe J, editors. FSME in Deutschland: Stand der Wissenschaft. BadenBaden: Deutscher Wissenschaftsverlag: 2019. p. 23–38.
 43
Gray JS, Kahl O, Lane RS, Levind ML, Tsaoe JI. Diapause in ticks of the medically important Ixodes ricinus species complex. Ticks TickBorne Dis. 2016; 7:992–1003. https://doi.org/10.1016/j.ttbdis.2016.05.006.
 44
Ostfeld RS, Jones CG, Wolff JO. Of mice and mast: Ecological connections in eastern deciduous forests. BioScience. 1996; 46:323–30. https://doi.org/10.2307/1312946.
 45
Clement J, Vercauteren J, Verstraeten WW, Ducoffre G, Barrios JM, Vandamme AM, Maes P, Van Ranst M. Relating increasing hantavirus incidences to the changing climate: the mast connection. Int J Health Geogr. 2009; 8:1.
 46
Nussbaumer A, Waldner P, Etzold S, Gessler A, Benham S, Thomsen IM, Jørgensen BB, Timmermann V, Verstraeten A, Sioen G, Rautio P, Ukonmaanaho L, Skudnik M, Apuhtin V, Braun S, Wauer A. Patterns of mast fruiting of common beech, sessile and common oak, Norway spruce and Scots pine in Central and Northern Europe. Forest Ecol Manage. 2016; 363:237–51.
 47
Konnert M, Schneck D, Zollner A. Blühen und Fruktifizieren unserer Waldbäume in den letzten 60 Jahren. LWF Wissen. 2016; 74:37–45.
 48
Deutscher Wetterdienst. Open Data Provided by the Climate Data Center of the German Weather Service. Offenbach am Main. 2019. https://opendata.dwd.de/climate_environment/CDC. Last access 11 June 2019.
 49
R Development Core Team. R: A Language and Environment for Statistical Computing, Version 3.3.2. Vienna, Austria: R Foundation for Statistical Computing; 2016. ISBN 3900051070, http://www.Rproject.org.
 50
Ripley B, Venables B, Bates DM, Hornik K, Gebhardt A, Firth D. Mass: Support Functions and Datasets for Venables and Ripley’s MASS (Modern Applied Statistics with S, 4th Edition, 2002). 2019. R package version 7.351.4.
 51
Dormann CF, Elith J, Bacher S, Buchmann C, Carl G, Carré G, Marquéz JRG, Gruber B, Lafourcade B, Leitäo PJ, Münkemüller T, McClean C, Osborne PE, Reineking B, Schröder B, Skidmore AK, Zurell D, Lautenbach S. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography. 2013; 36:27–46.
 52
Revelle W. Psych: Procedures for Psychological, Psychometric, and Personality Research. Evanston, Illinois: Northwestern University; 2018. https://CRAN.Rproject.org/package=psych. R package version 1.8.12.
 53
Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM. Mixed Effects Models and Extensions in Ecology with R. New York: Springer; 2009.
 54
Stull RB. An Introduction to Boundary Layer Meteorology. Dordrecht: Kluwer Academic Publishers; 1991.
 55
Kottek M, Grieser J, Beck C, Rudolf B, Rubel F. World map of the KöppenGeiger climate classification updated. Meteorol Z. 2006; 15:259–63.
 56
Rubel F, Brugger K, Haslinger K, Auer I. The climate of the European Alps: Shift of very high resolution KöppenGeiger climate zones 1800–2100. Meteorol Z. 2017; 26:115–25.
 57
Boehnke D, Brugger K, Pfäffle M, Sebastian P, Norra S, Petney T, Oehme R, Littwin N, Lebl K, Raith J, Walter M, Gebhardt R, Rubel F. Estimating Ixodes ricinus densities on the landscape scale. Int J Health Geogr. 2015; 14:23.
 58
Brugger K, Boehnke D, Petney T, Dobler G, Pfeffer M, Silaghi C, Schaub GA, Pinior B, Dautel H, Kahl O, Pfister K, Süss J, Rubel F. A density map of the tickborne encephalitis and Lyme borreliosis vector Ixodes ricinus (Acari: Ixodidae) for Germany. J Med Entomol. 2016; 53:1292–302.
 59
Dautel H, Kämmer D, Kahl O. How an extreme weather spell in winter can influence vector tick abundance and tickborne disease incidence In: Braks MAH, van Wieren SE, Takken W, Sprong H, editors. Ecology and Prevention of Lyme Borreliosis. Ecology and Control of Vectorborne Diseases, vol. 4. Wageningen: Wageningen Academic Publishers: 2016. p. 335–49.
 60
Kiffner C, Zucchini W, Schomaker P, Vor T, Hagedorn P, Niedrig M, Rühe F. Determinants of tickborne encephalitis in counties of southern Germany, 20012008. Int J Health Geogr. 2010; 9:42.
 61
Knap N, Avs̆ic̆Z̆upanc T. Correlation of TBE incidence with red deer and roe deer abundance in Slovenia. PLoS ONE. 2013; 8(6):66380.
 62
Reimoser S, Reimoser F. Habitat quality & hunting bag: culling densities of different wildlife species in Austria since 1955. Part 1: Roe deer (Capreolus capreolus) (in German). Weidwerk. 2005; 6/2005:14–5.
 63
Statistics Austria. Statistical Database STATcube. 2020. https://www.statistik.at/web_en/publications_services/statcube/index.html. Last access 04 Feb 2020.
 64
Rubel F, Brugger K. Tickborne encephalitis incidence forecasts for Austria, Germany, and Switzerland. Ticks Tick Borne Dis. 2020; 11:101437. https://doi.org/10.1016/j.ttbdis.2020.101437.
 65
Rubel F. Erklärende Modelle zur Dynamik der FSMEErkrankungen In: Rubel F, SchiffnerRohe J, editors. FSME in Deutschland: Stand der Wissenschaft. BadenBaden: Deutscher WissenschaftsVerlag: 2019. p. 243–60.
Acknowledgements
The authors are very grateful to Olaf Kahl for critical reading the manuscript and his valuable hints.
Funding
No funding was obtained for this study.
Author information
Affiliations
Contributions
FR, MW, JRV, and KB contributed to the study conception and design. Material preparation, data collection and analysis were performed by FR and KB. The first draft of the manuscript was written by FR and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
FR received speaker honoraria from Pfizer Deutschland GmbH. KB received unrestricted research grants from Pfizer Deutschland GmbH. All authors received travel grants and author honoraria from Pfizer Deutschland GmbH.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Additional file 1
Additional file S1. Frequency distributions (red bars) of the Austrian TBE incidence without vaccination N and the predictors used in the generalized linear models (GLMs): total human population N_{TOT}, net migration N_{MIG}, transformed scandinavian index SI, beech fructification index 2 years prior F_{year−2} and, annual sunshine duration in hours SD. the following rankorder correlations with N have been determined: r=0.29 (N_{TOT}), r=0.14 (N_{MIG}), r=0.52 (SI), r=0.38 (F_{year−2}), and r=0.27 (SD). maximum collinearity of r=0.59 (N_{TOT} vs. N_{MIG}). Additional file S2. Statistical features of GLM4. Additional file S3. Akaike information criterion (AIC) and explained variance (R\(^{2}_{{adj}}\)) values for stepwise developed models GLM1–GLM4.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Rubel, F., Walter, M., Vogelgesang, J.R. et al. Tickborne encephalitis (TBE) cases are not random: explaining trend, low and highfrequency oscillations based on the Austrian TBE time series. BMC Infect Dis 20, 448 (2020). https://doi.org/10.1186/s12879020051567
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12879020051567
Keywords
 Vectorborne disease
 Climate change
 Scandinavian index
 Beech fructification
 Mast seeding
 Prediction
 Ixodes ricinus