Quantifying the Northward Spread of Ticks (Ixodida) as Climate Warms in Northern Russia

: Climate change is affecting human health worldwide. In particular, changes to local and global climate parameters inﬂuence vector and water-borne diseases like malaria, dengue fever, and tick-borne encephalitis. The Republic of Sakha in northern Russia is no exception. Long-term trends of increasing annual temperatures and thawing permafrost have corresponded with the northward range expansion of tick-species in the Republic. Indigenous communities living in these remote areas may be severely affected by human and livestock diseases introduced by disease vectors like ticks. To better understand the risk of vector-borne diseases in Sakha, we aimed to describe the increase and spatial spread of tick-bite cases in the Republic. Between 2000 and 2018, the frequency of tick bite cases increased 40-fold. At the start of the period, only isolated cases were reported in southern districts, but by 2018, tick bites had been reported in 21 districts in the Republic. This trend coincides with a noticeable increase in the average annual temperature in the region since the 2000s by an average of 1 ◦ C. Maps illustrate the northward spread of tick-bite cases. A negative binomial regression model was used to correlate the increase in cases with a number of climate parameters. Tick bite case frequency per district was signiﬁcantly explained by average annual temperature, average temperature in the coldest month of the year, the observation year, as well as Selyaninov’s hydrothermal coefﬁcient. These ﬁndings contribute to the growing literature that describe the relationship between tick abundance and spread in Northern Latitudes and changes in temperatures and moisture. Future studies might use these and similar results to map and identify areas at risk of infestation by ticks, as climates continue to change in Sakha.


Introduction
The Arctic and neighboring northern regions are among the areas where climate change and its effects are predicted and have been felt most strongly [1,2]. These dramatic changes have sparked concern about the possible impacts of climate change on the economies, public health, and traditional way of life in the cold northern regions [3,4]. In response, a number of studies were conducted to identify and address the issues faced by northern communities [5]. From these, it became apparent that traditional communities in the cold north were among the most vulnerable to the effects of changing climate. There is a need for research to quantify and predict the effects that changing climate will have on dif-between temperatures and tick activity can explain the recently observed northward spread of ticks in various regions, as cold climates become more temperate.
Changing tick distributions, increased cases of tick bites, and tick-borne diseases observed in other northern regions should be a cause of concern for communities in Sakha since climate trends here are similar to the other affected regions. In Siberia, the majority of reported tick bites are from Taiga ticks (Ixodes persulcatus), but cases from other species like Dermacentor nutalli are reported [23]. About 25% of I. persulcutus collected in this region were infected by pathogens, including B. burgdorferi (the cause of Lyme disease), TBEV, and A. phagocytophylum, which causes anaplasmosis in livestock [23]. The presence of new disease vectors can expose communities and their native livestock to novel pathogens against which they have no resistance and for which counter-measures are not yet in place. New tick-borne diseases might be a threat to traditional reindeer and cattle husbandry-an important economic activity among communities in rural Yakutia [6,8].
In light of the mentioned points, this study aimed to quantify recent changes in the frequency and distribution of human tick-bite cases in Sakha. We correlate the annual tick-bite frequency per district (Administrative level 2) with meteorological variables of temperature and precipitation that have previously been identified to influence tick abundance. Quantifying which variables are contributing to the spread and increase of tick bite cases might contribute to understanding future tick-bite risks in different parts of the Republic, as the climate continues to change. Quantitative results can assist various stakeholders, including public health and agricultural entities, to adopt appropriate strategies to deal with emerging disease-related challenges. Finally, the large study area and long study period offers insight into the spread of Ixodida ticks in relation to climate changes at a regional scale [24].

Study Area
The Republic of Sakha, also known as Yakutia (Figure 1) is considered one of the coldest populated regions of the world, with very low temperatures in winter reaching −50 • C and lower. The towns Oymyakon and Verkhoyansk (Figure 1) are the coldest inhabited places on earth, with minimum recorded temperatures below −67 • C [25]. However, in recent years there has been a rise in the average annual temperature, which is currently registered everywhere in the world and in Russia. So, in regions of the Republic, temperatures of −60 • C and lower are becoming rarer, and in the central areas in wintertime, temperatures of −50 • C and lower are also becoming less common. Along with this, abnormally high temperatures in the summer period are often recorded, which are not usual for this region, reaching the level of 34 • C and higher, as well as the recording of changes in the average monthly temperature in the off-season periods [9].

Data
Data for the incidence of tick bites on humans were obtained from the Republic of Sakha's Federal Service for Surveillance on Consumer Rights Protection and Human Wellbeing (Rospotrebnadzor). We had access to total reported cases in the republic for the period 2000 to 2018, but district-level case data only between 2006 and 2018. The data is compiled annually by the health authority according to reports of tick bites from clinics in each district. The official state institution Rospotrebnadzor in the Republic of Sakha (Yakutia) compiles data from all the districts. The dataset does not distinguish tick species, but Taiga ticks (Ixodes persulcatus) are reported to be the most commonly encountered species [23].
Meteorological data from six weather stations (Figure 1) was provided by the Yakut Department of Hydrometeorology and Environmental Monitoring. For the final analysis, we considered only the annual average temperature measured at each station between 1961 and 2018. Four of the weather stations were located in central Sakha, in the cities of Yakutsk and Pokrovsk, and the villages Churapcha and Amga. The remaining two stations were in northern Sakha, in Verkhoyansk and Oymyakon.

Data
Data for the incidence of tick bites on humans were obtained from the Republic of Sakha's Federal Service for Surveillance on Consumer Rights Protection and Human Wellbeing (Rospotrebnadzor). We had access to total reported cases in the republic for the period 2000 to 2018, but district-level case data only between 2006 and 2018. The data is compiled annually by the health authority according to reports of tick bites from clinics in each district. The official state institution Rospotrebnadzor in the Republic of Sakha (Yakutia) compiles data from all the districts. The dataset does not distinguish tick species, but Taiga ticks (Ixodes persulcatus) are reported to be the most commonly encountered species [23].
Meteorological data from six weather stations (Figure 1) was provided by the Yakut Department of Hydrometeorology and Environmental Monitoring. For the final analysis, we considered only the annual average temperature measured at each station between 1961 and 2018. Four of the weather stations were located in central Sakha, in the cities of Yakutsk and Pokrovsk, and the villages Churapcha and Amga. The remaining two stations were in northern Sakha, in Verkhoyansk and Oymyakon.
We supplemented the weather station data with climate data from the ERA5 Climate Reanalysis Datasets [26], provided by the Copernicus Climate Data Store, and accessible through Google Earth Engine [27]. From the ERA5-Land monthly averaged-ECMWF climate reanalysis dataset, we extracted the monthly average ground temperatures at a resolution of 0.1 arc degrees. We used the Google Earth Engine platform to query and spatially aggregate the data into averages for each district (administrative level 2) in the Republic of Sakha. Additionally, we obtained daily air temperature and total daily precipitation between 2016 and 2018 at a resolution of 0.25 Arc Degrees from the ERA5 Daily Aggregated Climate Reanalysis. We spatially averaged the data to the district level.

Statistical Analysis
The data analysis was conducted initially in Microsoft Excel [28], but the final analysis was done with R 4.0.3 [29]. Figures were produced with the package ggplot2 [30]. Analysis and mapping of spatial data was done with the R packages sf and tmap, respectively [31,32]. Functions for generalized model fitting were obtained from the package MASS [33]. To describe long-term air temperature trends in Yakutia, we fit linear regression models to the annual average temperature data. We grouped together temperature observations from stations in the colder northern Yakutia (Verkhoyansk and Oymyakon), and we grouped together observations from the four weather stations in central Yakutia (Yakutsk, Pokrovsk, Churapcha, and Amga). This was to separately investigate the trend in colder northern and warmer central regions. We produced scatter plots and fitted regression lines using ggplot2.
To quantify the relationship between climate variables and the number of reported tick bite cases, we used a negative binomial regression model with a log-link function. We chose this model since our response variable counts data that showed over-dispersion. Our response variable was the number of tick bite cases reported per district per year. We considered the following predictor variables, all evaluated at the district level:

2.
Accumulated daily temperatures on days with an average temperature above 5 • C; this serves as a proxy for the warmth and length of the summer season, when ticks are most active.

3.
Average temperature in the hottest month of the year; this is also a proxy for the intensity of the warm season.

4.
Average temperature in the coldest month; the harshness of the winter might influence tick larvae survival and might dictate the presence of host species. 5.
The observation year (2006-2018); to account/control for possible non-climatic factors that change with time and influence the number of reported tick bite cases. These could include changing population behavior, or medical reporting quality. 6.
Selyaninov's Hydrothermal Coefficient; previous studies observed that Taiga tick abundance in Russia was well explained by this coefficient. It is defined as the sum of daily precipitation (mm) on days with average temperature above 5 • C, divided by the sum of daily average temperatures on the days with temperatures above 5 • C, then multiplied by 10.
In the model, we included all the variables as covariates, and systematically removed non-significant variables, until all remaining variables were significant, and the lowest AIC was obtained. However, collinearity between predictor variables could result in the significance of individual covariates to be misreported [34]. For example, correlation between the annual average temperature and temperature in the coldest month could result in one of the two being reported as 'non-significant' in a regression model that includes both as covariates. To account for this effect, we also fit separate models for each climate variable and year. We included a year as a covariate in each model to control for non-climatic factors that could be affecting tick bite reports independently of climate change. We then compared the models with the Akaike Information Criterion (AIC). We also considered, but rejected, models with other error distributions. Linear regression was inappropriate due to data being counts. Poisson regression was considered, but rejected due to over-dispersion in the response variable violating the assumptions of the Poisson regression model. Generalized linear mixed models were considered, to investigate the presence of random effects introduced by the observation year or district. However, the complex nature of these models resulted in non-convergence issues when multiple predictor variables and non-normal error distributions were considered. For this reason, we selected the negative binomial regression model as described above.

Climate Trends
During the considered period of observations from 1961 to 2018, we observed a significant positive upward trend in annual average temperatures at the weather stations in both central and northern Sakha ( Figure 2). The temperature increased on average by 0.058 • C per year at the central stations, and 0.049 • C per year at the northern stations.
Atmosphere 2021, 12, x FOR PEER REVIEW 6 of 16 climate variable and year. We included a year as a covariate in each model to control for non-climatic factors that could be affecting tick bite reports independently of climate change. We then compared the models with the Akaike Information Criterion (AIC). We also considered, but rejected, models with other error distributions. Linear regression was inappropriate due to data being counts. Poisson regression was considered, but rejected due to over-dispersion in the response variable violating the assumptions of the Poisson regression model. Generalized linear mixed models were considered, to investigate the presence of random effects introduced by the observation year or district. However, the complex nature of these models resulted in non-convergence issues when multiple predictor variables and non-normal error distributions were considered. For this reason, we selected the negative binomial regression model as described above.

Climate Trends
During the considered period of observations from 1961 to 2018, we observed a significant positive upward trend in annual average temperatures at the weather stations in both central and northern Sakha ( Figure 2). The temperature increased on average by 0.058 °C per year at the central stations, and 0.049 °C per year at the northern stations.     Figure 1. In Verkhoyanskiy, spring temperatures increased by between 3 • C and 4 • C, but winter temperatures did not increase. In Yakutsk gorsovet, average temperatures increased in all months except February.
Atmosphere 2021, 12, x FOR PEER REVIEW 7 of 16 two districts, Yakutsk gorsovet and Verkhoyanski gorsovet (Figure 4). The location of the central towns in these districts are shown in Figure 1. In Verkhoyanskiy, spring temperatures increased by between 3 °C and 4 °C, but winter temperatures did not increase. In Yakutsk gorsovet, average temperatures increased in all months except February.

Tick Bite Results
Aggregated at the republic level, the number of reported tick bite cases increased rapidly between 2000 and 2018 (Figure 1). A simple trend line reports an increase of 17.5 tick bite cases per year.

Tick Bite Results
Aggregated at the republic level, the number of reported tick bite cases increased rapidly between 2000 and 2018 ( Figure 1). A simple trend line reports an increase of 17.5 tick bite cases per year.

Spatial Distribution of Tick Bite Cases
In the year 2000, only isolated cases of tick bite reports on humans were recorded in Yakutia ( Figure 5), but by 2006, bites were registered in 6 districts of the Republic (Table 1; Figure 6). In 2012 and 2013, attacks were reported in 14 regions of the Republic, including not only the southern, but also the central regions. In 2015-2016, tick bites were recorded in 18 districts of the Republic. For the first time, a case was also reported in Verkhoyanskiy-north of the polar circle. By 2018, tick bite cases had been recorded in 21 different districts (Table 1). The cases showed a northward movement, as can be seen on the map in Figure 7.

Spatial Distribution of Tick Bite Cases
In the year 2000, only isolated cases of tick bite reports on humans were recorded in Yakutia ( Figure 5), but by 2006, bites were registered in 6 districts of the Republic (Table 1; Figure 6). In 2012 and 2013, attacks were reported in 14 regions of the Republic, including not only the southern, but also the central regions. In 2015-2016, tick bites were recorded in 18 districts of the Republic. For the first time, a case was also reported in Verkhoyanskiynorth of the polar circle. By 2018, tick bite cases had been recorded in 21 different districts ( Table 1). The cases showed a northward movement, as can be seen on the map in Figure 7.        Figure 7 shows maps of the best fitting model's residuals at 4 year intervals between 2006 and 2018. This was the model with year, average annual temperature and Selyaninov's HTC as predictor variables. The maps suggest no persistent spatial pattern in the residuals that would indicate some important but unaccounted for spatial driver of tick bite cases.

Discussion
Global warming is often cited as one of the important factors contributing to the range of expansion of ticks (order Ixodida) to the northern latitudes [4,24,35]. In the Republic of Sakha, temperature increases are particularly noticeable in both the central and northern regions (Figures 2 and 4). We observed an increase of about 0.5 °C per decade between the 1960′s and 2018. Other sources describe similar patterns in the Russian

Model Results
In the negative binomial regression models where each climate variable was considered separately, we found that all climate variables significantly explained the number of reported tick bite cases, except for the two variables that described the climate conditions during Summer, namely average temperature in the hottest month of the year and accumulated daily temperatures above 5 • C ( Figure 6). The other four predictor variables, namely, average annual temperature, temperature in the coldest month, Selyaninov's hydrothermal coefficient, and observation year were highly significant (p < 0.01). A model that included both average annual temperature and Selyaninov's HTC fit the data the best, and had the lowest AIC ( Figure 6). For this model, coefficient estimates and 95% confidence intervals are shown in Table 2. Shown estimates and intervals are transformed to show true response (i.e., to transform back from the log-response used in the negative binomial model).  Figure 7 shows maps of the best fitting model's residuals at 4 year intervals between 2006 and 2018. This was the model with year, average annual temperature and Selyaninov's HTC as predictor variables. The maps suggest no persistent spatial pattern in the residuals that would indicate some important but unaccounted for spatial driver of tick bite cases.

Discussion
Global warming is often cited as one of the important factors contributing to the range of expansion of ticks (order Ixodida) to the northern latitudes [4,24,35]. In the Republic of Sakha, temperature increases are particularly noticeable in both the central and northern regions (Figures 2 and 4). We observed an increase of about 0.5 • C per decade between the 1960 s and 2018. Other sources describe similar patterns in the Russian Arctic and also Sakha [1, 4,9]. The frequency of tick-bite cases in Sakha has increased with the increase of temperature ( Figure 5, Table 1). According to our regression model results, average annual temperature in a district was significantly and positively correlated to the frequency of tick bites in the district. However, it is known that tick activity and abundance is in particular influenced by temperatures and moisture during the warm season, since tick host seeking activity, survival, and developmental rates are closely related to the conditions in the summer when they are most active [20]. In correspondence with this, Opisova et al. [21] found that taiga tick activity was, among multiple meteorological predictor variables, best explained by temperature and precipitation in the periods of the year when temperatures reached above 5 • C. In contrast to this, however, our model reported that the sum of daily temperatures above 5 • C was not significant. The average temperature in the warmest month of the year, which was anticipated to relate to tick host-seeking activity in summer, was also not significant in our model. The only significant 'season specific' predictor variable was Selyaninov's hydrothermal coefficient. This corresponds to the earlier findings by Opisova et al. [21]. Selyaninov's HTC gives a measure of the precipitation in the warm part of the year. This suggests that increased summer precipitation should increase tick abundance. According to Gorkohov and Fedorov [9], patterns of changing precipitation are spatially variable in Sakha, with precipitation decreasing in the tundra landscapes, but increasing in 70% of the republic. The harshness of the winter is considered to be a constraining factor for the distribution of ticks, since nymph and larvae survival might reduce when winter temperatures drop below certain thresholds. Moreover, harsh winters might constrain the distribution of possible host species like deer [17]. Our analysis showed a significant relationship between average temperature in the coldest month of the year and tick bite frequency. More temperate winters might, over a number of years, enable host species to settle and tick populations to become established in new areas. Model predictions suggest a northward shift in the distribution of several Russian forest species in accordance with changing climate [36]. However, at the scale of one year, tick bite cases in an area are mostly related to tick host-seeking activity in the warm season [17,36], which is not largely influenced by the winter temperature variations. Host-seeking success at this time might be influenced by humidity. How frequently ticks need to return to water sources to rehydrate is dependent on local humidity, thus low humidity might reduce the time that ticks have to search for hosts [17]. In this study, we did not include a measure of humidity as a predictor variable of tick-bite frequency. Although Selyaninov's HTC gives a measure of moisture availability in the habitat, more explicit humidity variables like saturation deficit or relative humidity should be incorporated in the future.
The number of reported tick bite cases are potentially affected by multiple other non-climatic factors and need not directly relate to the abundance of ticks within an area. In recent years, there have been much public concern and focus in literature about the increased prevalence of tick-borne infections like Lyme disease and Encephalitis, and their association with climate changes. However, Sumilo et al. [19] observed that increased tick-borne diseases in the Baltics could not only be attributed to changing climate. There can be multiple complex factors and interactions that affect exposure to vectors and their diseases, including changes in socio-economic and public health conditions and improved reporting [14]. We attempted to account for the effect of non-climatic factors by also including observation as an explanatory variable in our model. However, when we controlled for observation year, variation in the climate variables remained statistically significant predictors of tick bite frequency in districts. One non-climatic factor that will have a certain influence on cases, but that we did not consider, is exposure due to population size. In our analysis, the negative binomial model did not converge when population size of the district was included as an exposure/offset term. This may be due to the added complexity, in an already complex model formulation. Therefore, we only considered the number of observed cases, irrespective of local population size. In Table 2, we do, however, report cases per 100,000 population in each district. Due to the sparse population of most of Yakutia, inclusion of population size inflated the magnitude of tick-bite increase, except for the most densely populated districts like Yakutsk gorsovet-home to the capital city Yakutsk.
The range shift of ticks to the north is noted in many countries, including other regions of the Russian Federation. In Sweden and Norway for the period of 1994-2008 the tick boundary has moved more than 200 km to the north along the Baltic coast [18]. In the mountains in the north of the Czech Republic, where the temperature has increased by 1.4 • C in four decades, ticks appeared at an altitude of 1300 m above sea level [37]. In Russia, a marked expansion of I. persulcatus ticks to the north occurred in the Arkhangelsk region and the Republic of Komi [38,39]. In recent years, in particular, the spread of tick populations to the northern territories of Yakutia has been observed [5,35]. In 2008, an outbreak of pyroplasmosis (babesiosis) was recorded for the first time among domestic reindeer in Central Yakutia. The outbreak was further complicated by the parasitic diseases. At the same time, a high mortality rate among the adult reindeer population was noted: Up to 30% at the "Tomtor" site, and up to 70% at the educational-production facilities of the Yakut State Agricultural Academy. It was established that the causative agent of the disease was Babesia bovis, and the parasite was transmitted by Ixodes ricinus ticks.
These reported cases of novel tick-borne disease outbreaks, coupled with the described pattern of northward spread of tick bite cases point again to the risk of infectious disease spread among livestock and communities in the near Arctic regions, as mentioned and predicted in recent years [4][5][6]. Therefore, it is necessary to expand ongoing research on integrated monitoring of ongoing climate change and forecasting its possible impact on traditional industries in the North. The climate variables that influence tick abundance can be used to map areas in the republic where risk of tick bites will become higher in the near future due to climate change. Similar tick-bite risk maps have been made in other areas, by considering remote sensing based indicators of possible tick abundance [40] and citizen science based reports of exposure to ticks [41]. Risk maps, informed by these climate variables, can help local authorities to identify areas where farmers should be supported and made aware of the current and future risks of tick-borne diseases in their livestock.

Conclusions
The warming climate in the Republic of Sakha corresponds with the increase of disease vectors in the region. We observed a clear northward spread and a general increase in reported tick bite cases in Sakha. Between 2000 and 2018, cases increased 40-fold. Additionally, cases were initially noted only in southern districts but by 2018 twentyone districts had been affected, including Verkhoyanskiy-which is above the Arctic Circle. Our negative binomial regression model predicts that winter temperature, annual temperature, and Selyaninov's hydrothermal coefficient are significantly and positively related to increases in tick-bite frequency in the Republic. We suggest that public health entities as well as agricultural organizations in the Republic of Sakha should continue to monitor the spread of tick-borne diseases, and consider possible interventions to reduce its impact on local health and livestock farming.  Data Availability Statement: The datasets generated in the current study are available from the corresponding author on request.