Integrated Social-Behavioral and Ecological Risk Maps to Prioritize Local Public Health Responses to Lyme Disease.

BACKGROUND: The risk of contracting Lyme disease (LD) can vary spatially because of spatial heterogeneity in risk factors such as social-behavior and exposure to ecological risk factors. Integrating these risk factors to inform decision-making should therefore increase the effectiveness of mitigation interventions.
OBJECTIVES: The objective of this study was to develop an integrated social-behavioral and ecological risk-mapping approach to identify priority areas for LD interventions.
METHODS: The study was conducted in the Montérégie region of Southern Quebec, Canada, where LD is a newly endemic disease. Spatial variation in LD knowledge, risk perceptions, and behaviors in the population were measured using web survey data collected in 2012. These data were used as a proxy for the social-behavioral component of risk. Tick vector population densities were measured in the environment during field surveillance from 2007 to 2012 to provide an index of the ecological component of risk. Social-behavioral and ecological components of risk were combined with human population density to create integrated risk maps. Map predictions were validated by testing the association between high-risk areas and the current spatial distribution of human LD cases.
RESULTS: Social-behavioral and ecological components of LD risk had markedly different distributions within the study region, suggesting that both factors should be considered for locally adapted interventions. The occurrence of human LD cases in a municipality was positively associated with tick density (p<0.01) but was not significantly associated with social-behavioral risk.
CONCLUSION: This study is an applied demonstration of how integrated social-behavioral and ecological risk maps can be created to assist decision-making. Social survey data are a valuable but underutilized source of information for understanding regional variation in LD exposure, and integrating this information into risk maps provides a novel approach for prioritizing and adapting interventions to the local characteristics of target populations. https://doi.org/10.1289/EHP1943.

Introduction

In Canada and elsewhere in North America, Lyme disease (LD) is caused by the spirochete Borrelia burgdorferi sensu stricto (for simplicity, B. burgdorferi will be used hereafter) and is transmitted by the blacklegged tick Ixodes scapularis in northeastern and midwestern North America (Burgdorfer et al. 1982). LD has three clinical phases: early LD in which most patients develop a characteristic erythema migrans (EM) skin rash along with nonspecific influenza-like symptoms; early disseminated LD with one or more manifestations of multiple EM lesions or neurological or cardiac problems; and late LD with neurological or arthritic manifestations (Wormser et al. 2006). The earlier antibiotic treatment is started, the more rapid and successful is patient recovery and the lower the impact on patients, their families, and the healthcare system (Wormser et al. 2016).

The increasing risk of LD in Canada is driven by the expanding geographic range of the tick vector I. scapularis in eastern and central Canada (Bouchard et al. 2015). The incidence of human LD cases has been increasing in southern parts of the country, with 987 cases reported in 2016 (PHAC 2017) compared to 338 cases in 2012, 682 cases in 2013, 522 cases in 2014, and 917 cases in 2015. In the province of Quebec, LD risk was first confirmed in the Montérégie region in 2008 (Bouchard et al. 2011; Ogden et al. 2008b, 2010), and now human LD cases are more frequently being reported in this region, reaching around 50 cases reported per year since 2013 compared to 13 cases in 2012 (MSSS 2016a).

There is no LD vaccine, and preventive actions by public health organizations have focused on public health campaigns aimed at communicating information about a) symptoms of LD so that the public seeks medical advice when symptoms of infection are observed; and b) measures to prevent LD infection, such as personal protection strategies and control of ticks in the environment (Ogden et al. 2015). LD preventive strategies include tick bite prevention measures such as wearing protective clothing, using insect repellent, and checking for and promptly removing ticks. Tick control interventions consist of the application of acaricides to vegetation and management of tick habitat through landscaping (reviewed in Piesman and Eisen 2008).

The risk of acquiring LD generally correlates positively with the density of host-seeking infected ticks in the environment. This is the product of tick density and the prevalence of infection in the ticks; both of which are spatially heterogeneous at a range of spatial scales (Pepin et al. 2012). This ecological risk for LD is spatially heterogeneous, and the pattern of LD emergence and spread in Canada is shaped by key ecological drivers affecting the geographic range and population dynamics of ticks and the dynamics of B. burgdorferi infection in ticks. The spatial distribution and local dynamics of I. scapularis and B. burgdorferi reflect seasonal and annual dispersal of ticks by animal hosts (Bouchard et al. 2013a, 2013b; Ogden et al. 2008a; Scott et al. 2008, 2012), meteorological and climatic conditions and climate warming that affect tick population survival (Leighton et al. 2012; Ogden et al. 2014), habitat suitability for tick hosts (Bouchard 2013; Bouchard et al. 2013a, 2013b; Lindsay et al. 1998, 1999; Ogden et al. 2006), B. burgdorferi reservoir host abundance (e.g., the white-footed mouse, Peromyscus leucopus) (Bouchard et al. 2011; Rogic et al. 2013; Roy-Dufresne et al. 2013), and other aspects of the host community structure, such as biodiversity, that influence B. burgdorferi transmission cycles (Bouchard et al. 2013a; Ogden and Tsao 2009; Werden et al. 2014).

The extent to which ecological risk factors influence the risk of people contracting LD also depends on factors that increase human exposure to infected ticks, such as location of residence and consequent local or peri-domestic LD risk exposure (Connally et al. 2009; Zeman et al. 2015), people’s occupations, and where they undertake their leisure activities (Quine et al. 2011). Preventive communication can be targeted to known populations with occupations or leisure activities leading to higher risk of exposure to infected ticks (Quine et al. 2011). A significant challenge that remains, however, is effectively communicating risk to the general public living in risk areas in a way that results in adoption of preventive behaviors (Connally et al. 2009; Cromley et al. 1998; Dister et al. 1997; Finch et al. 2014).

Many factors influence the adoption of preventive behaviors toward LD in the general population. Both knowledge and risk perception about LD have been associated with the degree of adoption of individual preventive behaviors in Canada (Aenishaenslin 2014, 2015a), the United States (Brewer et al. 2004; Hallman et al. 1995; Herrington et al. 1997, 2004; Hook et al. 2015; Valente et al. 2015), and Europe (Beaujean et al. 2013a, 2013b). Behaviors, influenced by knowledge and risk perception, are key risk factors that can be changed by appropriate risk communication. Therefore, public knowledge and risk perception of the disease should be considered in the design of regionally targeted prevention and control strategies. However, information about regional variation in this social-behavioral component of risk is rarely available to public health decision-makers, and its interaction with regional ecological risk factors in determining the geographic pattern of LD risk has not been assessed.

Objectives

The objective of this study was to develop a novel risk-mapping method for LD, integrating both social-behavioral and ecological risk factors in order to identify priority areas for LD interventions. To do so, we combined spatially explicit data on LD knowledge, risk perception, and preventive behaviors derived from a web survey conducted in 2012 in the Montérégie region of southern Quebec with annual surveillance data on ticks collected in the environment across the same region from 2007–2012. We then used these combined data to generate integrated LD risk maps and evaluated their ability to predict the spatial pattern of human LD cases at the municipal scale across the region. We use here the climate change adaptation terminology (Cardona et al. 2012) to define the risk of LD in the population investigated as the product of the population exposure to the ecological risk factor and their vulnerability (in this case, the capacity of the population to adopt LD preventive behaviors: the social-behavioral risk factor component).

Methods

Study Area

Montérégie is a region of southwestern Quebec with a population of (http://www.stat.gouv.qc.ca/statistiques/profils/region_00/region_00.htm) and an area of . Two geographic subunits were used to explore, describe, and analyze the data: a) the local health unit (Centre intégré de santé et de services sociaux; CISSS); and b) the census subdivision (CSD), with each CSD corresponding roughly to a single municipality. In 2012, Montérégie region was made up of 179 CSDs grouped in 4 local health units: CISSS Montérégie-Est (comprising the northeastern part of the Montérégie region), CISSS Montérégie-Centre (comprising the central part of the Montérégie region from the head of Lake Champlain and following the Richelieu river valley), CISSS Montérégie-Ouest (comprising the southwestern part of the Montérégie region), and CISSS Estrie (the southeastern part of the Montérégie region, which was subsequently incorporated into the neighboring Estrie region in April 2015).

Social-Behavioral Risk Factors

Data on LD knowledge, risk perception, and preventive behaviors were collected using a cross-sectional survey of 401 respondents conducted in Montérégie in 2012. The respondents were selected randomly from a panel administered by the survey firm Leger (http://leger360.com/en-ca/home). The sample was representative of the entire study region in terms of sociodemographic factors, including gender, age, education, income, and population size within the region. Questions were designed to measure levels of knowledge about LD, risk perception, and adoption of specific individual preventive behaviors. Further details on the survey design and administration, including the full questionnaire, can be found elsewhere (Aenishaenslin et al. 2014, 2015a). The survey was reviewed by the ethics committee for health research of the University of Montreal (Comité d’éthique de la santé, CERES, certificate number 12-050-CERES-D), and informed consent was obtained from all respondents.

Our risk maps were derived from three index scores computed from the survey data at the respondent level using methods described previously (for details on the questionnaire that was used for data collection, the scores, and their distribution, see Aenishaenslin et al. 2014, 2015a). The index scores were: a global preventive behavior score (GPB), a global knowledge score (GK), and a global risk perception score (GRP). The GPB represents a measure of the level of adoption of the three most commonly recommended preventive behaviors (i.e., performing regular tick checks, wearing protective clothing, and using tick repellents). The GPB score is 2 (high) if the respondent often or always applied two of the three preventive behaviors, 1 (moderate) if the respondent often or always applied one of the three behaviors, and 0 (nil) when preventive behaviors were not applied. The GK score is calculated as the sum of four knowledge questions about a) LD transmission, b) early symptoms, c) treatment, and d) the location of risk regions. The score ranges from 0 to 4, corresponding to the number of correct answers. Participants who declared that they had never heard of LD before the survey were automatically given a score of 0. The GRP score is the mean of four risk perception variables measured using Likert scales from 1 to 5 in the questionnaire: perceived severity of LD, perceived individual susceptibility to LD, perceived regional susceptibility to LD, and feeling of worry. The location of residence for each respondent was geocoded to their six-digit postal code to produce data point maps for the global index scores. The postal code conversion file of Statistics Canada was used to determine the centroid of each postal code.

Ecological Risk Factors

We used active field surveillance data collected from 2007 to 2012 to estimate the density of I. scapularis tick populations in the Montérégie region. The observed tick density in field surveys provides a measure of ecological risk related to factors such as host abundance and diversity, climate conditions, and habitat-specific risk factors that determine tick density and distribution. Active field surveillance involves collecting ticks using a standardized drag sampling method (trailing a square of white flannel cloth for at least three person-hours per site). Field surveillance was carried out at 183 different sites from 2007 to 2012, with each site visited at least once during this period. Tick sampling was carried out between May and October, encompassing the period of year when weather conditions are most favorable for tick activity. There was a total of 378 site visits, of which 104 sites were revisited two or more times (56.8% of the sites), and 79 sites were visited only once. Sampling sites were selected among suitable forested habitat in the target CSD based on the following criteria: deciduous (maple or mixed deciduous) woodland of minimal dimensions of , ease of access, and owner authorization. Most of these sites were private woodlots. In 2007–2008, the sites were visited one to three times per year. More details on the selection of sites, sampling method, and inclusion criteria for the active field surveillance and revisits can be found elsewhere (Bouchard et al. 2011, 2015; Ogden et al. 2010). The geolocation of each field site for analyses was the centroid of the sampling plot geolocated using a handheld Global Positioning System device.

The proportion of sites with blacklegged ticks and the density of ticks at these sites increased from 2007 (Bouchard et al. 2011) to 2012 (Ogden et al. 2014). To obtain a standardized tick abundance index for all sites that could be compared to the 2012 survey data, we modeled the increase in tick numbers from 2007–2012 and removed this temporal trend. In this model, the outcome variable was the log count of nymphs and adults combined at each site visit, and an identifier number for each site was included as a random effect to account for clustering of data within a site. Year of sampling was included as a categorical fixed effect. Month of sampling was not included in the model because the outcome variable did not distinguish between questing nymphs and questing adults, which are mostly active in spring months for nymphs and in the end of summer to fall months for adults in our study region. Sampling effort (area covered by drag sampling at each site visit) was included as an offset in our model. We produced a data point map based on each site location and the value of predicted tick density (PTD) per site as the measure used to interpolate the study area tick density as described below.

Human Population at Risk

The density of the human population at risk in each CSD was calculated based on the area of the CSD and the human population from the 2011 census carried out by Statistics Canada (http://www12.statcan.gc.ca/census-recensement/index-eng.cfm). The human population density measure was geolocated to the CSD centroid for interpolating population density for over the study area as described below.

Geoprocessing, Standardization, and Construction of Spatial Data Layers for Social-Behavioral and Ecological Risk Factors

In order to objectively compare and integrate data layers, all variables were first standardized by subtracting the mean and dividing by the standard deviation (SD) to obtain and (Quinn and Keough 2002). This resulted in a common scale without units for each variable. All mapping and geospatial methods were performed in ArcGIS (version 10.2; Esri) using the Statistics Canada Lambert projection and the North American Datum (NAD) 1983 datum.

We produced a standardized interpolated surface for a) the three social index scores (i.e., GPB, GK, GRP) given the resident locations of survey participants; b) the density of the population at risk (per ) given the CSD centroids; and c) the predicted value of the density of ticks given the locations of active surveillance sites (PTD). GPB was used to represent the social-behavioral risk of the population in analyses because we assume that overall preventive behaviors are a consequence of each individual’s level of knowledge (represented by GK) and risk perception (represented by GRP). Maps of GK and GRP are presented for comparison purposes.

Interpolation has been used in other studies of vector-borne diseases that predict surfaces for ecological variables (Bhunia et al. 2013; Epp et al. 2011; Healy et al. 2014). The ability of interpolated surfaces to represent the underlying variable in the study area depends on the assumption that the variable is spatially correlated. No significant global spatial autocorrelation was found, but significant local spatial autocorrelation was found for each social-behavioral risk factor (Figures S1 to S3). Interpolation was derived using the inverse distance weighting (IDW) method, which assigns a greater weight to nearby points than to distant points (resolution of , power of 2, , using a variable search radius, low-pass filter, and cubic convolution display). We set the search neighborhood for the IDW at a fixed number of points (i.e., 12) because given the uneven spatial distribution of respondents, a fixed-distance approach would have led to too many respondents driving the algorithm in densely surveyed urban areas and too few respondents driving the algorithm in sparsely surveyed rural areas. Furthermore, the mean search radius for 12-respondent neighborhoods given all respondents in the study area approximated the distance of the detected local spatial autocorrelation in the respondent data. To validate our IDW interpolation parameterizations, we used cross-validation statistics to confirm that training and validation datasets produced interpolated maps that varied within an acceptable range of the underlying variable (Table S1).

Interpolated surfaces were summed to create two types of integrated risk maps: a social/behavioral-ecological vulnerability index map that sums the negative product of GPB and the PTD raster map (i.e., ), and a prioritization index map that sums the social/behavioral-ecological vulnerability index map and the population-at-risk map (see Figure 1).

Figure 1.

Conceptual structure for the integration of social-behavioral and ecological risk to adapt local responses to Lyme disease risk.

Association between Social-Behavioral and Ecological Risk Factors with the Presence of Human Cases of Lyme Disease

We used logistic regression models to determine whether the estimated social/behavioral-ecological vulnerability index for each CSD and the estimated social/behavioral-ecological prioritization index for each CSD was a significant predictor of the presence/absence of LD cases at the CSD level. In addition, we modeled associations between LD cases and the average social-behavioral scores, PTD, and the human population density for each CSD. This was done to determine which, if any, underlying variables may be driving a significant relationship between LD presence/absence and the social/behavioral-ecological vulnerability or prioritization index.

We used available epidemiological data on confirmed and probable human cases of LD [according to the case definition reported in Montérégie at the CSD level in 2012 and 2013 (MSSS 2016b)]. A total of 69 locally acquired cases of LD were reported from 37 of 179 CSD within Montérégie in 2012 and 2013, and cases with known travel history and possible acquisition of LD outside of Quebec were not included in the analysis. In both model formulations, the response variable of human case presence, , occurred for one or more LD cases in a given CSD, while for no reported cases. In model formulation 1, we modeled the presence/absence of LD cases in each CSD as a function of individual predictors, including PTD, human population density, and either GK and GRP or GPB only (since it is a function of both GK and GRP) as independent variables (models 1A, 1B, and 1C). In model formulation 2, we modeled the social/behavioral-ecological vulnerability index for each CSD (model 2A) or the social/behavioral-ecological prioritization index for each CSD (model 2B) in two separated models as the only predictors of the presence/absence of LD cases. In each model, explanatory variables were tested individually with a liberal cutoff of in univariable model, and then we selected the most parsimonious multivariable model through a process of forward and backward substitution and elimination. The cutoff for keeping a variable in the final model was . Internal validity of the final models was verified by visually assessing the linearity of the predictors with the log odds of the outcome and the pattern of residuals. The area under the receiver operating characteristic curve was used as a measure of the predictive ability of the model. The presence of residual spatial autocorrelation was explored in posthoc analyses of residuals plotted as a function of geographic distance for both model formulations (Dormann et al. 2007). Finally, we used Pearson’s chi-square test to detect significant differences in the proportion of respondents with higher social-behavioral index scores at the CISSS level. Statistical analyses were carried out using R (version 3.2.4; R Development Core Team).

Results

Descriptive analysis of human LD cases (number and prevalence), human population density, social-behavioral risk measures (proportions of respondents with high GPB, GK, and GRP), and ecological risk (PTD) at the CISSS level are presented in Table 1. The geographic distribution of social-behavioral and ecological risk data, point maps of the participants recruited in the survey in 2012, and the field surveillance site locations visited in Montérégie from 2007–2012 are illustrated in Figures 2 and 3. The distribution of tick sampling sites provided good spatial coverage of the study area for the population survey (Figures 2 and 3). Social-behavioral risk, human population density, and ecological risk were all spatially heterogeneous (Figures 4–7 and Table 1). In all maps, darker areas represented higher intensity of the variable shown; however, the interpretation of the results in terms of risk is the reverse for the social-behavioral risk factors because a higher level of social-behavioral risk corresponds to lower risk for LD.

Table 1

Number of human Lyme disease (LD) cases in 2012–2013 and the proportion of respondents with higher social-behavioral index scores (lower LD risk) and higher ecological risk scores by health unit (Centre integré de santé et de services sociaux; CISSS) in the Montérégie region.

CISSS	No. of human LD cases				Human populationa (n)	Area (km^2)	Human population densitya (no./km^2)	Social/behavioral or ecological risk: proportion (percentage) of respondents with highest index scores (numerator) given the total number of surveyed participants at CISSS level (denominator)				Ecological risk (normalized value)	No. of field site locations
	2012	2013	2012–2013	No./100,000				GPBb	GKc	GRPd	PTDe	PTDf
Estrie	6	34	40	27	147,284	2,572	57	9/20 (45.0%)	7/36 (19.4%)	11/36 (30.6%)	7/36 (19.4%)	High (0.3)	46
Montérégie Ouest	5	9	14	3	410,921	3,748	10	11/48 (22.9%)	10/107 (9.3%)	38/107 (35.5%)	2/107 (1.9%)	Low (−0.3)	49
Montérégie Est	2	8	10	3	372,299	3,519	106	15/70 (21.4%)	25/131 (19.1%)	37/131 (28.2%)	46/131 (35.1%)	High (0.1)	46
Montérégie Centre	0	5	5	1	509,187	1,443	353	7/43 (16.2%)	13/82 (15.8%)	21/82 (25.6%)	35/82 (42.7%)	High (0.3)	42

Note: GPB, global preventive behavior index; GK, global knowledge index; GRP, global risk perception index; PTD, predicted tick density index.

Human population size and density data per CISSS were derived from the human population from the 2011 census carried out by Statistics Canada (http://www12.statcan.gc.ca/census-recensement/index-eng.cfm).

Respondents with values of 2 (numerator).

Respondents with values of 3 or 4 (numerator).

Respondents with values at 75e centile to 100e centile (numerator).

Respondents from a region with a PTD normalized values from 0 to 1 were classified as higher risk (numerator).

Figure 2.

Point map showing the distribution of locations of residence of survey participants in Montérégie in 2012 ().

Figure 3.

Point map showing the distribution of the field surveillance site locations visited in Montérégie, 2007–2012 ( site visits).

Spatial Variation in Social-Behavioral and Ecological Risk

High predicted LD social-behavioral risk, corresponding to lower adoption of preventive behaviors, was very heterogeneous, but was particularly low within the CISSS Montérégie-Est and CISSS Montérégie-Centre (Figure 4A and Table 1). The levels of knowledge and risk perception were also heterogeneous across all CISSSs (Figures 4B and 4C and Table 1). The lowest level of knowledge was observed in CISSS Montérégie-Ouest, and the lowest level of risk perception was observed in CISSS Montérégie-Centre (Table 1). The highest human population densities occur in northwestern part of the CISSS Montérégie-Centre and the southwestern part of the CISSS Montérégie-Est.

Figure 4.

Spatial variation in the level of (A) adoption of individual preventive behaviors (i.e., global preventive behavior score; GPB); (B) knowledge about the disease (i.e., global knowledge score; GK); and (C) risk perceptions (i.e., global risk perception score; GRP). Darker areas represent higher index scores for these social-behavioral drivers but a lower risk for Lyme disease risk transmission.

The highest tick population densities were estimated to occur in CISSS Montérégie-Centre, CISSS Estrie, and the southern part of CISSS Montérégie-Est (Figure 5 and Table 1). Higher social/behavioral-ecological vulnerability index values (Figure 6) and prioritization index values (Figure 7) were also found in CISSS Montérégie-Centre, CISSS Estrie, and the southern part of CISSS Montérégie-Est, although the spatial extents were different given that the prioritization index accounted for the density of the human population at risk.

Figure 5.

Spatial variation in the level of (A) population density based on the 2011 census, and (B) predicted tick density (PTD) based on ticks collected through active field surveillance from 2007–2012. Darker areas represent a higher Lyme disease risk.

Figure 6.

Social/behavioral-ecological vulnerability index map. Darker areas represent a higher Lyme disease risk based on the global preventive behavior score and the predicted tick density (PTD).

Figure 7.

Prioritization index risk map. Darker areas represent a higher Lyme disease risk based on the social/behavioral-ecological vulnerability index and the human population at risk.

Association between Risks Factors and the Presence of Human Lyme Disease Cases

Of the individual socio-behavioral risk scores, PTD, and human population density, CSD-level PTD was the only statistically significant predictor of CSD-level LD cases [Table 2, model 1C, ; 95% confidence interval (CI): 1.28, 2.52; and ]. The combined social/behavioral-ecological vulnerability index was also a significant predictor of LD cases (Table 2, model 2A: ; 95% CI: 1.25, 2.53; and ), but the social/behavioral-ecological prioritization index (that also accounts for human population density) was not (Table 2, model 2B: ; 95% CI: 7.85, 1.52; and ).

Table 2

Odds ratios (95% CI) and p-values for the presence/absence of LD cases in a census subdivision (CSD) as a function of individual predictors, including social-behavioral risk factors (either GK and GRP or GPB only), predicted tick density, and population density (model 1A–1C) or as a function of social-behavioral vulnerability index (model 2A) and prioritization index (model 2B).

Model formulation	Model 1A		Model 1B		Model 1C		Model 2A		Model 2B
CSD-level predictor	OR (95% CI)	p-value	OR (95% CI)	p-value	OR (95% CI)	p-value	OR (95% CI)	p-value	OR (95% CI)	p-value
Global knowledge score (GK)	0.87 (0.57, 1.29)	0.52	–	–	–	–	–	–	–	–
Global risk perception score (GRP)	0.92 (0.63, 1.35)	0.68	–	–	–	–	–	–	–	–
Global preventive behavior score (GPB)	–	–	087 (0.59, 1.26)	0.46	–	–	–	–	–	–
Predicted tick density (PTD)	1.75 (1.25, 2.51)	0.002	1.76 (1.25, 2.54)	0.002	1.77 (1.28, 2.52)	0.0009	–	–	–	–
Population density	0.66 (0.34, 1.06)	0.15	0.66 (0.34, 1.08)	0.16	–	–	–	–	–	–
Social/behavioral-ecological vulnerability index	–	–	–	–	–	–	1.76 (1.25, 2.53)	0.002	–	–
Social/behavioral-ecological prioritization index	–	–	–	–	–	–	–	–	1.08 (0.75, 1.52)	0.65

Note: Area under the curve (AUC) values were 0.70, 0.70, 0.67, 0.62, and 0.55 for models 1A, 1B, 1C, 2A, and 2B, respectively. –, is if a predictor was not included in the model formulation x; CI, confidence interval; OR, odds ratio.

The number of households surveyed ranged from 0–15 per CSD, with a total of 99 CSD out of 179 CSD having at least one respondent. The number of field site visits per CSD ranged from 0 to 11 site visits, with a total of 154 CSD out of 179 CSD having at least one site-visit (Table 1). No significant differences in the proportions of respondents with higher social-behavioral index scores were detected among CISSSs.

Significant local spatial autocorrelation was found for each social-behavioral risk index at the survey respondent level, as illustrated by spatial clusters and spatial outliers maps (Figures S1 to S3), although no significant global spatial autocorrelation was detected in any of the social-behavioral variables. We detected significant positive and negative clusters or outliers as explained in supplemental materials. Our logistic regression modeling approach adequately accounted for the spatial autocorrelation given that this process was not detected from an assessment of the residuals using the Moran’s I method (, , and for both models).

The cross-validation statistics for assessing our interpolated surfaces are presented in the supplemental materials (Table S1).

Discussion

To our knowledge, this study is the first to develop integrated social/behavioral-ecological risk maps that can help target LD interventions. Our findings illustrate heterogeneities in social-behavioral and ecological risks within a newly endemic region for LD that are significantly correlated with the spatial distribution of emerging human cases. Significant regional differences in levels of LD knowledge and preventive behaviors have been recently documented in Canada (Aenishaenslin et al. 2016, 2017), but our integration of this social-behavioral risk variation with ecological risk allowed the combined effects of these factors on risk to be directly assessed for the first time and mapped in relation to emerging hotspots for human LD cases. The identification of geographic hotspots of LD risk using this approach provides both useful targets for geographic prioritization of public health interventions to address the areas and populations at greatest risk, but also an assessment tool to examine of the relative contribution of different types of risk factors across the emerging landscape of LD risk.

In recent years, a number of new approaches have been applied to help prioritize LD interventions, including decision-theoretic approaches such as multicriteria decision analysis (Aenishaenslin et al. 2015b; Hongoh et al. 2011). Although such approaches allow the inclusion of multiple factors in the prioritization process, such models tend not to take into account the spatial characteristics and differences in LD risk across the targeted region. The current study provides first steps in this direction by proposing integrated social/behavioral-ecological risk maps that can be used in different ways to inform public health authorities. First, localities showing the highest levels of LD social/behavioral-ecological risk can be prioritized for the implementation of prevention or control strategies in a context of limited resources. For example, based on our results, Longueuil CSD, located in the southwestern part of the CISSS Montérégie-Est, as well as other CSDs located at the southern part of the CISSS Montérégie-Centre near Lake Champlain and the border with the United States, would be considered with a higher priority given the low preventive behavior score, higher population density, and higher predicted tick population density, which results in high social/behavioral-ecological and prioritization index values (Figures 4–7 and Table 1). Secondly, the social-behavioral and ecological risk maps can be interpreted individually to inform public health decision-makers about which interventions to target in each CSD. For example, if a CSD shows moderate levels of ecological risk but low levels of knowledge and preventive behaviors, then resources should be invested in risk communication to public and health professionals.

Previous studies have found significant spatial associations between ecological factors linked to tick population abundance and LD risk (Eisen et al. 2006; Pepin et al. 2012; Eisen and Eisen 2016). We also found a statistically significant association between tick density and LD case presence; however, we found that LD case presence was not significantly associated with any of the social-behavioral variables. We may not have had a sufficient number of respondents at the CSD level to adequately represent the social-behavior risk factors of the local population. Secondly, the administrative construct of CSD boundaries may not be an appropriate spatial scale or geographically positioned in a manner that captures the processes driving variation in social behavior and LD exposure and risk leading to LD reports. The impacts of this ecological bias and the related concept of the modifiable areal unit problem on epidemiological analyses were not explored in the current study. There may also be limitations relating to the web survey design that prevent adequate characterization of the population-level social-behavior (Aenishaenslin et al. 2014) and the fact that the social data were collected for other purposes. Furthermore, the uneven distribution in respondents (i.e., more respondents around Montreal area, Figure 2) can lead to error in modeling and interpretation when scaling up data from low-density areas compared to high-density areas if the survey data do not adequately capture the true variation across the study area. There can also be an issue of overinterpreting the interpolated maps. This method assumes that spatial variation is a consequence of the spatial distribution of respondent locations and does not take into account the potential for variation between their locations.

An alternative method would be to create regression models where LD cases are explained by social-behavioral, ecological (tick density), and human population (density) variables and then used the models to predict LD case presence given a denser and more homogeneous distribution of covariate values. A similar method was developed for the present analyses to produce the ecological risk index (PTD) for each site location, which was then interpolated and visually validated with unpublished and published maps (Unpublished data, 2018; Leighton et al. 2012). However, it was not possible to apply this approach for the social-behavioral risk data in this study given that collection and analysis of survey data for LD in the study region was only carried out in a single year using a cross-sectional design that does not allow validation of multiannual regression models. Moreover, given that there is evidence of rapid changes in LD awareness across Canada based on two national surveys conducted in 2014 (Aenishaenslin et al. 2016), the development of social-behavioral risk maps should be repeated and adapted over time.

LD is still considered as a newly endemic disease in southern Québec, and the levels of knowledge and preventive behaviors in the Montérégie region are low overall despite some geographical heterogeneity (under 20% of the population had a high level of knowledge, Table 1). Our results suggest that in our study area, ecological risk is a major driver of LD human cases, an observation that fits with its emerging epidemiological status, but that could be expected to change as epidemiological context evolves. CISSS Estrie had the highest number of human cases per capita, which could be explained by the high level of ecological risk, even though relatively high levels of adoption of preventive behavior and risk perception were observed. Further studies comparing the relative importance of social-behavioral and ecological risks in emerging vs. highly endemic regions would be needed to confirm these patterns and better understand the role of knowledge, risk perception, and preventive behaviors in modulating the impact of emerging ecological risk on regional case incidence.

We caution not to overinterpret our results from using only one level of adopting preventive behaviors to estimate the level of social-behavioral risk for LD. The types of social risk factors that may contribute to increase the risk of contracting LD are multiple. For example, risk behaviors, culture, and values all are known to have impacts on health issues and could also be geographically heterogeneous (Rosenstock et al. 1988), but were not measured in the survey. For LD, very few studies have addressed these factors and have tried to quantify their relative importance. Data collection and data access to these different layers of information can be a limit. Surveillance of human behavior to a particular infectious disease is not easily found, and this could limit and decrease the speed of future work investigating those social-behavioral risks and their influence on the risk of infection. Further research on the social determinants of this disease is needed in order to improve this work and methodology.

Conclusions

This current study is a first step towards an integrated approach for LD risk assessment, identifying social-behavioral risk factors that act jointly with ecological risk factors to influence management of this emerging tick-borne disease. This study highlights a new area for research on the development and validation of integrated social/behavioral-ecological risk maps. The approach developed here is widely applicable to other vector-borne or zoonotic diseases and to different epidemiological contexts where the geographic patterns of risk are driven by the interplay between social-behavioral and ecological factors.

Click here for additional data file.