INTRODUCTION
Crimean-Congo haemorrhagic fever (CCHF) has probably been present in Turkey since the late 1990s, but no confirmed cases occurred before 2002 [1]. During spring and summer of 2002, the first patients with clinical and laboratory findings compatible with CCHF were reported in the city of Tokat, located in the Middle Anatolia–Black Sea region of Turkey. In the spring and summer of 2003, patients with similar clinical and laboratory findings were identified in Tokat and neighbouring cities. CCHF has been listed as a notifiable disease in Turkey since December 2003 [Reference Yilmaz2]. Since that time, the Turkish Ministry of Health has tracked notified cases, and has maintained an adequately georeferenced dataset for places of infection.
The tick Hyalomma marginatum Koch, 1844, is widely distributed throughout the Mediterranean Basin [Reference Apanaskevich and Horak3] and has been incriminated as the main vector of CCHF for humans in the region, with a complex chain of reservoirs, hosts of immature stages, and social features that lead to close contact between ticks and humans [Reference Estrada-Peña4–Reference Ergönul7]. It is accepted that a series of social changes were responsible for previous epidemics in ex-Soviet regions and that reservoir hosts and their abundance played a pivotal role in the epidemiology of the disease. Further, it has been suggested that fragmentation of vegetation may dramatically increase appropriate habitat for ticks and their hosts, as well as resulting in a critical rise in the rate of contact between humans and vectors in some tick-borne systems [Reference Allan, Keesing and Ostfeld8, Reference Jackson, Hilborn and Thomas9]. This is very relevant to studies on CCHF, as previous reports, reviewed in [Reference Hoogstraal5], studied the effects of land abandonment and subsequent re-colonization as key factors in CCHF epidemiology. According to the cited authors, very large populations of H. marginatum hosts colonized farming land that had been abandoned for a number of years. The later entry of humans into these abandoned territories and the resumption of farm labour provided a good environment for increased transmission rates because of amplification of the virus in reservoir hosts followed by contact with humans.
However, the disease did not spread into all areas where the main tick vector, H. marginatum, is recorded in Turkey. A recent study hypothesized a link between climate change and increased transmission rates [Reference Purnak, Selvi and Altundag10], whereas another study stated that the trend in CCHF reporting depends only upon social changes [Reference Ergönul7]. However, no systematic studies have been performed on the roles played by the many abiotic factors involved in the biology of the tick vector during the years in which the disease was reported to be spreading across Turkey. In this paper we specifically aimed to address the question of whether the disease occupies a well-defined climate niche. We explored whether areas with autochthonous cases shared characteristic features, coherently observed in both time and space, and distinct from regions where the disease has not been reported. Such a niche would be reflected in monthly and seasonal composites of atmospheric and ground-derived climate, long-term climate trends, and landscape fragmentation. To address these questions, we used the dataset on case reports in Turkey between 2003 and 2008, together with remotely sensed information on temperature and vegetation, as well as assimilation data calculated at and below ground level.
MATERIALS AND METHODS
Sources of cases
Data on CCHF cases were compiled by the Turkish Ministry of Health, and the documentation includes details on the place of residence of each patient and the geographical location where each patient noticed the tick bite. All cases were further adequately georeferenced by assigning an unambiguous pair of coordinates to each case by reference to local gazetteers. For this study we used annual data between 2003 and 2008. All cases were confirmed by ELISA and/or PCR [Reference Yilmaz2]. The cases supported only by clinical findings were not included in this study.
Sources of climate and vegetation data
All the variables used in the study had an impact on tick biology, regulating both their abundance and involvement in the timing of seasonal activity. Temperature has an obvious effect on the regulation of the life-cycle of arthropods, and the Normalized Derived Vegetation Index (NDVI), a measure of photosynthetic activity [Reference Hoogstraal5], was used as a proxy for vegetation stress and water balance. Furthermore, we used both an air-derived (atmospheric) dataset and a ground-derived one. As the tick vector spends non-parasitic stages in the ground, overwintering at a variable depth below the surface, we used climate features at and below the surface. Our aim was to assemble the longest and most homogeneous dataset available, sometimes compiled from different sources, to produce reliable statistical results. For atmospheric mean temperatures and medium resolution NDVI information, we used the monthly dataset provided by NEO [NASA Earth Observations (http://neo.sci.gsfc.nasa.gov)]. These are layers obtained from the NOAA-AVHRR series of satellites, providing a homogeneous and calibrated set of data on environmental features, already corrected for atmospheric disturbance and orbital variations. The resolution of the dataset is 0·1°. Both temperature and NDVI values are calculated by NEO as the monthly mean of daily maximum values. Ground datasets were obtained from the GLDAS-DISC series of assimilation models available at http://agdisc.gsfc.nasa.gov/data/browse/GLDAS/GLDAS_NOAH025_M. These databases included mean monthly temperatures (in °C) at ground level and a depth of 4 cm, mean evaporation over the ground surface (in kg/m2.s), near-surface mean specific humidity (kg/kg), as well as mean moisture content 4 cm below the ground surface (in kg/m2). All data in GLDAS are at a nominal resolution of 0·25°. Both sets of variables, from NEO and GLDAS-DISC, span the interval 2000–2008.
NDVI mean monthly data at a resolution of 1 km were obtained from the SPOT satellite repository (http://free.vgt.vito.be). These data were exclusively used to produce an evaluation of landscape fragmentation, unattainable at the resolution of the NEO dataset. Monthly raster layers were obtained between 1999 and 2008.
Spatial and statistical procedures
Because of the different resolution of the datasets used in our study, we built a grid of hexagonal cells 0·15° in radius overlying the entire country. The procedure involved the choice of an adequate resolution and the preparation of a grid covering the whole territory to be studied. Then, the values from the original raster layers containing the environmental variables were transferred to the cells, each cell containing the mean value of the pixels covering that cell. Our geographical work unit was the hexagonal cell (hereinafter termed ‘cell’). Each cell was also marked as positive or negative according to CCHF case reporting for each year during the period 2003–2008. The number of positive and negative cells for each year was 71/1800 (2003), 123/1748 (2004), 155/1716 (2005), 204/1667 (2006), 284/1587, and 405/1466 (2008).
We first searched for significant differences in monthly variables between positive and negative cells, for both the year of case reporting (hereinafter termed the ‘current year’) and the previous year. Simple ANOVA tests were performed to compare monthly mean temperatures and NDVI values (dependent variables) with the case status of each cell (positive or negative), as the independent variable. Monthly variables recorded below the ground (humidity, evaporation, moisture, temperature) were processed in the same manner. We checked for coherence in the analysis, as many variables showed statistical variations between different years. Therefore, in order to assign a statistically significant effect to any variable, positive and negative cells must show consistently significant differences in every year of the study period. To restrict geographical extent and to avoid extreme climatic conditions that would bias the analysis, for negative sites we used cells where no cases had been reported during 2003–2008, but where the tick vector was known to be present. The distribution of the vector in Turkey was obtained from field surveys (compiled by Z.V.) and from a previous model of tick distribution [Reference Estrada-Peña11].
After studies on monthly variables, we analysed seasonal components of climate and vegetation data, as some of these variables had previously been identified as valuable markers for epidemiology of the disease [Reference Estrada-Peña4]. We used daily values of air temperature obtained from Fourier analysis of the NEO dataset to compute the length of growth and summer seasons [Reference Grimenes and Nissen12]. Total degree days accumulated in winter (defined as 365 days minus the length of the growing season [Reference Grimenes and Nissen12]) were also obtained. Similar analyses were performed with NDVI values in a search for the span of the NDVI anomaly (NDVIa). We defined the beginning of the NDVIa as the week of the year in which the NDVI signal became higher than the mean of the previous year [Reference Grimenes and Nissen12]. Thus, we computed the total number of days with an anomaly in the NDVI signal for each year. ANOVA analysis was used to compare positive and negative cells (within the range of expected tick presence) in a search for differences between the growth and summer seasons, accumulated degree days in winter, and NDVIa, for both the current and previous year.
The trends in climate and medium-resolution NDVI values were studied by means of Fourier decomposition, which seeks a periodical signal in a time-series of data. Such decomposition can give slopes of changing specific Fourier variables (e.g. amplitudes or phases of annual, biannual, triannual cycles) over a number of years. Of particular importance was calculation of a slope which showed the trend of a studied time-series (where a zero slope indicated no change). Significant differences in climate and vegetation trends between positive and negative cells were analysed by ANOVA, using the presence/absence of cases in each cell as the independent variable and the trend in temperature or NDVI as dependent variables.
Analysis of vegetation and landscape patterns
According to network theories of disease transmission and spread, it has been reported that highly fragmented habitat leads to overspread of disease [Reference Schimit and Monteiro13]. Our aim was to evaluate changes in habitat fragmentation in Turkey and to associate different patterns of fragmentation with disease reporting (presence/absence and case incidence). We used monthly NDVI data at a resolution of 1 km, captured by the SPOT satellite (S10 sensor) between 2000 and 2008. We performed a supervised classification of vegetation, separately for each year, using the baseline as prepared by the consortium GLC2000 (available at http://gem.jrc.ec.europa.eu/index.php), which was achieved with images captured by the same sensor between 1999 and 2000. To avoid noise in the classification, we removed every classified patch smaller than 9 pixels (i.e. 9 km2) which was then replaced with the surrounding larger patch category of vegetation.
The landscape reveals a new set of constraints for ticks (e.g. regarding dispersal), compared to the local constraints of survival and development derived from climate. The landscape structure integrates the spatial arrangement of habitat patches and has a deep influence on tick populations at patch scale. We thus needed to obtain a value relating the habitat fragmentation features with the abundance of ticks, since previous findings on the topic reported the correlation between CCHF case incidence and high fragmentation rates of areas suitable for the tick vector [Reference Estrada-Peña11]. We computed an index indicating the consequences of fragmentation on the habitat perception by tick-carrying hosts and therefore on the tick vector; this is the so-called recruitment index. The development of the recruitment index has deep roots in graph theory [Reference Urban and Keitt14] that relates the size of patches, their abundance and connectivity, with the probability of persistence of a local population of animals or plants. In this context, it is important to emphasize that a previous field study empirically demonstrated the correlation between recruitment values (computed from habitat fragmentation) and abundance of the tick Ixodes ricinus [Reference Estrada-Peña15, Reference Estrada-Peña16] in a wide area. Therefore, small and well-connected patches had high recruitment values (and therefore high local tick densities), whereas large, slightly fragmented patches, separate from other patches, tended to show low recruitment because of low connectivity (and thus low tick densities). Details of the derivation of the recruitment index, the statistical reasoning behind such modelling, and comparison of recruitment values with empirically obtained tick abundance data, have been provided previously [Reference Estrada-Peña17]. We expected an obvious difference in the structure of the I. ricinus meta-population for which the recruitment index was originally obtained, and that of Hyalomma ticks involved in the current study. However previous results agreed with general studies on patch occupancy and population dispersal (summarized in [Reference Wilson18]) accounting for a generalized application of the recruitment index to different tick species.
Cells were marked according to case reporting, separately for each year. Whereas case reporting spanned 2003–2008, recruitment values were included for the period 2000–2008 for further comparison of long-term trends. It is inadequate to ascribe simple case incidence to cells, because it is difficult to know the number of people at risk in any particular cell. It is well known that the rural population is at higher risk than urban dwellers, because rural residents are likely to have more contact with infective ticks. Therefore, a calculation of case incidence using the total population in each cell would bias the result. An adequate census at the district level exists for Turkey, compiled by the Turkish Government and updated to 2005, separating the rural and urban populations (http://www.turkstat.gov.tr/jsp/duyuru/upload/vt_en/vt.htm). Because of the relatively small and homogeneous size of district administrative divisions, we decided to calculate annual case incidence using rural population numbers and the number of cases reported for each district. Recruitment values were calculated for each district and compared with case incidence. The association between these values was estimated separately for districts with a case incidence higher or lower than 20 cases/100 000 rural inhabitants in the year 2008. This cut-off value separates districts with low and high incidence rates and was derived from the continuous frequency distribution of case incidence, since it marks a clear discontinuity in the distribution.
RESULTS
Analysis of monthly air climate features showed that many parameters were significantly different between positive and negative cells in various years. However, most changes lacked a coherent pattern (e.g. the differences were not consistent from year to year). The only group of variables that differed significantly every year between positive and negative cells comprised the mean temperatures from November of the previous year to February of the current year (P=0·001 in every calculation). In summary, mean air temperatures in the winter prior to the case-reporting year differed consistently every year, between positive and negative cells. No differences were found in ground climate features, including evaporation (P=0·356), humidity (P=0·411), mean ground temperature (P=0·121), and mean temperature below the surface (P=0·654). A similar lack of coherence was observed for NDVI values; many significant differences were detected in monthly variables, but we failed to observe any regular pattern being reproduced from year to year (Fig. 1). Both monthly mean and seasonal peak of NDVI values are equivalent in the years and between the positive and negative sites. No significant differences were found for the length of growing season, the summer season, or NDVIa. However, a coherent pattern was observed for winter accumulated temperatures (Fig. 2), in agreement with previous results above. Cells where CCHF cases had been reported always had a lower accumulated temperature than did sites where the tick existed but from which disease had not been reported (P=0·001). We further examined trends in monthly air temperatures and NDVI values in the positive and negative cells. Analysis showed that temperature was increasing and NDVI falling for most of the territory analysed. ANOVA analysis indicated that positive cells were associated with a smaller trend of increase in temperature (slope: 3·86 for positive cells, 4·12 for negative cells; P=0·004) and a higher trend towards decrease in NDVI values (slope: −0·009 for positive cells, −0·007 for negative cells; P=0·005).
Fig. 1. Monthly mean temperature values (°C) and NDVI values (unitless values) for positive and negative cells between 2002 and 2008.
Fig. 2. Values of winter accumulated temperatures (degrees/days), and duration of the growth and summer seasons (days), for 2002–2008, shown separately for positive and negative cells.
Significant differences were found in the rate of recruitment in positive and negative cells. CCHF cases were always associated with cells with a high recruitment value and there was a direct relationship between the presence of the disease and total recruitment in the cell. This unambiguously associated the disease with high tick recruitment values derived from habitat fragmentation features. Figure 3 shows the difference in recruitment values in negative and positive cells during the period 2003–2008. Positive cells were always associated with high recruitment values whereas negative cells always had low recruitment values of <10 (P=0·001 for every year). Figures 4 and 5 provide an overview of the spatial distribution of both CCHF cases and recruitment values for the period 2003–2008. These figures provide a visual link between the critical vegetation pattern detected by satellite and CCHF case distribution, and demonstrate a firm causal basis for retrospective mapping of the disease using only remotely sensed features. There was a direct relationship between case incidence (cases/100 000 rural inhabitants) and recruitment values. Figure 6 displays the evolution in time of case incidence and recruitment in positive cells for districts of either high or low case incidence (Fig. 6a, b, respectively). Recruitment experienced a sharp increase beginning in 2004, whereas case incidence began to dramatically rise after 3 years of continuous increments in recruitment values. Districts where the case incidence remained <20 in 2008 always showed small but increasing values of recruitment. It is interesting to note that, in districts with lower case incidence rates, a distinct rise in case frequency in 2008 was associated with a continuous increase in recruitment in the three or four previous years.
Fig. 3. Recruitment values (unitless) in negative and positive cells for 2003–2008.
Fig. 4. Spatial distribution in Turkey of cells reporting CCHF cases (left panels) and recruitment values (right panels) for 2003–2005. Darker grey indicates higher recruitment values (unitless values).
Fig. 5. Spatial distribution in Turkey of cells reporting CCHF cases (left panels) and recruitment values (right panels) for 2006–2008. Darker grey indicates higher recruitment values (unitless values).
Fig. 6. The temporal trend in recruitment value (▪) and case incidence (
) (cases/100 000 rural inhabitants) in districts of Turkey where CCHF cases were reported during 2000–2008 (recruitment) and 2003–2008 (case incidence). (a) Includes only districts with a case incidence >20 in 2008. (b) Shows the same parameters for districts where disease has been reported but where case incidence values were <20 in 2008.
DISCUSSION
The recent, and still spreading, epidemic of CCHF in Turkey, increasing from around 20 cases annually in 2002 and 2003 [Reference Karti19] to over 1300 in 2008, has made it necessary to understand the factors responsible for disease spread. In this study we used an accurately georeferenced dataset of cases, including the geographical location of each tick bite. The nosocomial infections comprise a small proportion of CCHF cases in Turkey, and the principal transmission route is tick bites, the distribution of which is commonly restricted by both climate and vegetation. Although somewhat contradictory opinions on the possible effect of climate on epidemic spread have appeared [Reference Estrada-Peña4, Reference Purnak, Selvi and Altundag10, Reference Randolph and Ergönül20], a comprehensive study incorporating the environmental factors needed to define the hypothetical niche for the disease has never been undertaken. Coherence in observations across space and time is essential to such a statistical approach [Reference Randolph21]. Temperature and rainfall are both important environmental drivers of disease transmission by ticks, regulating tick development and population seasonality. Small variations in these features at critical points in the year may result in a radically different framework for disease transmission [Reference Randolph22]. Thus, the range of climate and vegetation features across a large area must be analysed before any effect attributable to a specific feature can be unambiguously identified. Any environmental niche for the disease must be characterized, cross-checked for persistence in space and time, and unambiguously ascribed to areas reporting CCHF cases, as performed in studies of other diseases (i.e. [Reference Randolph21]). Previous work on the recent rise of tick-borne encephalitis has shown that climate may play a secondary role driving human activity changes [Reference Randolph23].
CCHF outbreaks have developed against a background of climatic and environmental changes advantageous for the survival of large numbers of ticks and hosts [Reference Hoogstraal5]. In the former Soviet Union, environmental changes included wartime neglect of agricultural lands, wide-scale collectivization of agriculture, and changes in pasture patterns, converting flood plains to farmland. During World War II, after the occupation of Crimea, normal agricultural activities were disrupted and the common sport of hunting European hares was abandoned. When Soviet troops reoccupied the hilly Crimean steppes in 1944, hares had become superabundant and neglected pastures were overgrown with weeds providing a good environment for the tick hosts [Reference Hoogstraal5]. A similar explanation was suggested as the cause of the initial outbreak in Turkey [Reference Karti19]. However, this hypothesis fails to explain why the disease continues to spread in Turkey, given that abandoned lands have already been re-occupied by farmers and that the disease is continuously spreading westwards and southwards.
Current theories of disease spread seek to explain dissemination through a network of connected foci [Reference Schimit and Monteiro13, Reference Wilson18]. The physical characteristics of such a network define the ability of a disease to extend across geographical space. Using remotely sensed information, we obtained fragmentation patterns of vegetation over Turkey, and then derived recruitment values, which are a surrogate of tick abundance according to habitat fragmentation and connectivity plus local climate suitability [Reference Estrada-Peña15–Reference Estrada-Peña17]. These features lead to three important results. First, the sites of highest habitat fragmentation within the known range of the tick were spatially correlated with regions of case reporting. Second, the districts in central Anatolia with highest case incidence showed a trend towards increasing recruitment values, linearly related to case incidence rates in the region. Third, districts where the disease is reported but case incidence rates have remained relatively low displayed consistently low values of recruitment, but these were always higher than areas where disease remained unreported.
It is of interest that neither changes in mean NDVI values nor seasonal changes of this figure in positive and negative sites are significant. The trend in the whole country, as detected in the current study, is towards a small and not significant decrease in NDVI values. As we were more interested in long-term trends, focused on a wide geographic area, our source imagery for NDVI values lacked the adequate temporal and spatial scales necessary to detect the reasons behind changes in habitat patterns. We thus ignore the reasons driving such a fragmentation of the habitat into smaller patches, which appeared as the main driver of the increase of case incidence in historically endemic zones. We must consider social factors, all of which are difficult to quantify. A continuous population exodus from rural sites to urban areas has been recorded for the latest years in the Anatolia region. Such migration leaves former cultivated areas, and favours the growth of bush and undergrowth, which is a good habitat for wild hosts that can amplify the infected tick populations. Further, there is promotion of ecological regeneration in the region, which also encourages formation of dense shrub and undergrowth as the first steps in forest recovery. Both a dramatic decrease in sheep population (according to the annual census of domestic ruminants in Turkey) and a ban on goat grazing in areas of young forest growth, seem to be behind bush growth in central Anatolia, with the consequent increase in densities of wild animals. Outside the region of central Anatolia, the trend in the satellite imagery series shows degradation of forests, probably due to farming activities, leaving small patches of undergrowth between cultivated land and nearby urbanized areas. This unprecedented increase in habitat fragmentation leads to greater contact between tick and reservoir, and therefore to a virus amplification cycle. The increased accessibility of these foci to humans (shrub areas where ticks are common remain within cultivated areas), further leads to closer contact of humans and infected ticks. The mechanism of an increased habitat fragmentation and a higher contact with humans is similar to that reported for Lyme disease in USA [Reference Allan, Keesing and Ostfeld8, Reference Jackson, Hilborn and Thomas9].
At this stage, two hypotheses remain to be tested. The structure of a growing contact network may allow fast spread of the infected reservoir hosts and/or vectors, by both short- and long-distance invasive events. Short-distance exchange of infected ticks through a network of highly fragmented and well-connected habitat patches would operate in districts with high case incidence, contributing to a massive increase in case incidence. Long-distance invasive events might act to spread infection to places suitable for tick establishment, as the ability of birds to carry feeding tick immatures has been demonstrated [Reference Hoogstraal24]. Such long-distance spread to new areas was observed after 2006, when new foci with high recruitment appeared far from central Anatolia. However, it is likely that infected ticks and reservoir hosts have always been present at a level below the critical threshold for infection of humans.
We found that climate variables derived from both air and ground were poor descriptors of any environmental niche of the disease. Further, no single variable served as an accurate discriminator of areas where disease had been reported, except for winter accumulated temperature. Areas where the disease is present were consistently colder than sites where the disease had not been reported, within the known area of tick distribution over the country. This was a reliable association, from year to year, over the whole territory of our study. Although this is not an unambiguous definition of a disease niche, the finding supports previous reports on the effect of winter temperatures on tick populations. Earlier climate characterization studies in areas of Crimea and Armenia where the tick H. marginatum is present showed that winter was a critical period for population survival [Reference Kondratenko25]. Pioneering research on the topic (reviewed in [Reference Hoogstraal5]) pointed towards a critical minimum winter temperature as a natural control factor limiting survival of overwintering ticks. A recent reappraisal of the issue [Reference Estrada-Peña26] showed that the most important factor for geographical restriction of H. marginatum populations in the area was the lower limit of late autumn and early winter temperatures. These limits affected the ability of engorged nymphs to molt before the sudden fall in winter temperatures and to survive the coldest period as flat overwintering adults. It has been shown for other tick species that interruption of molting by low temperature in winter is usually fatal [Reference Materna, Daniel and Danielova27].
It is thus reasonable to hypothesize a critical range of winter temperatures, within regions preferred by the tick for survival, as a secondary driver of disease emergence. The critical range of temperatures observed for areas of disease, where temperatures were slightly lower than in CCHF-negative sites, could be interpreted from the perspective of the activation temperature of questing adults. As immatures commonly molt in late autumn in the region, a medium to high temperature immediately after molting and cuticle hardening would allow adults to begin to quest. We hypothesize that lower temperatures associated with CCHF-positive areas would result in a lack of activity by adults after molt. These data might support the idea that more abundant overwintering tick populations would be synchronously activated in the following spring, thereby increasing the chance of contact between animal hosts and ticks and amplifying the infection. Additional work on selected field sites would provide new insights about this extreme. It is of interest that neither the length of the growing season nor summer length, both of which would allow a longer transmission period because of an extended tick activity season, was of any significance in the niche definition of zones with disease.
In summary, there are no consistent indicators defining a climate niche for CCHF in Turkey, although winter temperatures may have a influence on the synchronous activation of ticks in the following year. Whereas longer periods of temperatures above a given activation threshold for ticks (i.e. in the growing and summer seasons) have been reported to potentially increase the rate of tick contact with humans [Reference Purnak, Selvi and Altundag10], this can no longer be accepted as the main driver of disease spread. However, increase in habitat fragmentation is unequivocally associated with active foci of disease, geographical disease spread, and case incidence rate. A network-based study of the importance of high tick recruitment zones and the epidemiological significance of the resulting disease network is necessary to define further features of the spread of disease.
ACKNOWLEDGEMENTS
The authors thank the members of Communicable Diseases Unit of the Ministry of Health of Turkey, Mehmet Ali Torunoglu, Bedia Turkyilmaz, and Ahmet Safran, for their tireless support. Parts of this work were supported through the network of the European Community's Seventh Framework Programme under the ARBO-ZOONET grant agreement no. 211757.
DECLARATION OF INTEREST
None.
