Polymorphism of TLR2 in bank vole populations in North Eastern Poland is not associated with Borrelia afzelii infection prevalence

Polymorphism in innate immune genes in host populations can structure spatial variation in the prevalence of infectious diseases. In Europe, Borrelia afzelii is an important tick-borne pathogen of small mammals including the bank voles (Myodes glareolus). The Toll-like receptor 2 (TLR2) is an innate immune receptor that is important for detecting Borrelia burgdorferi sensu lato pathogens. The TLR2 gene is polymorphic in bank vole populations and is classified into four distinct clusters: C1, C2, C3, and C4. The C2 and C4 clusters versus the C1 and C3 clusters are associated with lower versus higher infection prevalence, respectively. We detected three TLR2 clusters in 487 bank voles from 30 populations in NE Poland: 84.2% of the obtained sequences belonged to the C1 variant, 7.2% to C2, and 8.6% to C3. However, no clear spatial structure of TLR2 clusters among the populations was detected. B. afzelii infection prevalence across all studied individuals was 12.1% and varied from 0 to 37.5% among populations. There were no significant differences in B. afzelii prevalence among voles carrying alleles of different TLR2 clusters, or between individuals belonging to two mtDNA lineages. Most infected individuals were adults, and males were infected more often than females. There was no significant relationship between the prevalence of TLR2 clusters in the vole populations and climatic and environmental factors within the study area. We therefore could not confirm an adaptive role of the TLR2 C2 alleles in reducing B. afzelii infection prevalence in bank voles.


Introduction
Lyme disease (LD) is the most common tick-borne disease in the human population across Europe, Asia and North America (Kilpatrick et al. 2017;Estrada-Peña et al. 2018). The causative agents of LD are some species of the bacterial Borrelia burgdorferi sensu lato complex (Hanincová et al. 2006;Franke et al. 2013;Hofmeester et al. 2016;Estrada-Peña et al. 2018). Members of this complex can infect various vertebrate species, including humans (Coipan et al. 2018). However, individual Borrelia species are specialized on different vertebrate hosts (Kurtenbach et al. 2002;Piesman and Gern 2004). For example, Borrelia garini and B. valaisiana are specific to birds (Hanincová et al. 2003b) while B. afzelii is specific to small mammals, including rodents, which are the main reservoir and source of B. afzelii (Hanincová et al. 2003a;Hofmeester et al. 2016).
B. afzelii is the most common cause of LD in humans in central Europe (Estrada-Peña et al. 2018). This bacterium is transmitted via infected ticks, with the most important tick vector for B. afzelii in Europe being the castor bean tick (Ixodes ricinus) (van Duijvendijk et al. 2015;Hofmeester et al. 2016;Strnad et al. 2017). One of the most important rodent hosts for B. afzelii in Europe is the bank vole (Myodes glareolus) (Tälleklint et al. 1993;Kurtenbach et al. 1995;Humair et al. 1999). The spread of Borrelia in tick populations is determined by the specific infectivity of voles, seasonality, the degree of acquired resistance to ticks and the age composition of the vole population Tälleklint et al. 1993).

Communicated by: Jeremy Herman
Electronic supplementary material The online version of this article (https://doi.org/10.1007/s13364-020-00518-7) contains supplementary material, which is available to authorized users.
Toll-like receptors (TLRs) are part of the vertebrates' antimicrobial innate immunity, which recognize a variety of microbial products such as peptidoglycan, lipoproteins, lipopolysaccharides, lipoteichoic acids, bacterial flagellin and CpG DNA (Hirschfeld et al. 1999;Medzhitov 2001;Wooten et al. 2002). Their presence is the first signal of infection for the organism. Such recognition of microbial products by TLRs leads to functional maturation of dendritic cells and to initiation of antigen-specific adaptive immune responses (Medzhitov 2001). The main function of TLRs is to recognize pathogens, thus TLRs evolved under parasite-driven selection (Kloch et al. 2018). This natural selection has shaped the genetic diversity of host populations.
Toll-like receptor 2 (TLR2) plays a key role in the response to Borrelia infection (Hirschfeld et al. 1999;Wooten et al. 2002;Garantziotis et al. 2008). TLR2 is specific to peptidoglycan and lipoproteins and cooperates with TLR1 and TLR6 in the recognition of distinct sets of ligands (Hirschfeld et al. 1999, Medzhitov 2001. The interaction between lipoproteins and TLR2 activates cells involved in innate defence, prompting infiltration of professional phagocytes, phagocytosis or killing of spirochetes (Wooten et al. 2002).
Recent studies of the relationships between B. afzelii prevalence and bank vole populations have investigated parasitemediated selection, and the role of TLR2 in shaping infection prevalence within and across European populations of voles (e.g. Tschirren et al. 2013;Tschirren 2015;Morger et al. 2015;Gomez-Chamorro et al. 2019a). Population genetic analysis of the TLR2 genes of the bank vole may give insights into parasite-mediated selection, and could elucidate the role of parasites in creating and maintaining diversity of immune genes within and among host populations (e.g. Haldane 1949;Tschirren et al. 2011). It has been suggested that 68% of the TLR2 sequence evolved under strong purifying selection, 32% evolved neutrally, and 0.5-0.6% of all sites (3-5 amino acids) were shaped by positive selection (Tschirren et al. 2011).
In bank voles, four TLR2 clusters: C1, C2, C3, C4 (groups of similar haplotypes, which differ from each other by single nucleotide polymorphisms) are associated with variation in infection prevalence of B. afzelii (Tschirren et al. 2013;Morger et al. 2015). Field studies found that individuals with the C2C2, C1C2 and C1C1 genotypes had low (approximately 15%), intermediate (30%) and high (50%) prevalence of infection with B. afzelii (Tschirren et al. 2013;Tschirren 2015). This led the authors (Tschirren et al. 2013) to suggest that the C2 allele provides resistance to B. afzelii infection in bank voles whereas the C1 allele makes bank voles more susceptible to B. afzelii infection.
Haplotypes belonging to the TLR2 C4 cluster might increase resistance to B. afzelii infection in bank voles in a similar manner to TLR2 C2 . There is no apparent association between TLR2 C3 haplotypes and Borrelia infection (Cornetti et al. 2018). However, a recent study by Gomez-Chamorro et al. (2019a), who experimentally infected bank voles carrying different TLR2 variants (C1, C2 and C3) with Borrelia afzelii, showed no differences in the prevalence of B. afzelii infection among bank voles carrying alleles from different TLR2 clusters. To date, several hypotheses explaining the observed TLR2 genetic structure among populations have been proposed, concerning natural selection and the environment. For example, co-occurrence of three TLR2 clusters within the same population might reflect selection promoting polymorphism (Morger et al. 2015). Spatial distribution of different TLR2 clusters may be related indirectly to abiotic factors, which influence the presence of ticks and Borrelia infection prevalence. Cornetti et al. (2018) showed that in Switzerland, Borrelia infection risk was highest at low and medium elevations, and decreased with increasing elevation, which corresponded proportionally to host infection prevalence. They suggested that such small-scale spatial variation in pathogen-mediated selection shaped host gene evolution, leading to small-scale spatial variation in the immunogenetic composition of host populations. Additionally, differences in the strength of selection might lead to spatial heterogeneity in the resistance, tolerance, and life history of the host species (Cornetti et al. 2018). The Borrelia infection risk depends on the density of infected nymphs (DIN), which is shaped by both abiotic and biotic factors, and varied at a local spatial scale. There is a complicated network of factors influencing transmission of the pathogen between hosts and vectors (Ostfeld et al. 2006;Diuk-Wasser et al. 2012;Perez et al. 2016;Cayol et al. 2017). The number of nymphs depended on the density of rodents (Ostfeld et al. 2006). Host-parasite interactions also depend on the type of landscape and its fragmentation (Diuk-Wasser et al. 2012;Perez et al. 2016). The distribution and abundance of tick populations can also be shaped by both large-scale and local climatic conditions (Ostfeld et al. 2006;Kilpatrick et al. 2017). Winter survival of ticks is mostly determined by extreme temperatures and occurrence of snow cover (Karbowiak 2014;Medlock et al. 2013;Medlock and Leach 2015). Insufficient snow cover during frosts can lead to ground freezing, a consequence of which is reduction of the tick population (Medlock and Leach 2015).
European bank vole populations are classified into mitochondrial DNA (mtDNA) lineages, which are defined as temporal populations connected by a continuous line of maternal descent based on mitochondrial markers. Previous studies in Europe have shown that the TLR2 clusters are conserved among bank vole mtDNA lineages suggesting that the formation of the TLR2 clusters predated the split of the bank vole lineages and that these clusters were subsequently maintained within lineages. For this reason, we were interested in studying the association between the mtDNA lineages and the TLR2 clusters in our bank vole populations. In NE Poland, mainly two mtDNA lineages of the bank vole occur: the Carpathian and the Eastern (Wójcik et al. 2010;Stojak and Tarnowska 2019). These two mtDNA lineages form a contact zone that varies in width from 45 km in the south to 243 km in the north; the shape and width of this contact zone has been related to both abiotic and biotic factors (Tarnowska et al. 2016;Stojak and Tarnowska 2019). Genetic differentiation of those lineages has been shaped during isolation in different glacial refugia and postglacial recolonisation of Europe (Deffontaine et al. 2005;Kotlík et al. 2006). The spatial distribution of these two bank vole mtDNA lineages in NE Poland (Tarnowska et al. 2016) corresponds to areas of varying levels of Borrelia infection risk in Ixodes ricinus ticks (from 1 to 37% in different areas) as indicated by Sytykiewicz et al. (2012).
In our study we tested the following hypotheses: (1) frequencies of the TLR2 clusters differ between bank voles belonging to the two main mtDNA lineages in NE Poland, where the level of Borrelia infection in ticks varied spatially (Sytykiewicz et al. 2012), (2) bank voles carrying different clusters and genotypes of TLR2 differ in B. afzelii infection prevalence. The main aims of this study were to (1) estimate B. afzelii infection prevalence in bank voles, (2) describe the genetic diversity of TLR2 and estimate the frequencies of the three TLR2 clusters, (3) compare spatial TLR2 population structure with mtDNA lineage distribution in bank vole populations in North Eastern Poland and, finally, (4) identify which biological, climatic and environmental variables might have shaped the prevalence of B. afzelii infection in bank voles in NE Poland.

Study area and sampling
This study was conducted in a forested area of North Eastern Poland (52°21′-54°20′ N, 18°59′-23°53′ E), covering approximately 45,000 km 2 (Fig. 1). The mean forest cover varied from 22% in the western part to 34% in the eastern part of the study area (Zielony and Kliczkowska 2012). North Eastern Poland is under the influence of two types of climate: Atlantic in the West, and continental in the East (Starkel 1991). The mean temperature in January ranges from − 3.8°C in the West to − 5.8°C in the East, based on mean values for 1950-2000 (Hijmans et al. 2005, the World Climate database). The mean depth of snow cover varies from 7 to 14 cm, and the e average number of days with snow cover during the calendar year varies from 44 to 80 in the western and the eastern parts of the study area, respectively (source of data: vectorized maps of mean snow fall from years 1970-2000 for Poland; Institute of Meteorology and Water Management -National Research Institute). According to Sytykiewicz et al. (2012), infection prevalence of the tick Ixodes ricinus with Borrelia spirochetes varies among different tick populations in this region. However, those previously studied populations did not overlap with our study sites.
We used tissue samples from 461 bank voles from 29 populations, collected in 2006 and 2011-2013, during previous projects conducted by the Mammal Research Institute (MRI) as described in Niedziałkowska et al. (2010) and Tarnowska et al. (2016). The bank voles were live-trapped and~4-mm samples of the end of the tails of each of the captured individuals were taken for genetic analyses. In 2014, we sampled 26 bank voles from an additional population (no. 29) according to the protocol described in Tarnowska et al. (2016); small ear fragments of each of the captured specimens were collected for the purposes of the study. In total, we analysed 487 tissue samples from bank voles from 30 localities. All analyses were conducted using this same set of 487 individuals. The number of specimens studied per population ranged from 10 to 26, which was sufficient to reliably estimate frequencies of TLR2 genetic clusters and mtDNA lineages in each population (Morger et al. 2015;Tschirren et al. 2015;Tarnowska et al. 2016). A list of all sampling sites and numbers of collected specimens is presented in Table 1. All samples were collected in concordance with permissions 07/2004, 15/2006, 16/2011, 43/2011, 12/2013 and 26/2014 obtained from the Local Ethical Commission in Białystok, Poland. Trapping sites were approved by the administrators of the study area.

Detection of Borrelia afzelii infection
Genomic DNA was extracted using Qiagen DNeasy Blood & Tissue Kit (Valencia, CA, USA). DNA extraction from the B. afzelii control sample was carried out using the Syngen DNA Mini Kit (Syngen Biotech), according to the manufacturer's protocol. We detected infection in bank vole samples by targeting an 89-bp fragment of the flaB gene with primers specific to B. afzelii. We used a modified version of the qPCR protocol and primers described by Råberg (2012). PCR was carried out in 10 μL volumes: 5 μL HotStartTaq MasterMix (Qiagen), 0.25 pmol of each primer, 1 μL of DNA, and water to a 10 μL volume. The PCR profile consisted of initial incubation at 95°C for 15 min and then 45 cycles of 95°C for 15 s, 59°C for 30 s and 72°C for 30 s. An additional sample of the reagent mix used to amplify the flaB gene constituted the first negative control. The second negative control contained RNAse-free water (Qiagen). For the positive controls, we used DNA isolated from inactivated B. afzelii spirochetes, kindly provided by the National Institute of Public Health (https://www.pzh.gov.pl/). The results of each PCR were checked on 1% agarose gel. All samples were amplified twice, and only individuals with two flaB-positive results were considered as infected.
To confirm that our PCR was detecting B. afzelii, we amplified 25% of the flaB-positive samples and sequenced them using Fla5F primer and the Big Dye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems), according to the manufacturer's protocol. All obtained sequences were compared with sequences available in GenBank using BLAST, and thereby confirmed as belonging to the flaB gene of B. afzelii.
To determine the concentration of bacterial DNA, which allowed us to detect infection in voles, we prepared serial dilutions of extracted B. afzelii DNA (max concentration 2.2 ng/μl, min concentration 2 × 10 −4 ng/μl) and amplified the fragment of the flaB gene as described above. It was possible to detect B. afzelii infection in bank voles when the concentration of bacterial DNA in a sample was 2 × 10 −4 ng/ μl and higher.

TLR2 analyses
To determine the TLR2 genotype, we used PCR to amplify an 1173-bp fragment of the TLR2 gene. We used two sets of primers, one that were previously described by Tschirren et al. (2011) and one pair of primers designed for this study (TLR2F_int: 5′ TATTTGAG A A A C T C A G C C T G T G 3 ′ a n d T L R 2 R _ i n t : 5 ′ TTGTTACTAACATCCAGCACC 3′). PCRs were set in 10 μL volumes. The reaction mix contained 5 μL HotStartTaq MasterMix (Qiagen), 0.25 pmol of each primer and 60 ng of DNA. The PCR profile consisted of initial incubation at 95°C for 15 min and then 37 cycles of 94°C for 30 s, 57°C for 30 s, and 72°C for 150 s, with final elongation step at 72°C for 10 min. PCR products were purified using CleanUp kit (A&A Biotechnology).
Next, we used MglTLR2F and TLR2F_int primers and the Big Dye Terminator (BDT) v3.1 Cycle Sequencing Kit (Applied Biosystems) to obtain TLR2 DNA sequences. The reaction mix contained 1.5 μL BDT, 1.5 μL 10X BDT buffer, 25 pmol primer, 10 ng of amplified and purified DNA. Sequencing was performed on ABI PRISM 3130xl sequencer (Applied Biosystems). Obtained sequences were read and checked manually in ChromasLite v2.4 (http://technelysium.com.au/) and aligned in BioEdit 7.0.5 (Hall 1999). For reference sequences, we used three haplotype sequences described and deposited in the GenBank database by Tschirren et al. (2012) (accession numbers: JN674541, JN674538 and JN674535). Haplotypes found in only one specimen were checked against the GenBank database and repeat sequenced if necessary, to avoid genotyping mistakes. MtDNA sequences of 461 individuals (from 29 populations) were taken from our previous study (Tarnowska et al. 2016). We used 424-bp fragment of the cytochrome b target to determine the mitochondrial lineage of each bank vole. The cytochrome b analysis of bank vole samples from one new population (population number 29; n = 26 bank voles) targeted the same 424bp fragment of the cytochrome b gene and was conducted according to the protocol by Tarnowska et al. (2016).

Statistical analyses
Variability of TLR2 TLR2 haplotype reconstruction was performed in DnaSP v 6.11.01 (Rozas et al. 2017) using the PHASE algorithm with default settings of a thinning interval of one, 100 burn-in iterations and 100 main iterations. The obtained sequences were submitted to the GenBank database (see Table S1 for accession numbers). Summary statistics of genetic diversity of the TLR2 gene (number of polymorphic sites-S, number of haplotypes-h, haplotype diversity-Hd, nucleotide diversity-π, average number of nucleotide differences-k, and number of amino acid variants-A), were calculated for each population and for the whole sample set using DnaSP v 6.11.01.
The phylogenetic relationships of all TLR2 sequences were reconstructed using Median Joining network analysis in PopArt v 1.7 (French et al. 2014), and compared with mtDNA lineage frequencies. We conducted a phylogenetic analysis of TLR2 haplotypes using the data obtained in this study. A maximum likelihood phylogenetic tree was constructed in Mega 7.0.26 (Kumar et al. 2016) using 10,000 bootstrap replications, with model K2 + G + I previously indicated as the best-fitting model by Mega 7.0.26 (having the lowest value of Bayesian Information Criterion), using 95 TLR2 sequences from the GenBank database (see Table S2 for accession numbers), and with Myodes rutilus sequence as an outgroup (GenBank accession number: HM215590). These sequences were aligned as previously described and cropped to the same length (1149 bp). Mega 7.0.26 was used to prepare a protein relationship reconstruction. We applied model JTT+G+I as the best-fitting model.
We compared the frequencies of different TLR2 clusters (C1, C2 and C3) to the frequencies of the two mtDNA lineages (the Carpathian and Eastern; Tarnowska et al. 2016) across the 30 studied populations (frequencies of TLR2 clusters and mtDNA lineages are shown in Table S3), using the rcor.test implemented in the "lme4" package R (Bates et al. 2015; http://www.R-project.org.).

Environmental data
Climate data, including mean temperature values, were obtained from the WorldClim database (Hijmans et al. 2005); available online (http://www.worldclim.org). Values of the climate variables for each the studied populations were calculated as the means across all sampling plots occupied by each population. The percentage cover of forest, water, grassland, and urban area were calculated for each of the populations, including a 2-km radius buffer around each area, using ArcGIS 10.2 (Redlands, CA: ESRI 2012) and the CORINE Land Cover 2006 dataset (European Environment Agency, Copenhagen 2007). The number of days with snow Factors affecting B. afzelii prevalence As an initial step in variable selection for generalized linear mixed model (GLMM) analysis, we tested for correlations between a range of genetic, climatic and environmental variables (Table S3), using the rcor.test function of the "lme4" package in R. To avoid correlated variables (Table S5) and to reduce their numbers we ran PCA analysis for climatic and environmental variables to find the subset of those variables, which best explained the observed variation among populations. The optimum PCA result included factors such as forest cover, winter temperature, thickness of snow cover, and snow retention time, which are factors that are known to infect the DIN and therefore the B. afzelii infection risk (Karbowiak 2014;Medlock and Leach 2015). The first two axes of the PCA explained 80.6% of the variation (Fig. S1). The PC1 reflected winter conditions: mean January temperature, snow depth and mean snow retention time, and the PC2 reflected forest cover. We identified 7 factors, which potentially could have affected the Borrelia afzelii prevalence among bank voles: (1) TLR2 genotype, (2) mtDNA lineage, (3) age, (4) sex, (5) winter conditions (PC1), (6) forest cover (PC2) and (7) trapping season (Table S4). We did not include the year of sampling, because 24 of the 30 populations were sampled in the same year (2012). We tested these selected factors with results from 393 bank voles, excluding 94 juveniles due to lack of sex identification (none of the excluded individuals were infected with B. afzelii).
Infection status (0 = uninfected, 1 = infected) are binomial data and therefore we used GLMMs with binomial errors to test the impacts of factors which may have affected B. afzelii prevalence among bank vole individuals. GLMMs include both fixed and random effects. We modelled the bank vole population as a random effect (30 populations -number of populations limiting total number of fixed effects in the model to seven variables), and compared GLMMs. We used data concerning the mtDNA lineage (Carpathian or Eastern), sex, and age class (distinguished based on body weight-adults ≥ 18 g and juveniles < 18 g; Pucek 1984). Winter conditions and the forest cover were included in the models as the values of the two first axes of PCA. Information on the timing of trapping was included, using four sampling time periods (I-the first half of July, II-the second half of July, III-the second half of August, IV-the first half of September). We did not trap bank voles in June and the first half of August. Information about TLR2 clusters was coded in two ways to compare these approaches. In the first way, we used information about TLR2 cluster genotypes, and in the second way, we used the number of C2 alleles (which are predicted to reduce susceptibility of B. afzelii infection).
We applied an information theoretic approach to model selection, using a second-order correction for small sample size, Akaike's information criterion-AIC c (Burnham and Anderson 2002), implemented in the 'MuMIn' R package (Bartoń 2014). All statistical analyses were completed in R (R Development Core Team 2016).

Borrelia afzelii infection prevalence in bank voles
Among 487 analysed samples of bank voles, we detected 59 Borrelia afzelii positive cases, meaning the overall infection prevalence of the study area was 12.1% (Table 2), which varied at the local level from 0 to 37.5% between populations ( Fig. 2 lower panel and Table S6). The results of a chi-squared test did not show any differences in infection prevalence between 30 studied bank vole populations (χ 2 = 35.7, df = 29, p = 0.18). Fifteen of the 59 bank vole samples that tested positive were sequenced for an 89-bp fragment of the flaB gene of B. afzelii. All 15 sequences matched to the reference sequence of B. afzelii (GenBank accession number: MG944963, positions: 400-488 bp).

MtDNA lineage distribution in NE Poland
The mtDNA lineage frequencies of studied populations were obtained in a previous study by our research group (Tarnowska et al. 2016). All specimens from the new population investigated in the current study belonged to the Eastern mtDNA lineage (haplotypes H2, H77, H48 and H78), described previously by Tarnowska et al. (2016). In total, 158 (32.4%) of the 487 individuals belonged to the Carpathian lineage and 329 (67.6%) to the Eastern lineage. We identified three pure Carpathian populations (in the westernmost part of the study area), eight pure Eastern populations (located in the eastern part of the study area) and 19 mixed populations of both lineages (Fig. 1). These two lineages formed a contact zone, which varied in width from 45 km (southern part of the study area) to 243 km (in the northern part) (Tarnowska et al. 2016).

Genetic variability of the TLR2 in bank voles
In total, we obtained 109 haplotypes of TLR2, with 75 polymorphic sites, including 87 haplotypes detected for the first time (see Table S1). Among all identified haplotypes, 76 (70%) belonged to cluster TLR2 C1 , 17 (15%) to TLR2 C2 and 16 (15%) to TLR2 C3 . The most frequent haplotype H4.24, belonging to cluster TLR2 C1 , was found in 26% of all individuals (Fig. S2). Ten of the most frequent haplotypes from all three TLR2 clusters were shared by the Carpathian  and Eastern mtDNA lineages (Fig. S2). After translation, we obtained 58 amino acid variants (Fig. S4), 41 for cluster TLR2 C1 , 9 for TLR2 C2 , and 8 for TLR2 C3 . The number of amino acid variants ranged from 5 to 12 among populations, while the numbers of private amino acid variants per population varied from 0 to 3. All considered parameters of genetic variability of the TLR2 gene (number of haplotypes, haplotype and nucleotide diversities and average number of nucleotide differences) were highly diverse among populations (Table 1). Among the analysed samples (487 individuals and 974 sequences) we identified the TLR2 variants belonging to three clusters (Figs. 2, S2 and S3): TLR2 C1 was found in 84.2% of sequences, TLR2 C2 in 7.2%, and TLR2 C3 in 8.6%. Most of the individuals (71.9%) were homozygotes of TLR2 C1 . Roughly equal numbers of specimens were heterozygotes of TLR2 C1 /TLR2 C2 (11.7%) and TLR2 C1 /TLR2 C3 (12.9%), and only 1% of individuals were heterozygotes of TLR2 C2 / TLR2 C3 . We also identified homozygotes of TLR2 C2 (0.8%) and TLR2 C3 (1.6%). The frequencies of the various TLR2 clusters did not significantly correlate with frequencies of mtDNA lineages in any of the 30 populations studied (r s from − 0.35 to 0.35; p from 0.06 to 0.98; Table S5).

Factors influencing B. afzelii infection prevalence in bank voles
Age and sex were indicated as the most important factors determining B. afzelii prevalence in bank voles according to the GLMM analysis (Table 3). Adult males were more frequently infected than adult females and juveniles (age: p < 0.001, sex: p < 0.01). AIC-based model selection indicated the negligible importance of environmental factors, TLR2 genotypes or number of TLR2 C2 alleles and distinct mtDNA lineages of bank vole in explaining B. afzelii prevalence in voles. The next-best 10 models were within 2 AIC units from the top model (Table 3). All of these models contained age and sex variables indicating that these were the critical explanatory variables. The rest of the explanatory variables in these models were not significant predictors of B. afzelii prevalence in bank voles.

Relationships between biological, climatic and environmental variables and B. afzelii infection prevalence
One of the goals of our study was to determine which biological, climatic and environmental variables determined the infection prevalence of B.afzelii in bank voles. We found that B. afzelii infection prevalence was lower in young than in old animals and higher in males than in females. These results are in agreement with previous studies (Tälleklint and Jaenson 1995;Hofmeister et al. 1999;Bunikis et al. 2004). The effect of age could be the result of a longer exposure to ticks and higher activity in older animals than in younger ones (Tälleklint and Jaenson 1995;Tschirren et al. 2013). A further explanation for low infection prevalence in juveniles was provided by Gomez-Chamorro (2019b), who demonstrated that, in a laboratory context, maternally transmitted antibodies may reduce the susceptibility of young bank voles to B.afzelli infection. The described mechanism may effectively limit the propagation of infection among young bank voles. Additionally, due to their higher dispersal rate (Kozakiewicz Table 3 The averaged confidence set of generalized linear mixed models (GLMMs; ΔAIC c ≤ 2, confidence set 0.95) explaining the spatial distribution of the Borrelia afzelii infection in 393 bank voles from 30 populations in NE Poland. Factors included in models: age, sex, frequency of TLR2 clusters or number of TLRC2 alleles (0, 1 or 2), winter conditions, mtDNA lineage, forest cover, sampling season. K the number of estimated parameters, AIC c Akaike's information criterion with a second-order correction for small sample sizes, ΔAIC c difference in AIC c between a given model and the most parsimonious model, ω i weight of the model  Perkins et al. (2003) showed that adult males of yellow-necked mouse (Apodemus flavicollis) carried more ticks than other combinations of sex and age class. In contrast to the results of Tälleklint and Jaenson (1995), we did not detect differences in Borrelia infection prevalence in bank voles captured during different times of the trapping seasons. However, we trapped voles only from July until early September, whereas the highest infection prevalence in bank voles has been observed in autumn (Tälleklint and Jaenson 1995). It is also possible that the 12.1% infection prevalence revealed in our bank vole populations was too low to track seasonal changes.
In contrast to previous studies (Karbowiak 2014;Medlock et al. 2013;Medlock and Leach 2015), our study did not find any relationship between B. affzelii prevalence in the bank vole and either winter conditions or forest cover. We may not have detected such associations, because we used averaged values of climate data representing a time span of 30-50 years, whereas most earlier studies have suggested year-toyear influence on the tick populations (Karbowiak 2014; Medlock et al. 2013;Medlock and Leach 2015). Moreover, any impact of climatic conditions on tick abundance might have been difficult to detect based on weather station data, because ticks could have actively avoided extreme temperatures by seeking microclimate refugia (Kilpatrick et al. 2017). Finally, the study by Hauser et al. (2018) found that climate variables that were measured at the study site were better predictors of tick abundance than climate variables measured at weather stations.
We tested the hypothesis that bank voles representing different clusters and genotypes of TLR2 differ in the prevalence of B. afzelii infection. However, we did not find any differences in infection prevalence between adult voles carrying different TLR2 genotypes. The lack of a significant association between infection prevalence in voles and TLR2 genotype frequencies supports the latest findings of Gomez-Chamorro et al. (2019a), whose results stand in opposition to the conclusions of other previously published studies (Tschirren et al. 2013). Under laboratory conditions, Gomez-Chamorro et al. (2019a) did not find any evidence that bank voles carrying different TLR2 genotypes (C1/C1, C2/C2, C3/ C3) varied in susceptibility to B. afzelii infection.
It is more likely that the most important factors influencing differences in infection prevalence among voles carrying different TLR2 clusters are the density of infected nymphs and the risk of exposure (Gomez-Chamorro et al. 2019a), these being potentially highly variable even on a small spatial scale (Vourc'h et al. 2016). Other mechanisms which may decrease B. afzelii infection prevalence in bank voles under natural conditions are grooming and acquired resistance to ticks Mysterud et al. 2019). In addition, infection with other parasites such as helminth worms may modify the immunological response to bacterial infection in rodents (Hübner et al. 2013). Thus, further studies on the immune genes conferring resistance to pathogens should not only distinguish variation in exposure from variation in susceptibility as indicated by Gomez-Chamorro et al. (2019a), but should also take into account parasite burden in the study organism.

B. afzelii infection detection
For the majority of bank voles analysed in our study, we used DNA extracted from a fragment of tail to detect infection with B. afzelii. It was possible to detect infection in an animal, when the concentration of bacterial DNA in a sample was at least 2 × 10 −4 ng/μl. Previous studies have shown differences in spirochete load between tissues and organs of voles: the smallest bacterial loads were observed in the heart and bladder, while the highest were observed in the skin and joints (Gomez-Chamorro et al. 2019a;Zhong et al. 2019). Our study showed that the tip of a tail (~4 mm, skin and tail bones) may also be used successfully for detection of B. afzelii infection.

Variability of the TLR2 gene and spatial distribution of TLR2 clusters
One of the aims of this study was to describe the genetic diversity of the TLR2 gene and to detect differences in the frequency of the TLR C2 cluster in populations of bank voles from NE Poland. Within the study area, over 80% of analysed sequences belonged to the most common TLR2 C1 cluster, while the other two clusters (TLR2 C2 and TLR2 C3 ) occurred in fewer than 10% of all analysed specimens. We found a discrepancy in TLR C2 cluster frequency between our study and an earlier study by Tschirren (2015). According to the Tschirren study, the frequency of the protective TLR C2 cluster in Poland was 29%, based on 12 individuals. Our results, based on a much larger data set (487 individuals), indicate that the frequency of this cluster is in fact much lower (7.2%). We found that the distribution of the TLR C2 cluster did not have a clear spatial structure, and varied from 0 to 25% among populations in NE Poland. Therefore, the study of Tschirren (2015) is likely to reflect only the TLR C2 cluster frequency of one location in Poland.
Of the three TLR2 clusters found in bank vole populations in NE Poland, TLR2 C3 was reported for the first time in this geographic region in the current study. This variant has already been described in Western Europe, with the highest recorded frequency in Spain (Morger et al. 2015). The genetic diversity of the TLR2 gene was high across all studied populations, similar to results obtained by Morger et al. (2015) for Central Europe, and those of Kloch et al. (2018) for two populations from NE Poland. Moreover, the number of private alleles was high, which indicates high population differentiation, despite the regional scale of the study. Our results are in concordance with previous research performed by Kloch et al. (2018). We showed that there was no spatial genetic structure of different TLR2 clusters at the regional scale, while Cornetti et al. (2018) showed clear spatial structure of TLR2 clusters, which depended on the risk of B. afzelii infection. We observed dominance of the TLR2 C1 cluster (in most populations its frequency was over 70%) and shifts in the frequency of the other two clusters (from 0 to 36%). Analyses performed on a continental scale have shown more distinct spatial patterns of TLR2 clustering (Morger et al. 2015;Tschirren 2015), but it was likely due to low sample size per country (< 20).
Correlation tests showed that there was no structuring of TLR2 clusters similar to genetic and spatial structure described for the mtDNA marker. A comparison of TLR2 genetic diversity with the variability of neutral markers in the same study area (10 microsatellite loci, Tarnowska et al. 2019), showed similar values in the mean number of alleles/ population, but not in the number of private alleles, which were higher for TLR2 (from 0 to 3, while the number of private alleles in microsatellite loci per population was from 0 to 0.6).

Parasite-mediated selection acting on the TLR2 gene
Parasite-mediated selection can be an important driver of evolution of the innate immune systems of hosts. Kloch et al. (2018) showed that some variants of the bank voles' TLR2 gene are associated with resistance to Babesia sp., which may reduce the fitness of voles and cause higher mortality, and to Heligmosomus mixitum. Similarly, in field vole (Microtus agrestis) polymorphism within TLR2 gene was associated with resistance to Cestode burden (Turner et al. 2011). Finally, Tschirren et al. (2013) indicated associations of TLR2 haplotypes with variation in susceptibility to B. afzelii.
However, in our study, we did not detect any association between TLR2 variability and differences in B. afzelii prevalence among voles, and we concluded that the impact of this pathogen in the studied populations was negligible, but we cannot exclude differences in virulence among B. afzelii. On one hand, our findings are supported by previous studies (Moody et al. 1994;Schwanz et al. 2011;Chambert et al. 2012;Voordouw et al. 2015;Cayol et al. 2018). For example, Moody et al. (1994) showed that only infants of the whitefooted mouse (Peromyscus leucopus) suffered from symptoms associated with B. burgdorferi infection, such as arthritis, carditis and neurological abnormalities. Capture-mark recapture studies in the field showed that Borrelia does not reduce survival of rodents and sea birds (Chambert et al. 2012;Voordouw et al. 2015). Experimental infections showed that Borrelia infection did not affect the health or behaviour of the white-footed mouse (Schwanz et al. 2011). Voordouw et al. (2015 point out that selection promotes B. burgdorferi genotypes which cause chronic, avirulent infection in their rodent host, because the transmission of bacteria depends on host mobility. On the other hand, B. afzelii infection did not influence bank voles' survival, but it had an impact on the reproductive success and spacing behaviour of bank voles (Cayol et al. 2018).

Different evolutionary history of mtDNA and TLR2 genes
We did not confirm our hypothesis that frequencies of TLR2 clusters differ between individuals belonging to two different mtDNA lineages of the bank vole in NE Poland. We detected the most common, deepest rooted (early arisen, considered as evolutionary older) haplotypes of TLR2 in both mtDNA lineages (see Fig. S2). Only the rarest haplotypes of TLR2, located at the edges of the TLR2 phylogenetic network (newly arisen) were associated with only one lineage. Similar to results presented by Morger et al. (2015), all three TLR2 clusters detected in North Eastern Poland were identified in bank vole individuals belonging to both the Carpathian and Eastern mtDNA lineages. Morger et al. (2015) and Tschirren (2015) did not find any relationship between the distribution of the mtDNA lineages and the TLR2 clusters in bank vole populations across Europe. Only one such association between mtDNA has been detected, between the Italian mtDNA lineage of bank vole and the TLR2 C4 (Morger et al. 2015). However, a recent study by Cornetti et al. (2016Cornetti et al. ( , 2018 showed the presence of the TLR2 C4 cluster within the range of the Western mtDNA lineage of bank voles. Therefore, our study supports the hypothesis of Morger et al. (2015) that TLR2 clusters are evolutionarily older than the divergence of mtDNA phylogenetic lineages in bank voles.

Conclusion
Question of the role of the TLR2 gene in maintaining immunity to B. afzelii is still open. In our study, we did not confirm that the C2 cluster of the TLR2 gene was involved in shaping resistance to B. affzelii infection. The mechanism of resistance to B. afzelii may be somehow associated with the TLR2 gene, but a deeper analysis of other immune genes, interactions between different parasites co-infecting a host, acquired immunity and behaviour of the host species is now required to reveal the exact role of this gene in immune response to Borrelia infection, as already proposed by Dizij and Kurtenbach (1995) and Mysterud et al. (2019). Our study indicates that B. afzelii, at least in NE Poland, is not a predominant driver of parasite-mediated selection in bank vole populations. under the National Strategy Reference Framework, which was cofinanced from the European Regional Development Fund of the EU).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval All procedures in this study involving animals were performed in accordance with the ethical standards of the Local Ethical Commission in Białystok, permissions nos. 07/2004, 15/2006, 16/2011, 43/2011, 12/2013 and 26/2014. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.