Lyme disease is the most commonly reported vector-borne disease in the United States (Hinckley et al., 2014; CDC, 2015b). Human infection occurs primarily in the North-eastern and North Central United States, but is not exclusive to these areas due to constant geographic expansion of the principal vector Ixodes scapularis (CDC, 2015a; Mead, 2015). Human risk of exposure is uncertain in many areas where Lyme disease is considered an emerging issue due to variability in disease manifestation (Sanchez et al., 2016; Steere et al., 2016), limitations of diagnostic tests (Aguero-Rosenfeld and Wormser, 2015) and variable infection prevalence ascertained through local tick surveys (Guerra et al., 2001; Duncan et al., 2005; Hamer et al., 2009). Tick surveys among companion animals have been identified as a potentially representative surveillance methodology for estimating human risk within geographic areas where Lyme disease is an emerging concern (Millen et al., 2013). Animals such as dogs, cats and horses are potentially effective sentinel populations due to increased likelihood of tick infestation and close association with their human owners (Anderson, 1989). However, few studies to date have evaluated the utility of companion animal sentinel programmes as public health tools with spatial analytic approaches to study the spatial structure of tick-host disease transmission dynamics (Killilea et al., 2008; Abdullah et al., 2016; Tulloch et al., 2017).
Humans and companion animals are dead-end hosts for Borrelia burgdorferi, the causative agent of Lyme disease (Little et al., 2010; Radolf, 2012). Exposure to tick-borne pathogens occurs when a host comes in contact with actively questing ticks in the environment. Presence of infected ticks are influenced by environmental factors such as availability of competent reservoir hosts and climatic factors that drive tick survival and pathogen propagation (Eisen et al., 2012; Brinkerhoff et al., 2015). Domestic pets, such as dogs, are significant risk factors for disease acquisition, and six times more likely to be exposed to infected ticks due to the increased potential exposure time in tick habitat (Jones et al., 2002; Hamer et al., 2009). Sampling companion animal populations is an effective approach to estimate human risk in areas where Lyme disease may be newly endemic.
Popular concern for Lyme disease has prompted many animal sentinel studies across the United States and the United Kingdom (Goosens et al., 2001; Johnson et al., 2004; Duncan et al., 2005; Glickman et al., 2006; Hamer et al., 2009; Smith et al., 2012; Abdullah et al., 2016). These studies primarily focused on ticks and associated pathogens collected from domestic dogs because of their ability to produce antibodies to B. burgdorferi, attainable travel history information, and frequency of outdoor exposure (Burgess, 1986; Eng et al., 1988; Lindenmayer et al., 1991; Mather et al., 1994; Walker et al., 1998; Duncan et al., 2005; Little et al., 2010). Methodologies for pathogen detection are restricted to serologic tests such as enzyme-linked immunosorbent assay (ELISA) and Western blot, as well as direct polymerase chain reaction (PCR) of ticks collected from domestic dogs (Goosens et al., 2001; Johnson et al., 2004; Duncan et al., 2005; Glickman et al., 2006; Hamer et al., 2009; Smith et al., 2012). Studies suggest that serology screening is an effective approach in Lyme-endemic areas (Bowman et al., 2009; Mead, 2011) but may not be as sensitive an indicator for estimating prevalence of infected ticks as direct screening of ticks via PCR within areas of low tick density (Hamer et al., 2009). Previous studies also suggest variable sample sizes may result in potentially bias results in localities with few participating veterinary practices (Johnson et al., 2004).
Some studies have shown that success in detection of B. burgdorferi in dogs or ticks collected from dogs correlates with human Lyme disease incidence and tick feeding behaviour within their respective localities (Eisen et al., 2006; Little et al., 2010). This conclusion, along with recognition that tick survival and disease transmission are spatial processes (Eisen et al., 2012), recommend analytic spatial approaches to identify areas of increased human risk. To date, no known study has investigated geographic associations between Ix. scapularis collected from sentinel companion animal surveillance and human Lyme disease cases comparing non-spatial and spatial regression approaches particularly in rural states such as West Virginia. Additionally, few or no available studies have utilised spatial diagnostic tools to identify areas where increased recruitment may be necessary to optimise Lyme disease surveillance or considered sentinel potential for cats in regards to Ix. scapularis and Lyme disease. The objectives of our study were to 1) conduct spatial regression to identify geographic associations between animal species-specific Ix. scapularis sentinel and human Lyme disease case data and to 2) investigate local clustering of regression residuals to identify counties were additional effort may be needed to optimise surveillance efficacy.
Materials and Methods
West Virginia is in a mid-latitude, temperate area, with a humid continental climate that is warmer in the lower elevation parts of the state in the south-western and eastern counties. Temperatures are near 0°C in winter and 25°C in the summer, with colder temperatures at higher elevations (National Climate Data Center, 2017). The state receives 80 to 140 cm of precipitation annually, most of which along the western side of the Allegheny Front. West Virginia is a mountainous state, averaging about 460 m above the mean sea level (MSL), with the lowest areas along the Ohio River on the western border and the Potomac River along the north-eastern border. The Allegheny Front runs from the southern to the north-eastern part of the state, and has the highest elevations, reaching 1,482 m above MSL at Spruce Knob (West Virginia GIS Technical Center, 2011). The landscape is dominated by temperate forests (West Virginia GIS Technical Center, 2017). The state is predominantly rural, with the largest cities being Charleston at 49,138 residents and Huntington at 48,113 residents (US Census Bureau, 2017). Annual confirmed human case counts within West Virginia have increased two-fold between 2005 (n=161) and 2016 (n=388) (CDC, 2015a; WVDHHR, 2016). Geographically, West Virginia is bordered on the East by the highly endemic (for Lyme disease) states of Pennsylvania, Maryland and Virginia and on the west by the low-incidence states of Ohio and Kentucky. States which have high incidence status have had ≥10 confirmed cases of Lyme disease per 100,000 persons for the last three reporting years (CDC, 2017). In 2017, West Virginia met this criterion based on 2014 to 2016 human surveillance data. In addition to human cases reported, companion animal sentinel surveillance data have been collected from local licensed veterinary practices within the state since 2013 (WVDHHR, 2017).
Data acquisition and management
Human Lyme disease case data in West Virginia were obtained for 2014-2016 from the West Virginia state health department (WVDHHR, 2016). Human surveillance data were limited to confirmed cases (based on the 2011 Council of State and Territorial Epidemiologists Lyme disease case definition) with reported infections acquired within their home county. Briefly, cases assigned to a confirmed status need to have either 1) Erythema migrans (EM) and known exposure; 2) EM rash, laboratory evidence, and no known exposure; or 3) a case with at least one late stage manifestation and laboratory evidence (CDC, 2011). Known exposure is defined as having been in a tick habitat defined as a wooded, brushy or grassy landscape within a county endemic for Lyme disease within 30 days of EM onset (CDC 2011). Laboratory evidence for diagnosing a confirmed case may consist of 1) a B. burgdorferi positive culture; 2) a positive IgM response within 30 days of symptom onset or a positive IgG response at any time; 3) a singletier positive IgG response; or 4) cerebrospinal fluid (CSF) antibody positive enzyme immunoassay or immunofluorescence assay when the concentration of antibody is higher than it was in serum (CDC, 2011). Restricting raw case counts to confirmed status-limited biases associated with travel related disease acquisition (Szonyi et al., 2015) led to consistency of medical, laboratory and epidemiologic evidence (Li et al., 2014) and included recent exposure to a tick habitat (CDC, 2017). Crude measures of human Lyme disease frequency were calculated by dividing the sum of 2014-2016 county-level case counts by 2015 U.S. census West Virginia county-level population estimates (US Census Bureau, 2015) to reflect county-level incidence throughout the state during the period 2014-2016 (Gertsman, 2013).
Invitation letters and information regarding our project were sent to every licensed veterinary practice in the state. Participating practices received quarterly updates using ArcGIS story mapping methodologies (ESRI, Redlands, CA, USA). No practices were incentivised other than providing yearly participation certificates. Companion animal surveillance data were obtained from the West Virginia state health department from 2014-2016 and contained yearly county level information relating the diversity and abundance of different tick species removed from animals by veterinarians participating in the convenience sample. Ticks removed by veterinarians were mailed in secured envelopes with a survey noting animal and practice location, prior tick borne disease testing results, use of tick prevention, animal travel history and relative exposure to a tick habitat (WVDHHR, 2017). Ticks were identified to species, sex and life-stage by trained staff using appropriate taxonomic keys (Keirans and Litwak, 1989). Observations were sorted by county and Ix. scapularis removed from animals were summed across individuals in SAS 9.4.1 (SAS Institute, Cary, NC, USA) to produce four variables representing the total and average number of Ix. scapularis collected from dogs and cats by county. Zeros were assigned to any county with no Ix. scapularis collection data, as a fundamental assumption in our analysis was that all practices received equal opportunity to participate. This assumption was in fact key in maintaining a sample size of counties large enough to conduct the exploratory analysis.
Human population and case estimates were paired with average number of Ix. scapularis estimates from dogs and cats by respective reporting county in Microsoft Excel spreadsheets that were eventually joined in ArcMap 10.4.1 (ESRI) to a U.S. Counties and Localities shape file obtained U.S. Department of Agriculture National Resources Conservation Service Geospatial Gateway (USDA:NRCS) (USDA NRCS, 2015).
Regression and spatial analysis
Separate ordinary least squares (OLS) regression models were performed using average number of Ix. scapularis collected on dogs and cats as the independent variable and the total sum of county-level human Lyme disease incidence per 100,000 population over the study period as the dependent variable. Rook and queen spatial contiguity weights matrices were visualised within connectivity histograms to characterise neighbour structures (Hendricks and Mark-Carew, 2017). Contiguity based spatial weights were utilised within our spatial analysis to satisfy regularity conditions and to reduce the potential for introducing heteroscedasticity (Anselin, 2002). A queen contiguity weight was chosen to increase the number of neighbours evaluated during the spatial analysis (Anselin et al., 2006). Univariate Local Moran’s I and local indicators of spatial autocorrelation (LISA) maps were calculated to identify the extent of spatial autocorrelation within the residuals (Tiefelsdorf, 2000). Lagrange multiplier values were calculated during the OLS regression utilising the previously defined queen spatial weights matrix as a model selection tool in comparing spatial lag and error regression models (Anselin et al., 2006; Matthews, 2006). Residuals from the spatial model were visualised using univariate Local Moran’s I LISA maps to display remaining spatial dependence. Akaike information criteria (AICs) between OLS and spatial autoregressive models were compared to examine model goodness of fit. All regression and spatial analyses were conducted in GeoDa 1.8 with a significance level of 0.05 (GeoDa Center, Tempe, AZ, USA; https://edirc.repec.org/data/gcasuus. html).
Regression and spatial analysis
County-level incidence of confirmed human Lyme disease for the entire period ranged from 0-388 with an average rate of 30.3 [standard deviation (SD)=65.4] per 100,000 persons. Of the 172 total veterinary practices located in West Virginia, 36% (n=62) participated in the study in at least one of the three years. The total number of unique identification numbers recorded for dogs in the study was 1,305, of which 26% (n=349) had at least one Ix. scapularis submitted. The total number of unique identification numbers for cats in the study was 363, of which 59% (n=213) had at least one Ix. scapularis submitted. These estimates serve as a crude estimation of the animal population size, but do not reflect the true population as they do not control for multiple site visits per animal. County-level raw ticks counts and incidence of human Lyme disease per 100, 000 persons are displayed together in Figure 1. County-level Ix. scapularis collections ranged from 0-93 for dogs with a mean number removed of 11.2 (SD=18.8) and 0-37 for cats with a mean number of 5.6 (SD=9.30) removed over the duration of the study period.
OLS regression parameters were significant (F=3.91, P<0.001) for the dog-specific model, and not significant (F=0.91, P=0.36) for the cat-specific model. Univariate Local Moran’s I indicated significant clustering among OLS dog (I=0.27, pseudo P=0.002) and cat (I=0.50, pseudo P<0.001) regression residuals. Lagrange multiplier (LM) values were significant for dog spatial lag (F=24.0, P<0.01) and dog spatial error (F=8.19, P=0.04) as well as cat spatial lag (F=30.1, P<0.01) and cat spatial error (F=27.2, P<0.001) models. The spatial lag autoregressive model was chosen based upon higher robust LM values (30.0 vs 14.0 for dogs and cats (12.8 vs 9.92) (Anselin et al., 2006) and because the spatial lag model adjusts for substantive rather than nuisance spatial dependence (Matthews, 2006). Average number of Ix. scapularis collected per dog remained significant (F=3.04, P=0.002) in the spatial lag model and identified significant positive spatial dependence (rho=0.66, P<0.001). The average number of Ix. scapularis collected per cat remained non-significant (F=-0.14, P=0.89) within the spatial lag model. However, a significant (rho=0.73, P<0.001) strong positive spatial dependence was identified. Univariate Local Moran’s I indicated significant dispersal among canine (I=-0.14 pseudo P=0.02) and non-significant dispersal among feline (I=-0.04 pseudo P=0.40) spatial lag residuals. Parameter estimates for animal specific regression models are summarised in Table 1. Model goodness of fit, measured by a reduction in AIC values, was improved for dog (AIC=582 vs AIC=605) and cat (AIC=591 vs 618) models utilising spatial lag regression. Figure 2 shows LISA maps based on OLS and spatial lag model residuals for dogs and cats. High-high ranking indicates counties where human Lyme disease cases per 100,000 persons are higher than what would be expected after adjusting for covariates. Conversely, low, high-low and low-high rankings indicate counties where human Lyme disease cases per 100,000 population are lower; higher surrounded by lower; or lower surrounded by higher than what would be expected after adjusting for covariates.
Exploratory statistical approaches are routinely applied as methods to evaluate associations between animal sentinel and human disease data. Our study is the first known animal sentinel and human Lyme disease surveillance study to compare and contrast OLS and spatial autoregressive techniques using a combination of visualisation and statistical methodologies, such as Moran’s I, in an exploratory approach. Descriptive analysis identified a higher frequency of Ix. scapularis among veterinarian cat submissions, 59% vs 26% for dogs. Despite this difference, more Ix. scapularis were removed from dogs. Significant associations between Ix. scapularis collections from dogs were identified in OLS and spatial lag models. This finding is consistent with previous studies which identified dogs as potentially effective sentinel animals for monitoring human Lyme disease due to their close association with pet owners and therefore increased probability for tick exposure (Anderson, 1989; Jones et al., 2002; Hamer et al., 2009). Conversely, no statistical association was identified between cat Ix. scapularis collections and human disease in either OLS or spatial lag autoregressive models. This finding is in contrast to recent work by Tulloch et al. (2017), which identified cats as potentially effective sentinel populations for tick-borne diseases. Discrepancy among study findings is potentially attributable to a substantially lower number of cat submissions compared to dog submissions, or an artefact of our decision to exclude non- Lyme disease tick vector species within the analysis. Further research is warranted as interactions between ticks, disease agents, and cat hosts is poorly understood (Krupka and Straubinger, 2010).
Further analysis of model residuals via local indicators of spatial autocorrelation permits additional insight into surveillance system effectiveness by providing a sound diagnostic methodology for identifying spatial dependence within the regression models (Tiefelsdorf, 2000; Burt et al., 2009). Spatial dependence, also termed spatial autocorrelation, occurs when the assumption of independence among dependent variables is violated, or when the regression residuals themselves display spatial dependence (Burt et al., 2009). In either case, parameter estimates from OLS regression are inefficient. Local indicators of spatial autocorrelation (LISA in Figure 2), calculated for OLS regression residuals, identified significant clustering of model over- and under-estimation throughout West Virginia. Incorporation of the spatial lag model resulted in a reduction in the number of counties regarded as highhigh (model under-estimation) and low-low (model over-estimation) among dogs and cats. However, residual spatial dependence was identified within spatial lag model residuals in the form of low-high and high-low spatial outliers. These spatial outliers are potentially an artefact of border bias (Hendricks and Mark-Carew, 2017) resulting from a shared border with the Lyme disease endemic states of Pennsylvania, Maryland and Virginia (CDC, 2015b). Additionally, they could also be spurious outliers resulting from low or variable county-level veterinary participation, potentially remedied by increasing veterinary practice recruitment efforts in West Virginia north-eastern counties.
The study findings are based upon a state-wide county level ecological regression between human Lyme disease incidence per 100,000 persons and animal sentinel data collected as a convenience sample from 2014-2016. Spatial lag models are often conducted for data aggregated at the county level, but are potentially prone to the ecological fallacy if data incorporated are at differing spatial scales (Anselin, 2002). Additionally, some studies call for the presentation of vector-borne disease surveillance data at scales finer than the county level (Eisen, 2007). We attempted to adjust for the ecological fallacy and misclassification biases by incorporating confirmed cases of human Lyme disease obtained within the patient’s home county of residence and animal sentinel data aggregated at the same county level scale. Our analysis was limited to the county-level scale as this is the finest scale for which travelassociated infection status can be noted. Incorporation of the assumption that all veterinary practices received equal opportunity to participate and would participate uniformly was employed within our analysis to conserve sample size. This is a potentially important study limitation, while other potential limitations are associated with differing reporting biases associated with unequal Lyme disease awareness among physicians or other public health officials and varying presentation of Lyme disease signs and symptoms among patients infected (Hayes and Piesman, 2003; Bacon et al., 2008; Brett et al., 2014; Sanchez et al., 2016; Steere et al., 2016). Unfortunately, we were unable to characterise the effect of differing reporting biases due to the aggregate nature of our retrospective investigation.
Despite limitations stated above, our findings fill significant gaps in the literature regarding the spatial biases of public health surveillance data. Results confirm the efficacy of dog sentinel surveillance programmes in association with human Lyme disease. Interestingly, residual spatial dependence was identified among spatial models and potentially point to other geographic biases not considered. In particular, West Virginia is prone to border biases resulting from shared boundaries with states both high and low incidence of Lyme disease (CDC, 2015b; Hendricks and Mark- Carew, 2017). In addition to border bias, residual spatial dependence is potentially the result of geophysical, vegetation, soil or other ecological variables known to influence tick from previous studies (Guerra et al., 2001; Diuk-Wasser et al., 2012). Future research is warranted to incorporate these ecological variables as well as behavioural variables into spatial models to produce more accurate visualisation of human Lyme disease risk, thus including parameters such as dog ownership in rural vs urban counties where the level of exposure to tick habitats may differ. Results highlight the unique geographic position that West Virginia shares with its neighbouring states, and emphasise the necessity of spatial analysis to adjust for inherent biases associated with collection of public health surveillance data.