Analyzing the Potential Risk of Climate Change on Lyme Disease in Eastern Ontario , Canada Using Time Series Remotely Sensed Temperature Data and Tick Population Modelling

The number of Lyme disease cases (Lyme borreliosis) in Ontario, Canada has increased over the last decade, and that figure is projected to continue to increase. The northern limit of Lyme disease cases has also been progressing northward from the northeastern United States into southeastern Ontario. Several factors such as climate change, changes in host abundance, host and vector migration, or possibly a combination of these factors likely contribute to the emergence of Lyme disease cases in eastern Ontario. This study first determined areas of warming using time series remotely sensed temperature data within Ontario, then analyzed possible spatial-temporal changes in Lyme disease risk in eastern Ontario from 2000 to 2013 due to climate change using tick population modeling. The outputs of the model were validated by using tick surveillance data from 2002 to 2012. Our results indicated areas in Ontario where Lyme disease risk changed from unsustainable to sustainable for sustaining Ixodes scapularis (black-legged tick) populations. This study provides evidence that climate change has facilitated the northward expansion of black-legged tick populations’ geographic range over the past decade. The results demonstrate that remote sensing data can be used to increase the spatial detail for Lyme disease risk mapping and provide risk maps for better awareness of possible Lyme disease cases. Further studies are required to determine the contribution of host migration and abundance on changes in eastern Ontario’s Lyme disease risk.


Introduction
The impact of climate change on vector-borne disease systems in North America has been subject to much debate and speculation over the last decade [1].Highly sensitive to changes in environmental conditions and patterns, vectors of tick-borne diseases may be particularly vulnerable to climatic changes, which could play a role in the vector's overall development, abundance and activity.Over recent years, the northern distribution limit of tick vectors of the causative agent of Lyme disease, Borrelia burgdorferi (B.burgdorferi), has advanced from the northeastern United States to southeastern Ontario [2,3].This has led to a significant increase in Lyme disease cases throughout Ontario, confirmed through a combination of thorough active and passive surveillance programs [4,5].In 2004, there were only 40 cases of Lyme disease in Canada [4], while in 2016 there were 343 cases in Ontario alone [6].
Ixodes scapularis (I.scapularis), commonly known as the black-legged tick, is the primary vector of B. burgdorferi in Canada.Populations of this tick have been identified in several provinces across Canada including Ontario, Manitoba, Nova Scotia, New Brunswick and Quebec [7].Depending on climatic conditions, the lifespan of the black-legged tick is approximately two to two-and-a-half years and includes three blood-feeding instars: larva, nymph and adult.While many previous studies have focused on analyzing the impact of climate change on I. scapularis ticks in the United States, the extent and rate at which tick populations are colonizing new regions in Ontario is only recently being studied [8].It is difficult to predict the impact of climate change on tick-borne diseases due to the high complexity and dynamic nature of the tick-borne pathogen systems [9].
A combination of several factors may influence the establishment and abundance of I. scapularis populations in southeastern Ontario.Firstly, the greatest densities of endemic ticks are located within close proximity to the northeastern US border, due to the relatively high population density in this region [10].Secondly, large mammalian hosts such as white-tailed deer and migratory birds carry ticks across a range of distances from the USA into Canada every year as ticks are unable to move great distances on their own [11].Thirdly, warming temperatures due to climate change increase the likelihood of immigrating ticks establishing endemic populations in their new territory [3].
The effect of climate on tick population survival is thought to be the effect of temperature on rates of interstadial development and thus on the duration of the lifecycle; if temperature conditions over the multi-year lifecycle (often expressed in terms of cumulative degree days) fall below a threshold, the probability that a larva survives to become a mated adult falls below unity [12].While both host abundance and climate are important factors in the survival of black-legged tick populations, 98% of the black-legged tick's two-year life cycle exists off the host where environmental conditions are particularly important to the survival and development of the tick [2,13].As such, some studies have suggested that climatic differences may be more limiting than host abundance in the spatial distribution of ticks [14].Therefore, an understanding of the influence of climate change on I. scapularis ticks is critical to predicting the geographical distribution of endemic tick populations and therefore to assessing future Lyme disease risk in Ontario.
Multiple studies have suggested that the warming of temperatures as a result of climate change may impact the abundance and survival rates of ticks by opening up previously uninhabitable territory for vectors, increasing their reproductive rates [15,16].Because the tick spends most of their lifetime at the ground level, land surface temperature (rather than air temperature) is a better predictor of survival [17].The temperature measure of environmental suitability for tick populations is typically expressed in terms of accumulated degree days >0 • C [18].In 2005, a climate suitability model of I. scapularis was used to obtain projections for possible spread of the geographic range of the tick in North America.Their results revealed that sustainable habitat regions for black-legged ticks could increase 213% by 2080 based on projected future climate [2].Another study using a climate model predicted that the theoretical range of I. scapularis establishment would move northwards by approximately 200 km by 2020 and 1000 km by 2080 from the current range limit of the tick [3].
A clear distinction between weather and climate should be made when analyzing the impact of climate change on vector-borne disease systems.Weather is most commonly measured as temperature, precipitation, cloudiness, brightness, visibility, wind, and atmospheric pressure over a particular area during a short period of time.Climate is the average measure of weather conditions over long periods of time.The World Meteorological Organization's standard for climate norms are measured by averaging 30 years of measurement for a climate element [19].Climate models are useful in predicting the future impact of these long-term trends.Climate models use existing knowledge of physical relationships between observed environmental variables and unknown variables.Observed values are then fed in to these models to produce the unknown variables as outputs.In this context, observed variables are those measured by researchers such as current and past climatic conditions.Unknown variables are those that predicted based on observed variables, such as future climatic conditions.
In order to most effectively utilize these predictive models, large-scale, high-resolution climate data are needed.Weather stations provide spatially discrete data for fixed points.The limited number of weather stations, particularly in rural areas, means that spatial interpolation is required to create a continuous spatial coverage of weather data.Since temperature varies greatly both spatially and temporally, weather station data can be inadequate for epidemiological studies of climate change [20].Remote sensing can more accurately measure spatial variation in temperature particularly at the regional and local scale where interpolation may introduce error, particularly in regions with wide variations in elevation.Advancements in remote sensing technologies have made it possible to retrieve large-scale climatological data at high resolution, and have been previously used for studying vector-borne diseases [21,22].
One such example of the collection of remote sensing data is the Moderate Resolution Imaging Spectroradiometer (MODIS), an instrument on board the Terra and Aqua satellites, operated by NASA [23].The MODIS instruments capture data on 36 spectral bands that range in both wavelength and spatial resolution.MODIS is particularly useful because it provides global coverage with accurate calibration in thermal infrared bands that were designed to measure land surface temperature [24].MODIS provides land surface temperature and emissivity products, created by a physics-based land surface temperature algorithm from pairs of daytime and nighttime observations in seven MODIS thermal infrared bands (bands 20, 22, 23, 29, and 31-33) from both Terra and Aqua MODIS data [24].The temperature measurements of MODIs are accurate to within 1 • C when measuring temperatures are between −10 and 50 • C [25].
The objective of this study is to use temperature data to determine (a) how climate has changed in Ontario between 2000 and 2013; and (b) whether any observed climate change could be associated with changes to the spatial distribution of Lyme disease risk in eastern Ontario and thus a possible cause for the increased number of Lyme disease cases [4,5].

Study Area, Data and Preprocessing
Southeastern Ontario and the surrounding region was the area of interest for this study.Figure 1 outlines the boundaries of the study area.Temperature data for this study region was obtained from the MODIS product MOD11C3.This is a monthly mean average of land surface temperatures, with a spatial resolution of 0.05 decimal degrees, or 5.6 square kilometers in the region of study.This product was created by taking the daily land surface temperature and emissivity product MOD11B1, reprojecting and averaging the values to 6 km equal-area grids in the sinusoidal projection, and then averaging the values over a month [19].The mean temperature for each month from 1979 to 2013 was extracted for each pixel within the study area in ArcMap (10.3).
Public Health Ontario provided tick surveillance data for the province of Ontario from 2002 to 2012.I. scapularis removed from both humans and animals were submitted between 2002 and 2008.From 2009 to 2012, only ticks attached to humans and removed by a medical professional were submitted to Public Health Ontario.In 2009, Lyme disease became a reportable disease in Ontario, meaning that medical practitioners were obligated to report all cases or suspected cases to their local public health unit [26].The date is reported as the month and year of tick collection, and the location is reported as the location of the medical practitioner who removed the tick.For this study, an assumption was made that the location where the person was bit by a tick is the same as the submitter city.This assumes that a person bitten by a tick will travel to the closest medical office to have it removed.This is a potential source of uncertainty, as in some cases individuals may prefer their family doctor regardless of their relative proximity to the location where transmission occurred.After data cleaning (sampling bias, data with the wrong location), there were 184 unique submitter cities from the Public Health Ontario data where matching locations were found.The Canadian Place Names dataset of the Atlas of Canada [27] was used to locate the submitter cities.Using this dataset, 175 of the 184 unique submitter city entries were found.Locations that could not be found in the place names dataset were found using the Natural Resources Canada query online by geographical name, which returned the approximate latitude and longitude coordinates of the submitter city.Astra was removed from the dataset because it could not be found in either the place names dataset or the Natural Resources website.
Names were identified as points, which is not representative of the geographic area covered by that location's name.Therefore, all of the point locations and names were used to find the corresponding geographic area.The place names were then used to match with three main data sources that contain boundaries outlining areas.The first was Statistics Canada's census shapefiles, which delineate the census subdivisions, consolidated census subdivisions, designated places, and population centers.The remaining places were then first matched with the federal polling station dataset.Any place that could not be matched with the federal polling stations was then matched using Canada Post's forward sortation areas and local delivery units.

Analysis
The analysis includes two parts, both conducted in ArcMap (10.3).The first is to determine whether Ontario's climate has changed.If Ontario's climate has not changed, then it cannot be a possible factor in causing changes in Lyme disease risk.However, if it has changed, then it is a possible factor for causing changes in Lyme disease risk.The second part of the analysis focuses on finding areas of changing Lyme disease risk when the only input variables used for the population model was mean monthly temperature for every grid point within the study area.

Climate Data Analysis
The annual degree days above 0 °C is calculated by cumulatively adding the temperature for each day of the year when the temperature is above 0 °C.Although this study was focused on finding After data cleaning (sampling bias, data with the wrong location), there were 184 unique submitter cities from the Public Health Ontario data where matching locations were found.The Canadian Place Names dataset of the Atlas of Canada [27] was used to locate the submitter cities.Using this dataset, 175 of the 184 unique submitter city entries were found.Locations that could not be found in the place names dataset were found using the Natural Resources Canada query online by geographical name, which returned the approximate latitude and longitude coordinates of the submitter city.Astra was removed from the dataset because it could not be found in either the place names dataset or the Natural Resources website.
Names were identified as points, which is not representative of the geographic area covered by that location's name.Therefore, all of the point locations and names were used to find the corresponding geographic area.The place names were then used to match with three main data sources that contain boundaries outlining areas.The first was Statistics Canada's census shapefiles, which delineate the census subdivisions, consolidated census subdivisions, designated places, and population centers.The remaining places were then first matched with the federal polling station dataset.Any place that could not be matched with the federal polling stations was then matched using Canada Post's forward sortation areas and local delivery units.

Analysis
The analysis includes two parts, both conducted in ArcMap (10.3).The first is to determine whether Ontario's climate has changed.If Ontario's climate has not changed, then it cannot be a possible factor in causing changes in Lyme disease risk.However, if it has changed, then it is a possible factor for causing changes in Lyme disease risk.The second part of the analysis focuses on finding areas of changing Lyme disease risk when the only input variables used for the population model was mean monthly temperature for every grid point within the study area.

Climate Data Analysis
The annual degree days above 0 • C is calculated by cumulatively adding the temperature for each day of the year when the temperature is above 0 • C.Although this study was focused on finding climate changes during the period 2000 to 2013, long-term climatological trends are measured by 30-year trends to encompass annual variability.Therefore, the annual degree day above 0 • C was calculated from 1979 to 2013.Although this encompasses climate change beyond the study's target dates (2000 to 2013), a minimum of 30 years still shows the direction of change over time.
For each pixel, the linear trend of its time series annual degree days was calculated using the least squares fit to the data.The slope was then used to determine whether the location was warming (a slope greater than 0), cooling (a slope less than 0), or unchanged (a slope equal to 0).The slope at each grid location was then mapped to create a raster covering the study area that showed areas of warming, cooling, and no change.

Tick Population Modeling
A periodic ordinary differential equation model developed by Wu et al. simulates the development of ticks by using twelve differential equations to describe the proportion of ticks growing from one developmental stage to the next [28].Each equation of the system represents the change in number over time of a specific stage of the whole tick life cycle (egg-laying adult females, eggs, hardening larvae, questing larvae, feeding larvae, engorged larvae, questing nymphs, feeding nymphs, engorged nymphs, questing adults, feeding adult females and engorged adult females) as shown in Equation ( 2).Each d i (t) (i = 1, 2,...,12) has a temperature-dependent component and a non-temperature-dependent component, both of which are derived from field data [28].The latter is life stage-specific and includes factors such as mortality, number of hosts, and rate of attachment to host.In each equation, these are summarized into a single constant.For instance, represents the rate at which eggs hatch into larvae.T represents temperature in degrees Celsius.The remaining d i values follow the same general format.For this study, temperature data was extracted from MODIS product MOD11C3.All other parameters were obtained from [28].
The role of this model is to calculate the basic reproductive number (R 0 ) for tick populations under changing environmental conditions, particularly in southern Canada.In mathematical biology and population ecology, R 0 represents the number of female offspring produced by a female individual in its whole life under ideal conditions.Due to the seasonality and population structure implicit in the model developed by Wu and colleagues [28], R 0 calculations were conducted based on a methodology developed for time-periodic systems [29,30].
Using a matrix algebra-driven approach, these twelve differential equations were consolidated into the equation: where F(t) and V(t) are 12 × 12 matrices, representing ticks born into the population and the progression of ticks from each developmental stage to the next, respectively.The solution to (3) is linked to the population growth or decay rate by a linear operator referred to in epidemiology as the next-generation operator [28][29][30].This value indicates the average number of offspring surviving until the next generation per tick, and how these offspring are distributed across the population structure.R 0 is calculated as the dominant eigenvalue of a next-generation operator defined from the linear system (3).The value of R 0 indicates the average number of offspring surviving until the next generation produced per tick, which is linked to population growth if R 0 > 1 or decay if R 0 < 1. Full calculations are available in [28].All R 0 calculations were performed in MATLAB.

Model Validation
Model validation was done by comparing tick validation data against the model output for the approximate location of the tick bite to determine if the R 0 value at that location was greater or equal to one within the last two years.This cut-off was chosen because black-legged ticks have a lifespan of approximately two years.An R 0 value greater or equal to one within the last two years prior to the tick bite shows the climate could support a reproductive tick population.Using the spatial boundaries for the locations where ticks were submitted, the corresponding R 0 value was extracted.
The total number of tick submissions from 2002 to 2012 was normalized using the road length in 2014.The purpose of normalizing using the road length was to account for differences in population accessibility.That is, areas with higher road densities may have more tick submissions due to greater chances of tick bites in areas with higher population accessibility, while areas with low road length may have lower chances of tick bites.The road length in each area was used as a proxy for population accessibility.The number of submitted ticks was divided by total road length.The normalized tick validation data was plotted with the model output data in order to visualize the relationship between the two.

Changing Climate
The vast majority of Ontario within the study area exhibited warming trends between 1979 and 2013, although the very southernmost tip of Ontario has experienced a modest cooling trend (Figure 2).Overall, the climate within most of Ontario has been changing and become warmer, with slightly greater overall temperature increases in northern Ontario.
Remote Sens. 2017, 9, 608 6 of 12 linked to the population growth or decay rate by a linear operator referred to in epidemiology as the next-generation operator [28][29][30].This value indicates the average number of offspring surviving until the next generation per tick, and how these offspring are distributed across the population structure.R0 is calculated as the dominant eigenvalue of a next-generation operator defined from the linear system (3).The value of R0indicates the average number of offspring surviving until the next generation produced per tick, which is linked to population growth if R0 > 1 or decay if R0 < 1. Full calculations are available in [28].All R0 calculations were performed in MATLAB.

Model Validation
Model validation was done by comparing tick validation data against the model output for the approximate location of the tick bite to determine if the R0 value at that location was greater or equal to one within the last two years.This cut-off was chosen because black-legged ticks have a lifespan of approximately two years.An R0 value greater or equal to one within the last two years prior to the tick bite shows the climate could support a reproductive tick population.Using the spatial boundaries for the locations where ticks were submitted, the corresponding R0 value was extracted.
The total number of tick submissions from 2002 to 2012 was normalized using the road length in 2014.The purpose of normalizing using the road length was to account for differences in population accessibility.That is, areas with higher road densities may have more tick submissions due to greater chances of tick bites in areas with higher population accessibility, while areas with low road length may have lower chances of tick bites.The road length in each area was used as a proxy for population accessibility.The number of submitted ticks was divided by total road length.The normalized tick validation data was plotted with the model output data in order to visualize the relationship between the two.

Changing Climate
The vast majority of Ontario within the study area exhibited warming trends between 1979 and 2013, although the very southernmost tip of Ontario has experienced a modest cooling trend (Figure 2).Overall, the climate within most of Ontario has been changing and become warmer, with slightly greater overall temperature increases in northern Ontario.

Lyme Disease Risk
The range of potential sustainable tick habitat was found to have expanded (Figure 3).The R 0 values output by the tick population model were converted to a binary variable, with all values 1 or greater (indicating sustainable) reclassified as 1, and all values less than 1 (indicating unsustainable) reclassified as 0. The binary outputs from the 14 years included in the study period were then summed.A summed R 0 value of 14 means that a region was sustainable throughout the entire study period, while a value between 0 and 13 inclusive, 14 indicates that it became sustainable between 2000 and 2013.

Lyme Disease Risk
The range of potential sustainable tick habitat was found to have expanded (Figure 3).The R0 values output by the tick population model were converted to a binary variable, with all values 1 or greater (indicating sustainable) reclassified as 1, and all values less than 1 (indicating unsustainable) reclassified as 0. The binary outputs from the 14 years included in the study period were then summed.A summed R0 value of 14 means that a region was sustainable throughout the entire study period, while a value between 0 and 13 inclusive, 14 indicates that it became sustainable between 2000 and 2013.The vast majority of southern and eastern Ontario has been consistently suitable for ticks, while northern Ontario has been consistently unsustainable for ticks during the same time period.However, there is a band separating the two regions where tick suitability has changed.This includes regions within the Niagara Escarpment and the Madawaska Highlands of Algonquin Provincial Park.

Model Validation
Model output was validated using tick submission data from Public Health Ontario.The majority of R0 values for submitted tick locations were above one, with a few occurring below one as could be expected due to some ticks identified in surveillance being ticks dispersed out of reproducing populations by migratory birds [31].Nevertheless, for the locations reporting R0 values below one, the R0 value was above one during the two years previous to the year when the tick was reported.This indicates that while a tick may have been collected in a year with a climate that was not sustainable for sustaining a tick population, the two years prior were sustainable for sustaining a tick population.
This visualization demonstrates that the majority of submitted ticks by the end of the study period were collected in areas having a climate consistently sustainable for sustaining tick populations (Figure 4).Comparison was conducted qualitatively, as it is apparent that the majority of cases occur within the always-sustainable regions of the map.However, it is important to note that The vast majority of southern and eastern Ontario has been consistently suitable for ticks, while northern Ontario has been consistently unsustainable for ticks during the same time period.However, there is a band separating the two regions where tick suitability has changed.This includes regions within the Niagara Escarpment and the Madawaska Highlands of Algonquin Provincial Park.

Model Validation
Model output was validated using tick submission data from Public Health Ontario.The majority of R 0 values for submitted tick locations were above one, with a few occurring below one as could be expected due to some ticks identified in surveillance being ticks dispersed out of reproducing populations by migratory birds [31].Nevertheless, for the locations reporting R 0 values below one, the R 0 value was above one during the two years previous to the year when the tick was reported.This indicates that while a tick may have been collected in a year with a climate that was not sustainable for sustaining a tick population, the two years prior were sustainable for sustaining a tick population.
This visualization demonstrates that the majority of submitted ticks by the end of the study period were collected in areas having a climate consistently sustainable for sustaining tick populations (Figure 4).Comparison was conducted qualitatively, as it is apparent that the majority of cases occur within the always-sustainable regions of the map.However, it is important to note that there were a small number of ticks which were submitted to Public Health Ontario in regions that became newly sustainable for carrying tick populations over the study period.These submitted ticks coincided with the expanded climatological range for sustainable tick populations predicted by the tick population model.In 2002, there were very few ticks submitted in Ontario.By 2012, the number of ticks submitted to Public Health Ontario had greatly increased.Furthermore, there was a geographic concentration of submitted ticks in eastern Ontario in the areas around Kingston and Ottawa.
there were a small number of ticks which were submitted to Public Health Ontario in regions that became newly sustainable for carrying tick populations over the study period.These submitted ticks coincided with the expanded climatological range for sustainable tick populations predicted by the tick population model.In 2002, there were very few ticks submitted in Ontario.By 2012, the number of ticks submitted to Public Health Ontario had greatly increased.Furthermore, there was a geographic concentration of submitted ticks in eastern Ontario in the areas around Kingston and Ottawa.

Discussion
While it is apparent that the number of Lyme disease cases in Ontario rose in the last decade [4,5], whether or not climate change has contributed to the increase requires further study.Much of the current research on Lyme disease risk mapping is focused on finding factors that correlate well with current tick presence [32], or projecting future tick habitat distribution due to climate change [2,33].This study examined whether there is evidence that the geographic distribution of Lyme disease risk in eastern Ontario has changed between 2000 and 2013 due to climate change, where the capacity to sustain an endemic tick population was used as a proxy for Lyme disease risk.The capacity to sustain a tick population was defined by cumulative annual degree days above 0 °C, and by changes in tick populations calculated using a tick population model.
First, it was determined that a warming trend has occurred in Ontario between 1979 and 2013, except for a small pocket in southern Ontario which has experienced cooling.This finding is consistent with previous research using data from Canadian meteorological stations [34].Then, it was shown that the vast majority of ticks submitted to Public Health Ontario were found in areas where the climate was sustainable for ticks from 2001 to 2014.There was a clear increase in submitted ticks in eastern Ontario in the last decade in areas already sustainable for a tick population.The increase in submitted ticks suggests the existing tick populations in these areas were stabilizing, rather than spreading to new geographic areas due to climate change.However, a small number of ticks submitted to Public Health Ontario during the study period were from newly sustainable regions, coinciding with the expanded climatological range for sustainable tick populations.This supports the idea that climate change has played some role in the changed spatial distribution of Lyme disease risk.When combined with Lyme disease risk mapping studies projecting climate change forward

Discussion
While it is apparent that the number of Lyme disease cases in Ontario rose in the last decade [4,5], whether or not climate change has contributed to the increase requires further study.Much of the current research on Lyme disease risk mapping is focused on finding factors that correlate well with current tick presence [32], or projecting future tick habitat distribution due to climate change [2,33].This study examined whether there is evidence that the geographic distribution of Lyme disease risk in eastern Ontario has changed between 2000 and 2013 due to climate change, where the capacity to sustain an endemic tick population was used as a proxy for Lyme disease risk.The capacity to sustain a tick population was defined by cumulative annual degree days above 0 • C, and by changes in tick populations calculated using a tick population model.
First, it was determined that a warming trend has occurred in Ontario between 1979 and 2013, except for a small pocket in southern Ontario which has experienced cooling.This finding is consistent with previous research using data from Canadian meteorological stations [34].Then, it was shown that the vast majority of ticks submitted to Public Health Ontario were found in areas where the climate was sustainable for ticks from 2001 to 2014.There was a clear increase in submitted ticks in eastern Ontario in the last decade in areas already sustainable for a tick population.The increase in submitted ticks suggests the existing tick populations in these areas were stabilizing, rather than spreading to new geographic areas due to climate change.However, a small number of ticks submitted to Public Health Ontario during the study period were from newly sustainable regions, coinciding with the expanded climatological range for sustainable tick populations.This supports the idea that climate change has played some role in the changed spatial distribution of Lyme disease risk.When combined with Lyme disease risk mapping studies projecting climate change forward [2,35,36], these results indicate that Lyme disease risk mapping will only become more important in the future to aid in early diagnosis of Lyme disease cases.
The following conclusions can be made from this study: 1.
Climate in Ontario between 1979 and 2013 has been warming except for a small pocket in southern Ontario, which has been gradually cooling.

2.
The changing climate in Ontario has increased in areas in Ontario that are sustainable for carrying a reproductive tick population, given stable host finding for ticks.
This study has also shown that remote sensing data can be used to increase the spatial resolution for Lyme disease risk mapping.The MODIS climate data used has a resolution of 0.05 × 0.05 • , making the findings more applicable at a local and regional scale.Previous studies on the impact of climate change on Lyme disease risk in eastern Ontario and northeastern United States have used climate models that have coarser resolutions.Brownstein, Holford and Fish used the Canadian Global Coupled Model with a resolution of 3.75 × 3.75 • (interpolated and resampled to 0.5 • ) while Ogden, Maarouf and Barker used the Coupled Global Climate Model with a resolution of approximately 3.7 • [2,36].While valuable for identifying areas of concern for potential Lyme disease risk monitoring, the coarse resolution can be of limited value to government organizations responsible for monitoring and mitigating Lyme disease risk at a more local level.In another study by Ogden, Bigras-Poulin and O'Callaghan, weather station data was used to model R 0 values [3], which can create accurate but geographically sparse results [37].Interpolation of these weather stations can introduce error and uncertainty in studying Lyme disease risk over large geographic areas.This study demonstrates that MODIS remote sensing data can be used to yield high resolution, continuous geographic coverage of Ontario to determine changes in climate and its effects on Lyme disease risk.
Higher resolution Lyme disease risk mapping can better identify areas of concern at local and regional scales than temperature data using coarser resolution, which may better assist with Lyme disease risk monitoring.Monitoring of Lyme disease at higher resolution will become of greater interest as Lyme disease continues to progress northwards in Ontario.This allows for more accurate risk mapping, which is extremely valuable on a local scale.
These results coincide with previous work [3,28] indicating that the theoretical niche of ticks is increasing rapidly, which is important for public health practitioners to be aware of.Also, these outputs are in line with previous findings [4,8] indicating that the spatial distribution of Lyme disease risk has spread in a northwestern fashion in Ontario over the past decade, progressing northward from the United States.This supports climate change models indicating that Lyme disease risk will continue to move northwards [3,38].The findings of this study have implications in the future study of Lyme disease as it gives insight into which regions may experience increased risk in the future.Since the first identified Lyme-positive tick in Ontario was found in 1993 [39], the regions in which ticks can be found has continued to expand.Climate models have historically been an importance predictor of future conditions.This study contributes to the history of population and climatic models that were used previously to study vector-borne disease [3,18,28,38].This type of predictive modelling has significant clinical implications as common Lyme disease tests lack specificity, meaning that false positives are prevalent [40].Therefore, whether individuals live in or have travelled to a Lyme-endemic region remains an important diagnostic criterion.
There are a few limitations in this study which should be noted.Firstly, the ability for a region to sustain an I. scapularis population is treated as a proxy for Lyme disease risk.In some cases, adventitious ticks may be present due to deposition by migratory hosts, or established tick populations may not be infected with the B. burgdorferi bacterium [11].However, due to the large spatial range and small physical size of I. scapularis, it is logistically difficult to directly measure their habitat range.Thus, proxies such as sustainable climate are valuable.
Also, migration was omitted in the tick population model because ticks have little horizontal movement.However, host migration may play a role in introducing ticks into new areas where an endemic tick population does not exist but climatic conditions are sustainable [41].Furthermore, the tick population model assumes a constant host density throughout the study region.Further research is needed on the impact of changes in host distribution and migration as a possible factor for increased Lyme disease cases in southeastern Ontario.
Ticks collected by health practitioners and submitted to Public Health Ontario were used in this study as a proxy for locations of Lyme disease risk.This methodology relies on adherence by health practitioners to submit all ticks that were collected.However, it is possible that health practitioners who collect numerous ticks underreported the number of ticks collected due to the high frequency, and it is also possible that practitioners who rarely collect ticks were unaware of the need to report.Possible lack of reporting could manifest in differences in spatial distribution of collected ticks than what was available and mapped.
Given that the majority of ticks submitted to Public Health Ontario are from previously sustainable regions, this reinforces the importance of future research into other potential causes of a Lyme disease risk increase.This suggests the increased incidence of Lyme disease may be predominantly due to factors such as tick populations stabilizing within their existing habitats, or an increase in tick-human interactions.Several of these factors have been previously studied, in terms of how they correlate with I. scapularis presence.For instance, in the Thousand Islands region, Ontario, host-species abundance was found to be strongly positively correlated with tick abundance, while species richness, a measure of biodiversity, is inversely associated [32].Also, in the eastern United States, altitude and monthly mean vapour pressure were found to be predictors of tick density [38].In New York State, low precipitation, low forest cover, and high urban density are all positively correlated to tick abundance [42].These factors may play a role in the stabilizing tick populations observed during this study.
In order to build upon the findings of this study, further research is needed into these other potential determinants of Lyme disease risk distribution.For instance, greater insight is needed into vector host activity, including density and migration patterns.One of the potential reasons for the observed increase in Lyme disease prevalence is an increase in the density of host species such as white-tailed deer and white-footed mice throughout southern Ontario.Additionally, host density may itself be impacted by climate change, further complicating this relationship.Also, ticks are often introduced into new regions by migratory hosts, so greater insight into the migratory patterns of the main host species would be beneficial.Additionally, human development patterns should be modelled on smaller scales in order to assess the role of changes in the frequency of human-tick interactions to changing Lyme disease risk.

Conclusions
This study used MODIS remote sensing data to evaluate the changing niche of I. scapularis in southern Ontario.It was found that temperatures have steadily increased throughout the province, making climatic conditions sustainable for tick populations over a greater area.The northernmost boundary at which ticks can survive expanded further north between 2000 and 2013.This indicates a greater potential risk area for Lyme disease, particularly in the Niagara Escarpment and Algonquin Provincial Park.The majority of ticks submitted to Public Health Ontario over this period were from regions that were sustainable throughout the entire study period, indicating the importance of other factors such as host abundance and human development in assessing risk.Greater insight into the factors driving the rapid increase in Lyme disease risk in Ontario are valuable to public health practitioners, as they enable better predictions to be made about which regions may see an increase in risk.

Figure 1 .
Figure 1.The study area boundaries delineate the region of Ontario included for analysis.The northern boundary is defined by row 840 in the MODIS11C3 data, which encompassed all of southern Ontario and part of northern Ontario that coincided with ticks collected and submitted to Public Health Ontario.

Figure 1 .
Figure 1.The study area boundaries delineate the region of Ontario included for analysis.The northern boundary is defined by row 840 in the MODIS11C3 data, which encompassed all of southern Ontario and part of northern Ontario that coincided with ticks collected and submitted to Public Health Ontario.

Figure 2 .
Figure 2. Change in mean monthly temperature in southern Ontario between 1979 and 2013.Data was obtained from the MOD11C3 land surface temperature product.

Figure 2 .
Figure 2. Change in mean monthly temperature in southern Ontario between 1979 and 2013.Data was obtained from the MOD11C3 land surface temperature product.

Figure 3 .
Figure 3. Summed basic reproductive number for tick populations in southern Ontario from 2001 to 2014.The values were reclassified as binary, with 1 for sustainable and 0 for unsustainable.A value of 14 indicates that a region was sustainable throughout the entire study period.A value between 1 and 13 inclusive indicates that a region became sustainable during the study period.

Figure 3 .
Figure 3. Summed basic reproductive number for tick populations in southern Ontario from 2001 to 2014.The values were reclassified as binary, with 1 for sustainable and 0 for unsustainable.A value of 14 indicates that a region was sustainable throughout the entire study period.A value between 1 and 13 inclusive indicates that a region became sustainable during the study period.

Figure 4 .
Figure 4.The number of ticks submitted to Public Health Ontario (normalized according to road length to account for population accessibility differences) and plotted over the tick suitability output by the model.Red represents always sustainable, green regions became sustainable during the study period, and black areas are consistently unsustainable.The number of ticks was multiplied by a factor of 100 for clarity.

Figure 4 .
Figure 4.The number of ticks submitted to Public Health Ontario (normalized according to road length to account for population accessibility differences) and plotted over the tick suitability output by the model.Red represents always sustainable, green regions became sustainable during the study period, and black areas are consistently unsustainable.The number of ticks was multiplied by a factor of 100 for clarity.