Bacterial microbiota composition of Ixodes ricinus ticks: the role of environmental variation, tick characteristics and microbial interactions

Introduction

Microbial communities within organisms consist of symbionts, commensals, mutualists and pathogens that co-occur simultaneously and potentially influence each other (Petney & Andrews, 1998; Rigaud, Perrot-Minnot & Brown, 2010; Sofonea, Alizon & Michalakis, 2015). These microbial communities may be shaped by a range of factors and processes, including the environment, host and microbe genetics and the occurrence and abundance of other microbial species (Adair & Douglas, 2017). For example, certain microbial species might tolerate only specific abiotic conditions, which makes it more likely that species with similar requirements co-occur within a host (‘environmental filtering’, Dallas & Presley, 2014). Similarly, the host’s immune system can influence colonization success of microbes (Hawley & Altizer, 2011), with cross-immunity preventing the colonization of different microbes with similar antigenic properties (Durand et al., 2015). Furthermore, mutualistic interactions between hosts and microbes can influence the structure of bacterial communities within host individuals (Chu & Mazmanian, 2013; Lee et al., 2013). Finally, direct interactions among microbes might affect colonization, or replication success after colonization, through competition or facilitation processes. Competition may occur when different microbes use the same, limited resources within a host (Lello et al., 2004), whereas facilitation may occur directly through the production of public goods (West & Buckling, 2003) or indirectly through the modification of the host’s physiology (Abraham et al., 2017) or immune defense (Rodríguez et al., 1999).

Ixodes ricinus is the most common tick species in Europe and an important vector for a range of human, domestic animal and wildlife pathogens (Medlock et al., 2013). Its distribution and abundance are strongly influenced by environmental conditions, in particular temperature and humidity (Cortinas et al., 2002; Gatewood et al., 2009). Previous studies that characterized the bacterial community composition of I. ricinus ticks have found mostly environmental and free-living bacteria but also several endosymbionts and human, domestic animal or wildlife pathogens, including Borrelia (Mannelli et al., 2012), Rickettsia (Venclikova et al., 2014), Anaplasma (Jahfari et al., 2014) and Candidatus Neoehrlichia (Kawahara et al., 2004).

Differences in the bacterial community structure and composition of ticks across habitats (Estrada-Peña et al., 2018), geographical sites (Carpi et al., 2011), and tick life stages and sexes (Carpi et al., 2011; Vayssier-Taussat et al., 2013) have been documented. Large-scale biotic or abiotic factors such as vegetation structure, elevation, temperature or rainfall may influence tick microbial communities directly, or indirectly through effects on tick physiology or activity patterns (Van Treuren et al., 2015) or via influencing the distribution and abundance of tick hosts species (Randolph et al., 1999; MacDonald et al., 2017). Small-scale and/or tick-specific effects on microbial communities may be explained by microhabitat or microclimatic conditions experienced by individual ticks (Gern, Morán Cadenas & Burri, 2008), individual tick behavior or genetics (Hawlena et al., 2013), direct biotic interactions among microbes (Moutailler et al., 2016) or parallel acquisition from a host during a bloodmeal (Andersson, Scherman & Råberg, 2014; Belli et al., 2017; Swei & Kwan, 2017).

Currently, the relative importance of these factors in shaping tick microbial communities is not well understood, which hampers progress in our understanding of the processes shaping microbial communities in nature and predicting the occurrence of specific microbes (e.g., human pathogens). Elevational gradients are excellently suited to quantify the importance of large-scale ecological variation in shaping tick bacterial microbiota because they cover a large range of environmental conditions within a small geographical area. Furthermore, including replicated transects along gradients allow us to quantify the robustness of ecological associations within sites and along elevational clines on tick microbial communities.

Ixodes ticks are commonly found to be co-infected with several (human, domestic animal and/or wildlife) pathogens (Andersson, Scherman & Råberg, 2013; Michelet et al., 2014; Diuk-Wasser, Vannier & Krause, 2016; Moutailler et al., 2016). Currently, it is unknown whether these co-infection patterns are caused by similar environmental preferences of pathogens, parallel acquisition from host communities or direct microbe-microbe interactions within ticks. Yet, previous studies suggest that the latter process, (i.e., facilitation and competition processes among microbes) may play a role in shaping microbial communities (Haine, 2008; Bonnet et al., 2017). For example, it has been found that pathogenic Rickettsia species prevent co-infection with other Rickettsia species in Dermacentor variabilis ticks (Macaluso et al., 2002), whereas the presence of Francisella sp. endosymbionts increases the colonization success of pathogenic Francisella novicida in D. andersoni ticks (Gall et al., 2016). Facilitation has also been suggested to promote co-infection with different Borrelia afzelii strains in Ixodes ricinus ticks (Andersson, Scherman & Råberg, 2013). Most strikingly, dysbiosis in I. scapularis ticks (i.e., ticks with low microbiotal diversity) leads to a defective peritrophic matrix which decreases the colonization success of B. burgdorferi s.s., suggesting that the pathogen requires the presence of an intact microbiota to be able to invade ticks (Narasimhan et al., 2014). Thus, the microbial community may have a crucial impact on vector competence of ticks and thereby on disease dynamics.

Yet, as outlined above, co-occurrence of microbes can be due to environmental filtering or direct microbial interactions, and distinguishing between these processes is non-trivial. Indeed, previous studies that have documented pathogen co-occurrence in ticks have not accounted for potential confounding variables such as shared ecological requirements, and are thus limited in their ability to differentiate between co-occurrences due to shared environmental niches, and co-occurrence shaped by facilitation or competition among microbes.

To address these gaps, we exploited the substantial environmental heterogeneity along replicated elevational gradients in the Swiss Alps to quantify the relative importance of environmental factors, tick characteristics and direct microbial interactions in influencing the structure of bacterial communities in I. ricinus ticks in general, and the (co-)occurrence of pathogens in particular, using a combination of 16S sequencing and joint species distribution modelling (JSDM) (Warton et al., 2015). By taking shared environmental preferences into account, JSDMs allows to identify residual co-occurrence patterns among microbes that can result from unaccounted environmental effects or direct microbial interactions. However, the correct spatial scale with regards to the focal biological processes is of importance, as well as the type of the hypothesized biotic interaction (Araújo & Rozenfeld, 2014; Zurell, Pollock & Thuiller, 2018) when interpreting JSDM patterns (Dormann et al., 2018).

Specifically, we ask (i) how do large-scale abiotic factors and small scale tick-level variables affect tick microbiota composition, (ii) which large-scale abiotic and small-scale tick-level variables predict pathogen occurrence, and (iii) are there patterns of non-random microbial co-occurrence that cannot be explained by environmental responses and might be due to unmeasured variables, such as microbial interactions.

Material and Methods

Tick sampling

Questing Ixodes ricinus ticks were collected at three locations in the Swiss Alps (Kanton Graubünden). At each location, one site at low (630–732 m above sea level, masl), one at medium (1,094–1,138 masl) and one at high (1,454–1 673 masl) elevation were identified (Fig. 1, Table 1, N = 9 sampling sites). At each site, questing ticks were sampled thrice, once in June, once in July, and once in August 2014 by dragging a white blanket (1 m × 1 m) over the ground vegetation as described previously (Lemoine, Cornetti & Tschirren, 2018). Ticks were collected from the blanket and stored in 95% ethanol. Tick species, sex and life stage were verified by morphological features following (Hillyard, 1996) using a stereomicroscope.

Table 1:

Tick sampling sites in the Swiss Alps.

Location	Site	Coordinates		Elevation	Succesfully sequenced Ixodes ricinus ticks
		North	East	(masl)	nymphs	males	females
1	Sagogn	46.783	9.233	693	0	9	15
	Flims	46.827	9.280	1138	3	5	3
	Ruschein	46.795	9.169	1454	0	1	1
2	Rodels	46.760	9.425	630	2	5	4
	Tomils	46.772	9.453	1144	3	6	4
	Feldis	46.789	9.453	1673	1	1	0
3	Passug	46.840	9.538	732	0	5	6
	Castiel	46.826	9.569	1094	0	3	3
	Praden	46.817	9.589	1582	1	0	1

DOI: 10.7717/peerj.8217/table-1

Results

Tick microbiota variance partitioning

Variance partitioning revealed that most of the variation in tick microbiota composition explained by our model related to tick ID: for the hundred most common OTUs, tick ID accounted for 64.1% of the variation explained by the model. Fixed effects (e.g., tick life stage, elevation; see Table S1) accounted for 12.5% (tick-level: 7.3%, site-level: 5.2%) and spatial and temporal random effects (i.e., location, site and month) explained 23.3% (Fig. 2). This suggests that there is extensive tick-level variation which cannot be accounted for by tick-specific characteristics included in our model (i.e., sex, life stage, genetic diversity) or site-level environmental factors. The pattern differed slightly for endosymbionts and human, domestic animal or wildlife pathogens: while tick ID was still the most important variable explaining 49.9%, fixed effects explained 31.8% (tick-level: 17.5%, site-level: 14.3%) and spatial and temporal random effects explained 18.3% of the total variation explained by the model, when averaged over all pathogens and endosymbionts (Fig. 2). Thus, tick- and site-level fixed effects explained a larger proportion of the variation in the occurrence of obligate resident pathogens and endosymbionts than the occurrence of other (potentially non-resident) OTUs.

Patterns of microbial association within ticks

Numerous bacterial OTUs were either significantly more or less likely to co-occur within a tick than expected by chance after accounting for shared environmental preferences (Fig. 3A; Table S2). At the level of the individual tick, the occurrence of the tick endosymbiont Spiroplasma was negatively associated with the occurrence of the endosymbiont Lariskella and several tick-borne pathogens, namely Rickettsia sp., Ca. Neoehrlichia and B. miyamotoi (Fig. 3A). Associations among pathogens, if they occurred, were all positive (Fig. 3A), suggesting that ticks are more likely to be co-infected with several human, domestic animal or wildlife pathogens simultaneously than expected by chance or based on shared environmental preferences. Borrelia oligotypes showed positive co-occurrence patterns among each other, except for B. miyamotoi, which was not associated with other Borrelia sp., but negatively with Spiroplasma and positively with Lariskella. At the level of the sampling site, significant associations were sparser. Spiroplasma was more likely to co-occur with Lariskella and Rickettsiella across sites, whereas the latter two were less likely to co-occur across sites than expected by chance after accounting for shared environmental preferences (Fig. 3B). At the level of month or location, there were no significant associations.

Discussion

We used a JSDM framework to quantify the relative importance of large scale, site-level environmental variables, tick-level characteristics and interactions among microbes in shaping tick microbiota composition along elevational gradients in the Swiss Alps. We show that although our study captured extensive environmental variation, with sampling sites spanning across an elevational gradient from 630–1,580 masl, and a large number of ecological variables was considered in our models, the explanatory power of such large-scale ecological factors was comparably weak. In contrast, individual tick ID explained over 60% of the variation in microbiota composition. This substantial microbiota variation across individual ticks may be partly explained by some of the bacteria present in ticks being non-resident (i.e., bacteria that were by chance obtained from the environment through the mouth, the anal pore or spiracles or during blood-feeding; Horner-Devine & Bohannan, 2006; Engel & Moran, 2013; Zolnik et al., 2016; Zolnik et al., 2018; Ross et al., 2018). Indeed, there has been a debate whether ticks have a stable microbiota (Ross et al., 2018), mirroring the wider debate on how common resident microbiota is in arthropod hosts (Hammer et al., 2017).

However, also for endosymbionts and human, domestic animal or wildlife pathogens, which are obligate resident, tick ID accounted for half of the variation in occurrence, suggesting that microhabitat or -climatic conditions experienced by individual ticks, tick-specific traits and behaviors not included in our models, as well as microbial interactions within ticks such as facilitation and competition (Abraham et al., 2017; Gurfield et al., 2017), play a crucial role in shaping microbiota composition and the occurrence of endosymbionts and human or wildlife pathogens in I. ricinus. Focusing on such small-scale variables, rather than large-scale climatic or environmental factors as is usually done when modelling tick-borne pathogen prevalence (Norman, Worton & Gilbert, 2016; Rosà et al., 2018), is thus likely a more fruitful approach to advance our understanding of microbiota composition of natural populations as well as (tick-borne) disease dynamics.

Co-occurrence of human, domestic animal or wildlife pathogens in ticks has been documented previously, both in I. ricinus (Lommano et al., 2012; Michelet et al., 2014) and other tick species (Gurfield et al., 2017; Laaksonen et al., 2018). Yet, previous studies did not control for environmental filtering, which limited their ability to disentangle shared responses to the environment from direct microbe-microbe interactions. Our study revealed that when accounting for shared environmental preferences, associations among human or wildlife pathogens were often pronounced and mostly positive. These positive associations may result from direct facilitation among microbes or parallel colonization from co-infected tick hosts. Because our sampling unit was the whole tick, whereas bacteria inhabiting a tick can be situated in different organs, co-occurrence at the tick-level does not necessarily mean that there is direct interaction between co-occurring OTUs, although indirect interactions, via, e.g., host immune system, can still occur.

Within ticks, the significant positive associations among the Lyme disease-causing Borrelia genospecies (B. afzelii, B. garinii and B. valaisiana) were particularly striking. This positive co-occurrence is surprising because B. garinii and B. valaisiana are bird specialists (Hanincova et al., 2003b; Comstedt, Jakobsson & Bergström, 2011), whereas B. afzelii is a rodent specialist (Hanincova et al., 2003a). Thus, the parallel colonization from co-infected tick hosts cannot explain this pattern. Rather the positive co-occurrence is indicative of facilitation processes among Borrelia genospecies, as has been suggested previously (Andersson, Scherman & Råberg, 2013). Such facilitation, and the resulting co-infection of ticks with several Borrelia genospecies has implications for the severity, diagnosis, treatment and control of Lyme disease. Finally, the co-occurrence of these different Borrelia genospecies suggests that I. ricinus feeds on multiple, phylogenetically diverse host species during its life cycle and does not show pathogen-mediated host specialization as has been suggested previously (McCoy et al., 2005; McCoy, Léger & Dietrich, 2013).

Although associations among microbes were mostly positive, there were negative associations between the tick endosymbiont Spiroplasma and several human or wildlife pathogens, which may be explained by competition. The most common infection route for Spiroplasma is maternal (i.e., vertical) transmission (Herren & Lemaitre, 2011), indicating that horizontal or environmental transfer plays a minor role in its transmission. Protective effects of Spiroplasma have been previously described in Drosophila spp., where Spiroplasma is associated with a decreased probability of nematode and parasitoids infections (Xie, Vilchez & Mateos, 2010; Jaenike et al., 2013). Although the exact mechanisms mediating Spiroplasma-induced competition effects are currently unknown, this finding may stimulate further research into the potential of tick endosymbionts to manage tick-borne pathogens.

In contrast to the numerous positive or negative associations among microbes at the tick-level, little statistical support for positive or negative microbial co-occurrence was found at the site-level, with the exception of the associations among three endosymbionts. Interestingly, the pattern of co-occurrence of Spiroplasma and Lariskella at the site-level was opposite from what was observed at the tick-level. It suggests that Spiroplasma and Lariskella are more likely to co-occur at the same sites but not within the same tick. Co-occurrence at the site-level can be due to an environmental variable not included in our model, for which the three OTUs had similar responses. It has also been suggested that negative associations generate checkerboard patterns of co-occurrence that can be captured at finer spatial scales but that are lost with increasing scales, but positive associations can be captured across scales (Araújo & Rozenfeld, 2014).

Despite the large among-tick variation in microbiota composition, we identified a range of environmental variables that significantly predicted the occurrence of specific tick endosymbionts and human, domestic animal or wildlife pathogens. However, the predictor variables as well as their effect were typically OTU-specific rather than universal. For example, B. garinii was less likely to occur at higher elevations, whereas R. helvetica and R. monascensis were more likely to occur at higher elevations. Generally, the environmental factors shaping Rickettsia spp. distribution are poorly understood, as is their range of host species (Halos et al., 2010; Eremeeva & Dasch, 2015). Yet, it has previously been found that spotted fever incidence in humans, caused by R. ricketsii, is highest in areas or regions, where ticks are less common (Atkinson et al., 2013). This is in line with our findings and suggests that Rickettsia spp. are more likely to colonize ticks living under suboptimal conditions (e.g., at range edges).

The finding that B. garinii is less likely to occur at higher elevations is in line with previous observations (Jouda, Perret & Gern, 2004b; Cornetti et al., 2018) and may be explained by changes in vegetation structure and associated changes in host communities (Halos et al., 2010), in particular changes in the diversity and/or abundance of birds, the natural hosts of B. garinii (Comstedt, Jakobsson & Bergström, 2011). In contrast, the occurrence of the mammal specialist B. afzelii was not related to elevation, potentially because elevational clines in mammal diversity and/or abundance are less pronounced (McCain, 2005). Indeed, we did not observe an association between elevation and bank vole abundance across our study sites (ANOVA: F_1,8 = 0.357, p = 0.57, R² = 0.05).

Interestingly, temperature and precipitation, which vary strongly across elevational gradients (average temperature and precipitation: high elevation sites: 11.8 °C and 17.8 mm per month; in low sites: 16.5 °C and 12.1 mm per month), were not significant predictors of the occurrence of endosymbionts or human or wildlife pathogens, with the exception of precipitation correlating positively with the probability of Rickettsiella occurrence. This may be partly explained by the temperature and precipitation measures included in our models not fully capturing the microclimatic variation across sites and along elevational clines. Indeed, slope and aspect, which are important determinants of the topography, and thus microclimate (Bennie et al., 2008), were significant predictors of pathogen and endosymbiont occurrence. The probability of Rickettsia sp. occurrence was higher on steeper slopes. Furthermore, the probability of occurrence was higher on north-facing slopes for B. afzelii and Spiroplasma and higher on south-facing slopes for Rickettsiella (see also Stuen, Granquist & Silaghi, 2013). Microclimatic conditions may affect microbial occurrence directly, or indirectly via affecting tick behavior or host community composition (Swei, Meentemeyer & Briggs, 2011; Lawson et al., 2014). Furthermore, topography can affect population connectivity and dispersal in metapopulation networks (Swei & Kwan, 2017).

Previous work has found that tick abundance is a strong predictor of Borrelia spp. prevalence, potentially because larger tick populations facilitate co-feeding transmission (Jouda, Perret & Gern, 2004a). No relationship between Borrelia spp. occurrence and tick abundance was observed in our study. However, both Ca. Neoehrlichia and Rickettsiella were more common at sites where ticks were more abundant, suggesting that co-feeding transmission may also play a role in the life cycle of these microbes.

Finally, differences in host competence can lead to dilution effects and thus affect the prevalence of tick-borne pathogens (Keesing, Holt & Ostfeld, 2006). Whereas for some tick-borne pathogens the vertebrate hosts are known or suspected (e.g., small mammals for B. afzelii (Hanincova et al., 2003a) and Ca. Neoehrlichia (Jahfari et al., 2012), birds for B. garinii and B. valaisiana (Hanincova et al., 2003b), both for Anaplasma (Keesing et al., 2012) and R. helvetica (Sprong et al., 2009), for others the host species range is less well understood (e.g., B. miyamotoi; Wagemakers et al., 2015). The bank vole is a common tick host at our study sites and their abundance was a significant negative predictor of R. monacensis and R. helvetica occurrence. Interestingly, bank voles are not known hosts for either (Burri et al., 2014). Most likely, the relation is thus indirect, explained by an unmeasured biotic or abiotic variable that correlates with bank vole abundance. No evidence was found that the proportion of bank voles to other rodents affects the prevalence of tick-borne pathogens.

A limitation of our sampling design is the uneven sample distribution across sites. We collected ticks up to the upper elevational limit of tick distribution, which leads to a large variation in environmental variables included in our models, but at the same times means that we have a limited number of samples from the high elevation sites. Yet, adequate model fit suggests that this uneven sample distribution did not compromise model performance. Furthermore, although JSDM is a powerful approach to model community structure, it has a number of limitations. First, it assumes that interactions among microbes are similar across environments (but see Tikhonov et al., 2017). This is not necessarily the case as both environmental factors and tick host community may shape microbial interactions (Elliot, Blanford & Thomas, 2002). Second, the model assumes that the explanatory variables affect the microbial community composition (or rather, the presence or absence of individual OTUs), 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 explanatory variables poses a problem to any correlative modelling approach. Our model is built on two distinct variable sets to aid in handling such covariation: the full variable set includes elevation, whereas the variables with the strongest covariation (i.e., temperature and precipitation) are included in the variable selection set. Fourthly, the inferred residual associations between focal taxa are assumed to be symmetrical. If there are asymmetric interactions (e.g., predator–prey-relationships), the sum outcome can be seen as either positive or negative correlation (Zurell, Pollock & Thuiller, 2018). However, in our study, the expectation was facilitation or competition, which are symmetric positive or negative interactions, respectively. Thus, given sufficient signal, we expect that the focal interactions can be captured by our modelling approach.

Conclusions

Our study demonstrates that a JSDM framework can contribute to a better understanding of the factors shaping bacterial communities in natural populations as well as patterns of co-occurrence among microbes. Overall, our study highlights the role of small-scale, tick-level characteristics rather than large-scale ecological variation in shaping microbial communities of I. ricinus. We identified a number of ecological variables that predict the occurrence of specific tick endosymbionts and human, domestic animal or wildlife pathogens with strong statistical support, but these effects were typically microbe-specific rather than universal. This highlights that environmental change can have different, even opposite effects on different human pathogens, and thus disease risk. Furthermore, by accounting for shared environmental preferences, our approach identified patterns of microbial co-occurrence that are consistent with microbe-microbe interactions, which result in pathogen co-infections within ticks, as well as competition between Spiroplasma and a number of human, domestic animal or wildlife pathogens. The latter opens up new and exciting avenues for the control and management of tick-borne diseases in regions with high human disease incidence.

Supplemental Information