Human pathogen co-occurrence in Ixodes ricinus ticks: effects of landscape topography, 1 climatic factors and microbiota interactions 2

15 The factors shaping microbial communities within organisms are still poorly understood. Besides 16 ecological factors and host characteristics, direct interactions among microbes may shape the 17 occurrence of microbes and the structure of communities. In the past it has been difficult to 18 disentangle if patterns of microbial co-occurrence are due to facilitation or competition effects, or 19 shaped by shared ecological preferences (i.e., environmental filtering). Here we use a joint species 20 distribution model to characterize the bacterial microbiota composition of an important human 21 disease vector, the sheep tick Ixodes ricinus, along ecological gradients in the Swiss Alps, and to 22 test for facilitation or competition effects among human pathogens and tick endosymbionts. We 23 identify a number of ecological variables that significantly predicted the diversity of tick microbial 24 community and the occurrence of specific tick endosymbionts and human pathogens. However, 25 ecological associations were generally microbe-specific rather than universal. We also found 26 evidence for significant microbe interactions, in particular widespread facilitation among pathogens, 27 which promotes pathogen co-infection within ticks, as well as competition between the tick 28 endosymbiont Spiroplasma and a number of human pathogens. These findings highlight that direct 29 interactions among microbes can affect the vector competence of ticks and thereby tick-borne 30 disease dynamics. 31 large among-tick variation in microbiota composition, we found significant associations between specific environmental variables and OTUs. Particularly pronounced was the association

We quantified tick bacterial community composition by sequencing the hypervariable V3-V4 1 4 2 region of the 16S gene. Negative controls (N=5) were processed alongside the tick samples.   The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/559245 doi: bioRxiv preprint 1 9 4 pipeline [63]. Oligotyping uses all the sequences, which form an OTU, and performs Shannon 1 9 5 Entropy Analysis to regroup sequences based on within-OTU variation. This results in higher-1 9 6 . CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/559245 doi: bioRxiv preprint 1 0 resolution grouping than OTUs as the different oligotypes might differ only by a single nucleotide 1 9 7 [63]. We used the standard operation procedures of the oligotyping pipeline software [64].
Including a large number of explanatory variables in statistical models is inherently challenging. To 1 9 9 reduce the number of variables, while maintaining their information value, we used two variable specific OTU is an endosymbiont and/or a pathogen as traits [65]. This allowed us to test if OTUs. The set of explanatory variables under variable selection included additional information on 2 0 7 the environmental conditions of the sites (namely the number of days above 7 C° during the year, 2 0 8 monthly precipitation, mean monthly temperature, forest coverage, slope, aspect, bank vole 2 0 9 abundance, the proportion of voles to other rodents and expected tick heterozygosity) (Table S1).

1 0
We considered all parameter estimates, including associations among bacterial OTUs, having strong  Although JSDM is a powerful approach to model microbial community structure, it has a number of interactions [68]. Second, the model assumes that the explanatory variables affect the microbial 2 1 7 community composition (or rather, the presence or absence of individual OTU), but not vice versa. However, this is a valid assumption for most environmental (e.g. elevation and temperature) and tick-related variables (e.g. tick sex, life stage) included in our models. Thirdly, covariation among 2 2 0 variables poses a problem to any modelling approach. Our model is built on two distinct variable We observed a significant association between tick microbiota alpha diversity and elevation with 2 4 5 lower bacterial diversity observed at higher elevations (F 1,78 =3.98, p = 0.05, conditional R 2 = 0.07, 2 4 6 Figure 2). No other environmental variables were significantly associated with tick microbiota 2 4 7 alpha diversity (Table S2).

4 8
The analysis of tick microbiota beta diversity revealed a significant interaction effect between tick 2 4 9 life stages/sex and elevation, indicating that differences in tick microbiota beta diversity across tick 2 5 0 life stages/sex are shaped by elevation (pseudo-F 1 = 1.59, p = 0.01, R 2 = 0.03, Figure S1; Table S3). For other variables no significant association with microbial beta diversity was observed ( Table   2 5 2 S3).There was a significantly larger group dispersion at lower elevation sites (F 8,73 = 19.25, p < 2 5 3 0.001, Figure S2) suggesting that among-individual variation in bacterial community composition is 2 5 4 higher at lower elevation. For other variables, no significant heterogeneity in bacterial community 2 5 5 composition was observed (Table S4). Variance partitioning revealed that most of the variation in tick microbiota composition was 2 6 0 explained by tick ID: for the hundred most common OTUs, tick ID accounted for 83.5% of the 2 6 1 variation explained by the JSDM model. Fixed effects accounted for 9.3% and other random effects 2 6 2 (location, site and month) explained 7.2% ( Figure 3). This suggests that there is extensive among- individual variation which cannot be accounted for by environmental variation. Notably, the 2 6 4 situation differed for endosymbionts or pathogens: while tick ID was still the most important 2 6 5 variable, fixed effects explained 21.9% and random effects explained 11.2% of the total variation  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/559245 doi: bioRxiv preprint 1 3

Tick-specific and environmental factors associated with OTU occurrence
The occurrence of some OTUs was strongly predicted by specific explanatory variables (Table 3).

7 0
Slope, for example, was a strong predictor of the occurrence of a number of endosymbionts and 2 7 1 pathogens with a higher probability of R. helvetica, Rickettsia sp., Anaplasma and B. afzelii and a 2 7 2 lower probability of Ca. Neoehrlichia occurrence on steeper slopes. Females were less likely to 2 7 3 harbour the endosymbionts Spiroplasma, Lariskella and Rickettsia spp. (Table 3) and ticks at higher 2 7 4 elevations had higher probability to harbour R. helvetica and R. monacensis, but were less probable 2 7 5 to harbour B. garinii (Table 3). Ticks from sites facing northwards had a higher probability of harbouring Spiroplasma and B. afzelii, but a lower probability to harbor Anaplasma (Table 3).

7 7
Higher tick density was associated with a higher probability of Rickettsiella and Ca. Neoehrlichia  (Table 3), while more extensive forest coverage was associated with a higher probability 2 7 9 of R. helvetica and B. afzelii occurrence (Table 3). Associations between tick life stage, 2 8 0 temperature, the number of days > 7 C°, precipitation or relative vole abundance and the occurrence 2 8 1 of specific OTUs were not strongly statistically supported. The effect sizes of strongly statistically supported associations varied substantially ( Figure S5a-i). corresponded to a twofold increase in Neoehrlichia prevalence from 20% to 40% ( Figure S5e). Numerous bacterial OTUs were either significantly more or less likely to co-occur within a tick 2 8 8 than expected by chance after accounting for shared environmental preferences ( Figure 4; Table   2 8 9 S5). The occurrence of the tick endosymbiont Spiroplasma was negatively associated with the The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/559245 doi: bioRxiv preprint 1 4 they occurred, were all positive (Figure 4), suggesting that ticks are more likely to be co-infected 2 9 3 with several human pathogens simultaneously than expected by chance or based on shared 2 9 4 environmental preferences. Borrelia oligotypes showed positive co-occurrence patterns among each 2 9 5 other, except for B. miyamotoi, which was not associated with other Borrelia sp., but negatively 2 9 6 with Spiroplasma and positively with Lariskella. As our approach allowed us to control for shared 2 9 7 environmental preferences, the significant negative or positive associations among endosymbionts 2 9 8 and pathogens can be considered as indications of competition and facilitation. Using a Joint Species Distribution Modelling approach, our study provides evidence that a range of 3 0 2 tick endosymbionts and human pathogens are more or less likely to co-occur within I. ricinus ticks 3 0 3 than expected by chance or due to shared environmental preferences. Although the co-occurrence of 3 0 4 microbes within ticks has been documented before, both in I. ricinus [16,69] and other tick species filtering) and were therefore unable to disentangle shared responses to the environment from direct 3 0 7 microbe-microbe interactions. Our study reveals that when accounting for shared environmental suggesting that facilitation processes significantly contribute to endosymbiont-pathogen co-3 1 0 occurrence as well as pathogen co-infections within ticks. Such facilitation processes will affect the 3 1 1 structure and dynamics of microbial communities [9,71]. At the same time, the resulting co-3 1 2 occurrence of pathogens within ticks has implications for the severity, diagnosis, treatment and  Although associations among microbes were mostly positive, there were negative associations 3 1 5 between the tick endosymbiont Spiroplasma sp. and several human pathogens, which is indicative show pathogen-mediated host specialization as has been suggested previously [77,78]. In addition to testing for facilitation and competition among microbes, we used a JSDM framework to quantify the role of tick-specific and environmental variables in explaining spatial variation in 3 3 0 tick microbiota composition. Even though environmental variation across our sampling sites was and tick-related variables was considered in our models, overall they accounted for only 9.3% of the 3 3 3 variation in bacterial microbiota composition. This finding highlights that unexplained variation in 3 3 4 microbiota composition across individual ticks is substantial, which is in line with the idea that at further supported by the finding that the explanatory power of tick-related and environmental 3 3 8 variables was higher for tick endosymbionts and human pathogens, which are obligate resident.

9
Despite the large among-tick variation in microbiota composition, we found significant associations The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/559245 doi: bioRxiv preprint 1 6 between tick microbiota diversity and elevation. Across three independent elevational gradients we 3 4 2 found that alpha diversity was consistently higher at lower elevations. Although a decrease in plant 3 4 3 and animal diversity with increasing elevation is a widely documented pattern in ecology, it has 3 4 4 been suggested that the biogeographical patterns exhibited by bacteria may be fundamentally unlikely to affect human disease risk.

5 6
We identified a range of tick-specific and environmental variables that significantly predicted the 3 5 7 occurrence of various tick endosymbionts and human pathogens. However, these associations were 3 5 8 typically OTU specific rather than universal, with the same variable being both positively and was less likely to occur at higher elevations whereas R. helvetica and R. monascensis were more 3 6 1 likely to occur at higher elevations. Generally, the environmental factors shaping Rickettsia spp. distribution are poorly understood, as is their range of host species [44,81]. Yet, it has previously 3 6 3 been found that spotted fever incidence in humans, caused by R. ricketsii, is highest in suboptimal  contrast, the occurrence of the mammal specialist B. afzelii was not associated with elevation, 3 6 9 potentially because elevational clines in mammal diversity and/or abundance are less pronounced 3 7 0 [83]. Indeed we did not observe an association between elevation and bank vole abundance in our  Interestingly, temperature and precipitation, which vary strongly across elevational gradients clines. Indeed, slope and aspect, which are important determinants of the topography, and thus probability of pathogen and endosymbiont occurrence was higher on steeper slopes. Furthermore,  Host community composition and density are crucial for pathogen and endosymbiont occurrence 3 8 6 and dispersal [41]. Previous work has found that tick abundance is a strong predictor of Borrelia  No association between Borrelia spp occurrence and tick abundance was observed in our study.

8 9
However, both Ca. Neoehrlichia and Rickettsiella were more common at sites where ticks were 3 9 0 . CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/559245 doi: bioRxiv preprint 1 8 more abundant, suggesting that co-feeding transmission may also play a role in the life cycle of 3 9 1 these microbes.

9 2
Finally, differences in host competence can lead to dilution effects and thus affect the prevalence of proportion of bank voles to other rodents affects the prevalence of tick-borne pathogens. In conclusion, our results show how we can jointly model tick-associated bacterial communities, predict the occurrence of specific tick endosymbionts and human pathogens with strong statistical 4 0 7 support, but these effects were generally microbe-specific rather than universal. This highlights that  for Forest, Snow and Landscape Research WSL for providing spatial data and for help with spatial 4 2 6 analyses, and Otso Ovaskainen for discussions on the modelling.   The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/559245 doi: bioRxiv preprint