Abstract
Babesiosis is an emerging arthropod-borne infection that has been increasing in incidence for the last decade in the northeastern United States. Babesiosis may share features of its landscape epidemiology with other arthropod-borne infections transmitted by the same tick vectors in similar geographic spaces. This study examined 11 years of surveillance data in New York State to measure the relationship between forest fragmentation and the incidence of human babesiosis. Adjusted Poisson models showed that increasing edges of contact between forested land and developed land, as measured by their shared perimeters, was associated with a higher incidence of babesiosis cases (incident rate ratio [IRR]=1.015, 95% confidence interval [CI] 1.01–1.02; p<0.001), even after controlling for the total developed land area and forest density, and temperature and precipitation. Each 10-km increase in perimeter contact between forested land and developed land per county was associated with a 1.5% increase in babesiosis risk. Higher temperature was also strongly associated with increasing babesiosis risk (IRR=1.18, 95% CI 1.10–1.27; p<0.001), wherein each degree Celsius increase was associated with an 18% increase in babesiosis risk. While direct causal conclusions cannot be drawn from these data, these findings do identify a potentially important signal in the epidemiology of babesiosis and suggest that the underlying physical landscape may play a role in shaping points of contact between humans and tick vectors and the subsequent transmission of Babesia microti.
Introduction
Materials and Methods
All incident cases of babesiosis in humans were obtained from the New York State Department of Health (NYS DOH) by county for each of the 11 years from 2000 through 2010 (New York State Department of Health 2012). NYS DOH identifies human babesiosis cases by ongoing surveillance via mandated reporting. Cases were defined by positive blood smear, nucleic acid test, or Babesia-specific antibody titer >256 with an indirect fluorescent antibody (IFA) test for immunoglobulin G (IgG) or total antibody. The latter serology has demonstrated good validity and reliability (Ruebush et al. 1981). All cases require laboratory confirmation at both the New York State Wadsworth Center and the New York City Public Health Laboratory (New York State Department of Health 2010). County-specific estimates were not available for the 5 counties that make up the municipality of New York City, so these counties were not included in this analysis, leaving a total of 57 of the total 62 counties that constitute the state of New York.
Temperature data were obtained from a monthly surface data product from the National Climate Data Center, which manages meteorologic data products for the National Oceanic and Atmospheric Administration (NOAA 2012). Mean monthly temperatures were recorded for each of the 11 years from 2000 through 2010. Given the obvious problem of temperature variation by month, but the potential residual impact of higher temperatures overall throughout the babesiosis season (Spielman et al. 1985), this study took a unique approach to quantifying temperature. Rather than creating discrete panel temperature data for each month, a monthly-weighted cumulative measure of temperature was used to represent the human transmission season as a whole from May through October (Spielman et al. 1985, Vannier and Krause 2012). The mean temperature for each month from May through October was multiplied by its weighting factor, which consisted of the proportion (ranging from 0 to 1) of that month's temperature's estimated contribution to the overall occurrence of babesiosis in humans throughout the season. The relative importance of each month's temperature to the collective season's babesiosis occurrence was determined by the average distribution of heat across each month from May through October. The months of June through September were weighted higher because of their relative importance for exposure to ticks, and because of the demonstrated seasonality observed in annual babesiosis incidence (Ruebush et al. 1981; Spielman et al. 1985; Piesman et al. 1987; Vannier and Krause 2012). The specific weighting was as follows: The mean temperature for May was multiplied by 0.05, the mean temperature for June was multiplied by 0.15, and July, August, September, and October mean temperatures were multiplied by 0.3, 0.3, 0.15, and 0.05, respectively. These weights reflect the relative importance considered for each month's mean temperature to the overall occurrence of babesiosis for each year. The weighted mean temperatures for each month were then summed over all months to create a single seasonal temperature for each county for each year under study, which then has the same unit interpretation in degrees Celcius. To test the validity of this seasonal temperature measure, sensitivity analyses were conducted wherein the analyses were repeated using the mean temperature for each month alone, again from May through October. The associations between temperature and babesiosis, and the temperature-adjusted associations between forest fragmentation and babesiosis, were very similar between the seasonal temperature measure and the individual monthly mean temperatures. Therefore this investigation proceeded with the use of the monthly-weighted seasonal temperature, which is a more comprehensive measure of temperature, as well as being more easily interpretable on a year-to-year basis.
Hydrogeography data were obtained from the United States Geological Survey National Hydrography Dataset at the midpoint of the decade under current analysis (US Geological Survey 2012). The area, in square miles, of all lakes, ponds, swamps, and marshes was used to designate the total surface water per county. The area per county was obtained by calculating the geometry in a geographic information system (GIS). There were no differences observed across types of surface water, i.e., lakes and ponds versus swamps and marshes, so all freshwater bodies of surface water were combined for easier presentation and interpretation of results.
Forest fragmentation and land use was determined by data obtained from the Unites States Department of Agriculture's Forest Service (US Forest Service 2012). The Northeastern Forestry Inventory Analysis classified each 30-meter pixel as forested or developed, and the following fragmentation metrics were created. For each county, the proportions of the total land area that are forested and developed were computed; the length, in kilometers, of all edges of forested land parcels that are shared with developed land parcels was summed; the average forest perimeter patch length, in kilometers, average forest patch area, in hectares, and largest forest patch area, in hectares, were all summed; and the forest density was the average length of the edge of a forest land parcel divided by the parcel's area.
Statistical analysis
Random effects models were fit using Poisson regression for panel data, and the results were used to identify the associations between increasing contact between developed and forested land parcels and babesiosis cases per county across the 11 years of babesiosis surveillance, controlling for landscape factors, temperature, and precipitation. Poisson regression is a technique that models count data, which are assumed to follow a Poisson distribution. In the current study, the Poisson regression is modeling the counts of human infections per the total county population density (population/square kilometer) per year, and therefore estimates the incidence rate. As such, incidence rate ratios are reported in the results. These models were used to assess both the unadjusted bivariate associations between each landscape and climate factor and babesiosis (Table 1, below), and the adjusted associations between each factor demonstrating a significant (p<0.05) or borderline significant (p<0.1) bivariate association and babesiosis (Table 2, below). Precipitation showed no significant or borderline significant bivariate association with babesiosis, but was nevertheless included in the fully adjusted model in Table 2 to control for any potential confounding.
These associations are derived from Poisson models, with county and year as the panel levels. Each factor was entered into the Poisson model as the single independent predictor of babesiosis cases.
Temperature is the monthly-weighted seasonal cumulative temperature in degrees Celsius
IRR, incident rate ratio; CI, confidence interval.
The full model is derived from multiple Poisson regression, with county and year as the panel levels. The associations between babesiosis occurrence and the independent variables are adjusted for the other variables in the model.
Temperature is the monthly-weighted seasonal cumulative temperature in degrees Celsius.
IRR, incident rate ratio; CI, confidence interval.
Finally, a map was created to demonstrate the spatial prediction of babesiosis incidence by the landscape and climate factors. The map is comprised of 2 layers. The first layer, consisting of graduated symbols of observed babesiosis cases from 2000 to 2010, is superimposed over the second layer, which is a choropleth map of babesiosis cases predicted by forest fragmentation (shared parcel edges between developed and forested land, total proportion of developed land, and forest density), temperature, and precipitation. ArcGIS was used for all mapping procedures and for calculating the Moran Index to identify whether spatial clustering of the residuals was present (ESRI), and STATA v.11 was used for the Poisson regression modeling specific to panel data, with the xtpoisson command procedure (STATA Corp., College Station, TX).
Results
During the 11-year period from 2000 to 2010, 1610 human cases of babesiosis were reported across the state of New York, excluding the 5 counties of New York City. These cases were distributed across 26 of the remaining 57 counties, ranging from just single sporadic incident cases in 10 counties to hundreds of cases in Dutchess (285), Suffolk (936), and Westchester (232) counties. The median shared perimeter between forested and developed land was 624.0 kilometers and was greatest in Suffolk (4467.71 km) and Westchester (3667.554 km), and least in Hamilton (231.8569 km), Yates (144.1664 km), and Wyoming (138.374 km) counties.
Table 1 displays the unadjusted bivariate associations between babesiosis incidence and each of the landscape and climate factors. Forest fragmentation was associated with the occurrence of babesiosis. As hypothesized, increasing edges of contact between forested and developed land, as measured by the distance of the total shared edges between these land classes, was associated with increased incidence of babesiosis. Each 10-km increase in shared edges between distinct land parcels corresponded to a 1.6% increase in the risk of babesiosis (incident rate ratio [IRR]=1.016, 95% confidence interval [CI]=1.01–1.02; p<0.001). In addition, increasing forest density among forested land parcels was associated with increased risk of babesiosis, wherein each meter/hectare increase in density was associated with a 6% increase in the incidence of babesiosis cases. As expected, the overall percentage of developed land was associated with increased risk. Each 1% increase in the total developed land parcels per county corresponded to a 9% increase in the incidence of babesiosis, though this association was only of borderline significance (IRR=1.09, 95% CI 1.0–1.19; p=0.07). Finally, a strong relationship between seasonal temperature and babesiosis emerged, with each degree Celsius increase corresponding to a 20% increase in risk (IRR=1.20, 95% CI 1.12–1.28; p<0.001). The average forest patch perimeter and the average and largest forest patch areas were not associated with babesiosis. Neither precipitation nor surface water area were associated with babesiosis incidence over the study period.
Table 2 presents the adjusted associations between babesiosis and each of the landscape and climate factors. As above, increasing edges of contact between forested and developed land were significantly associated with increasing babesiosis incidence, even after controlling for the overall proportion of developed land, the density of forested land, and climate. Each 10-km increase in shared edges between forested land and developed land corresponded to an adjusted 1.5% increase in risk (IRR=1.015, 95% 1.01–1.02; p<0.001). Moreover, neither the total proportion of developed land (IRR=0.97, 95% CI 0.94–1.01; p=0.133) nor forest density (IRR=1.04, 95% CI 0.99–1.09; p=0.127) was still significantly associated with babesiosis. Finally, temperature remained significantly associated with babesiosis (IRR=1.18, 95% C.I. 1.10–1.27; p<0.001), while precipitation remained unassociated.
Figure 1 shows graduated symbols, representing observed cases of babesiosis, superimposed over a choropleth map of babesiosis cases predicted by forest fragmentation (shared parcel edges between developed and nondeveloped land, total proportion of developed land, and forest density), temperature, and precipitation (i.e., cases predicted by the adjusted Poisson model described in Table 2). The map demonstrates the high concentration of observed and predicted cases in the Hudson Valley and on Long Island. Furthermore, this map highlights the close agreement between the observed babesiosis cases at the county level and those predicted by forest fragmentation and climate. The Moran Index was used to identify global spatial clustering of the residuals from the full Poisson model. The index was equal to −0.033 (p=0.77), which indicated no substantive spatial clustering.

Map of the cumulative observed and predicted babesiosis cases in New York State from 2000 to 2010. Graduated symbols of observed babesiosis cases, with increasing size symbolizing increasing number of cases, are superimposed over a choropleth map of babesiosis cases predicted by forest fragmentation, hydrogeography, and temperature. In the choropleth layer of the map, darker shading indicates a higher concentration of predicted babesiosis cases. Color images available online at
Discussion
In this investigation, the incidence of babesiosis was significantly associated with increasing contact between developed and forested land parcels. Each 10-km increase in shared edges between developed land and forested land corresponded to a 1.5% increase in risk, even after adjusting for the total proportion of developed land, forest density, temperature, and precipitation. This is the first report to describe the landscape epidemiology of human babesiosis in the northeastern United States, and, more specifically, to describe the association between forest fragmentation and incidence.
Human babesiosis has expanded dramatically over the last 2 decades. Only 136 cases were reported across all counties in New York State between the years 1982 to 1991. And 134 of these cases came from just 1 county (Suffolk), whereas the other 2 were imported (Meldrum et al. 1992). In contrast, the current report shows that the number of cases has increased more than 10-fold in the period from 2000 to 2010, with cases now being reported in almost half of New York State counties. A recent report by the Centers for Disease Control and Prevention (CDC) published in the Morbidity and Mortality Weekly Report showed that human endemicity is now well established in New York State, with 361 cases in 2011, which is the highest number of cases out of all the states in which babesiosis is reportable (Herwaldt et al. 2012). This same report identified the northeastern and upper Midwestern states as the current areas of greatest incidence in the country, which follow the same geographic clustering as the hyperendemic areas of Lyme disease and, ultimately, the primary vector of both, Ixodes scapularis (Herwaldt et al. 2012, Marquardt 2004). Moreover, the incidence of babesiosis has been steadily increasing since 2000 throughout the lower Hudson Valley (Joseph et al. 2011), which is corroborated here.
While the relationship between forest fragmentation and human babesiosis may be newly described, the importance of forest fragmentation for Ixodes ecology in general, and for other tick-borne infections in particular, is certainly not new (Guerra et al. 2002; Bunnell et al. 2003). An investigation of a majority section of the state of Maryland identified such forest fragmentation as a key factor in predicting the incidence of Lyme disease (Jackson et al. 2006). This investigation showed that contiguous undisturbed forested land parcels are protective against the occurrence of human Lyme borreliosis, while fragmented forest landscapes, in which residentially developed parcels shared boundaries more extensively with nondeveloped land, were associated with increased occurrence of cases (Jackson et al. 2006). Forested land lost to development not only reduces the absolute area of forested land, but also increases the perimeter-to-area ratio of forested land parcels. Thus, boundary contact between developed and nondeveloped land parcels increases as does the potential contact between vectors and human hosts (Lambin et al. 2010). In fact, 1 of the earliest investigations of human babesiosis in the northeastern United States described a greater number of cases identified among residential communities with closer proximity to forest habitat (Ruebush et al. 1981). To provide more substance for the potential shared ecology suggested above between diseases that share the same tick vector, i.e., Lyme disease and babesiosis, and, subsequently, to demonstrate the robusticity of the association between forest fragmentation and tick-borne disease in New York State, this association was also analyzed specifically for Lyme disease incidence over the same 11-year period and across the same counties of New York State. The results show that this association was quite similar between the 2 different tick-borne infections. Each 10-km increase in shared edges between forested and developed land parcels corresponded to a 1.8% increase in risk for Lyme disease (IRR=1.018, 95% CI 1.01[–1.02; p<0.001) adjusted for the same land cover and climate factors included in the babesiosis model, which, as described above, identified an adjusted 1.5% increased risk.
The life cycle of the tick vector employs a complex incomplete metamorphosis that requires multiple mammalian hosts to complete. Larvae acquire their blood meal primarily from the white-footed mouse, which also serves as the sylvan reservoir for B. microti (Marquardt 2004, Vannier and Krause 2012). After overwintering, the nymphs emerge in the spring and also favor the white-footed mouse for their second blood meal. However, humans are often accidental hosts at this stage of tick development. Transmission to humans can be facilitated during this time due to the small size of the nymphs, which allows feeding ticks to go unnoticed, thus allowing sufficient time to transmit the parasite. After taking the blood meal, the nymphs molt into adults later in the summer or fall. The adults seek larger mammals, most typically white-tailed deer, for their final blood meals before females overwinter again and lay their eggs the following spring (Marquardt 2004, Vannier and Krause 2012). The adult stage represents another point of potential exposure and accidental transmission to human hosts. Because of the vector's extended life cycle with the exploitation of multiple hosts in sylvan settings, and because of the easy transition from sylvan to domestic hosts for I. scapularis, greater contact between nondeveloped and developed land permits a fairly efficient conduit to human infection from the natural reservoir (Keesing et al. 2010, Lambin et al. 2010). Favorable ecology for transmission to humans is further amplified by changes in the population density and range of larval/nymphal hosts (i.e., white-footed mice) and adult hosts (white-tailed deer), which is another consequence of human-altered landscapes (Spielman et al. 1985).
Temperature has also been demonstrated to influence exposure to the tick vector, as well as the frequency of tick bites, because warmer weather can lead to more time spent in tick habitat and can also lead to increases in tick population density (Spielman et al. 1985), both of which can lead to higher levels of human exposure. Although increased exposure due to increased time spent outdoors is very difficult to assess, increased exposure due to increased numbers of ticks in the environment can be assessed more readily. One particularly large investigation surveyed I. ricinus in 2082 distinct geographic sites across 18 different habitats in Spain and found that higher temperatures were strongly associated with higher tick densities (Estrada-Peña 2001). Nevertheless, such data are still limited in that only the tick populations themselves are measured, rather than direct exposure episodes between the tick and human populations. In addition, the effect of temperature on human infection can be mediated directly by the parasite itself. B. microti sporozoites mature faster in warmer conditions, and thus may require substantially less time after attachment for the tick to transmit the infection to the host (Spielman et al. 1985). Our results agree with these previous findings, with each degree increase in temperature corresponding to an adjusted 18% increase in babesiosis risk. However, it is important to recognize that neither the associations presented in this current study, nor those identified in earlier studies, can establish causation between temperature and babesiosis because none of these assessed direct human tick exposure experimentally with temperature.
This study has some important limitations that must be recognized. The outcome and all the landscape variables are aggregated at the county level for each year. Because of this, the specificity of the associations is diminished and the resolution of the ecologic and epidemiologic picture is necessarily coarse. Undoubtedly, more spatially and temporally localized microgeographic phenomena operate with respect to specific boundaries between developed and nondeveloped parcels in the landscape. Hyperlocal ecotonal ecology may be particularly relevant, but not discernible from the current investigation. This report did not address such localized features, thus the results are only suggestive, not conclusive. Furthermore, any direct causal claims would require a controlled experimental design involving human, animal host, and vector sampling for babesiosis georeferenced precisely with the land cover. Moreover, this investigation has used a single fixed measure of land cover to explore the geometry of forest fragmentation and its potential association with babesiosis incidence over an 11-year period. Under this study's construct, babesiosis occurrence (and climate) changed over time, whereas forest fragmentation did not. As such, the impact of the nuance of change in forest fragmentation over the study period could not be assessed and may, in fact, be relevant to the epidemiology of babesiosis. Finally, the results described derive from the state of New York, and so the interpretations are limited to the northeastern United States.
In conclusion, this report shows a dramatically increasing incidence of babesiosis in the state of New York between 2000 and 2010. Moreover, this investigation identified a strong association between this rising incidence and the extent of shared boundaries between adjacent developed and forested land parcels. This association was independent of other elements of land use and climate. Finally, temperature was also strongly associated with babesiosis, suggesting that multiple features of the physical landscape may be relevant to babesiosis occurrence in humans. Given the nature of the county-aggregated data, no causal conclusions can be made with respect to forest fragmentation, or temperature, on the transmission of B. microti. Instead, this study serves as a surveillance instrument that has been used to identify a potentially important signal with respect to land cover geometry and the landscape epidemiology of babesiosis. On the basis of these findings, this report does not posit definitive etiologic pathways for human babesiosis, but rather suggests more localized follow-up surveys of vector populations and human babesiosis cases in areas of both fragmented and homogeneous land cover to better understand whether specific geometries in the landscape provide active contact between otherwise disparate ecologies and, thus, facilitate vector transmission of this parasite.
Footnotes
Author Disclosure Statement
No competing financial interests exist.
