Elevational Changes in Bacterial Microbiota Structure and Diversity in an Arthropod-Disease Vector

Environmental conditions change rapidly along elevational gradients and have been found to affect community composition in macroscopic taxa, with lower diversity typically observed at higher elevations. In contrast, microbial community responses to elevation are still poorly understood. Specifically, the effects of elevation on vector-associated microbiota have not been studied to date, even though the within-vector microbial community is known to influence vector competence for a range of zoonotic pathogens. Here we characterize the structure and diversity of the bacterial microbiota in an important zoonotic disease vector, the sheep tick Ixodes ricinus, along replicated elevational gradient (630–1673 m) in the Swiss Alps. 16S rRNA sequencing of the whole within-tick bacterial microbiota of questing nymphs and adults revealed a decrease in Faith’s phylogenetic microbial alpha diversity with increasing elevation, while beta diversity analyses revealed a lower variation in microbial community composition at higher elevations. We also found a higher microbial diversity later in the season and significant differences in microbial diversity among tick life stages and sexes, with lowest microbial alpha diversity observed in adult females. No associations between tick genetic diversity and bacterial diversity were observed. Our study demonstrates systematic changes in tick bacterial microbiota diversity along elevational gradients. The observed patterns mirror diversity changes along elevational gradients typically observed in macroscopic taxa, and they highlight the key role of environmental factors in shaping within-host microbial communities in ectotherms. Supplementary Information The online version contains supplementary material available at 10.1007/s00248-021-01879-5.


Introduction
Because of the rapid changes in abiotic conditions along elevational gradients, radical elevational species community turnovers are often observed within short geographical distances [1]. Indeed, temperature and oxygen pressure decrease, whereas ultraviolet radiation increases with increasing elevation [2]. Other abiotic factors, such as precipitation, wind velocity, seasonality, soil formation processes, and disturbance, typically also show systematic changes along elevation clines, but also vary by geographical context [2]. All of these abiotic factors affect species community composition and diversity along elevational gradients [1].
In macroscopic taxa, such as plants, invertebrates, and mammals, biodiversity is typically found to either linearly decrease with increasing elevation or to show a humpshaped pattern where diversity is highest at intermediate elevations [3]. In addition, community composition and the strength and direction of biotic interactions have been found to vary along elevational clines [4,5]. Much less is known about how microbial communities change along elevational gradients, and the existing empirical studies suggest inconsistent elevation patterns that differ from patterns observed in macroscopic taxa [e.g., [7][8][9].
If diversity gradients along elevational clines are different in microbes compared to macroscopic taxa, the underlying factors affecting these gradients are likely to differ as well.
For macroscopic taxa, climatic variables seem to be the most important factors affecting diversity patterns [3], whereas for microbes, we do not yet have a good framework to understand trends in community composition [9]. To date, the best studied microbial communities are soil bacteria: their diversity seems to be mainly determined by the quality and composition of the soil, such as soil pH or carbon, without any systematic changes along elevational clines [6].
Compared to 'free-living' microbial communities, hostassociated microbiota are shaped by an additional key factor: the identity, ecology, life history, and quality of the host [10], which may affect microbial community composition along elevational gradients. However, patterns of diversity and community composition in host-associated microbiota do not seem to follow a consistent pattern either. For example, in pika (Ochotona curzionae) individuals living at higher elevation were found to have a higher alpha (i.e., the mean species diversity) and beta (i.e., heterogeneity in species composition) diversity in their gut microbial community compared to individuals living at lower elevation [11], whereas human skin microbiota shows a decrease in alpha diversity but an increase in beta diversity with increasing elevation [12].
Microbial communities of ectothermic hosts are expected to be most strongly affected by elevational gradients because ectotherms do not buffer ambient temperature as strongly as endothermic species. Furthermore, the microbiota of invertebrates is typically not stable [13]. Yet, changes in microbiota composition in invertebrate hosts along environmental gradients remain largely unexplored.
Ixodes ticks are common vectors for human pathogens including Borrelia sp. and Rickettsia sp., and recent studies suggest that the tick's commensal microbiota composition affects the probability of harboring human pathogens [14]. Ixodes ricinus has three life stages throughout its life cycle: larva and nymph stages, during which a blood meal is needed to progress to the next stage, and the adult stage, during which females require a blood meal to lay eggs. Usually the hosts for earlier life stages are small-sized mammals or birds, whereas adult ticks seek blood meals from larger mammals, such as ungulates. Due to climate warming, I. ricinus has expanded its distribution into both higher latitudes and higher elevations in many parts of Europe [15]. At the same time, the incidence of the diseases caused by tick-borne pathogens, such as Borrelia sp., continues to increase in many regions [16]. The abundance and distribution of ticks are strongly influenced by temperature and other climatic variables: low winter temperatures increase tick mortality, whereas warmer temperatures during summer months lead to a faster life cycle and a longer activity period [17,18].
While there is growing body of research on how changing environmental conditions might affect specific tick-associated microbes, specifically pathogens [19][20][21], there is a lack of studies on how the structure of the commensal tick microbiota (or the microbiota of any other disease vector) changes along elevational clines and how this may indirectly affect pathogen prevalence and disease dynamics (but see [21]). Abiotic variables may affect tick microbiota composition and diversity either directly through effects on microbial growth, competition, and/or transmission [22] or indirectly through changed tick behavior or life history [17]. Furthermore, ticks quest in the undergrowth, attach to their vertebrate host, and suck blood for a number of days. The tick microbiota is thus likely acquired from soil and plants but also from their host's skin and blood [23], which all likely vary along elevational clines and may thus shape tick microbial communities.
Here we exploit the rapidly changing environmental conditions along elevational gradients in the Swiss Alps to quantify changes in tick microbiota diversity and community structure along elevational cline. Specifically, we test if tick bacterial microbiota diversity and community structure vary along elevational gradient and if these patterns differ with season or across tick life stages and sexes. Furthermore, we analyzed how ecological processes influence community turnover by comparing phylogenetic relatedness within and between tick microbiota.
Based on previous findings, we predict systematic differences in tick bacterial microbiota composition along elevational gradients. Furthermore, because of differences in their behavior and ecology [17,18], we predict differences in microbiota structure and diversity between tick life stages and sexes.

Tick Sampling
We collected questing Ixodes ricinus ticks at three different locations in the Swiss Alps (Kanton Graubünden). Three sites per location were identified, one at low (630-732 m above sea level), one at medium (1094-1138 m) and one at high (1454-1673 m) elevation ( Fig. 1; Table 1). At each site, tick sampling was performed 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 [24] (see Supplementary Materials S1 for additional details).

Tick Microbiota
16S rRNA sequencing has been described previously [21]. In short, we randomly selected I. ricinus ticks from each sampling site (Table 1). We cut ticks in half with a sterilized blade to facilitate DNA isolation and then used the DNeasy Blood & Tissue kit (Qiagen; Hilden, Germany) to extract DNA. We processed negative controls (N = 5) alongside the tick samples. We characterized tick bacterial community composition by sequencing the hypervariable V4 region of the 16S gene using the primers 515FB and 806RB [25] and prepared sequencing libraries following the Earth Microbiome 16S Illumina Amplicon protocol (see Supplementary materials S1 for details). The libraries were sequenced on Illumina MiSeq at the Functional Genomic Center Zurich with a target length of 250 bp following the manufacturer's protocol.
Sequences were analyzed using the mothur pipeline with MiSeq standard operation procedures [26]  Using mothur, we purged unsuccessful contigs and preserved only contigs between 250 and 310 bp. The alignment was made against aligned SILVA bacterial references  (release 128; https:// www. arb-silva. de/ docum entat ion/ relea se-128/). We used 99% similarity to determine OTUs and classified them using SILVA taxonomy. Only samples with > 500 amplicons and OTUs which were present in at least two samples were included in the analyses. We rarified the samples to 500 amplicons to account for variation in amplicon numbers. The rarefaction approach was required because of substantial variation in amplicon numbers across samples.

Bacterial Alpha Diversity
Bacterial taxonomic alpha diversity (inverse Simpson index; [27]) and Faith's phylogenetic alpha diversity [28] were calculated with the R package vegan [29]. Additionally, we calculated two more phylogenetic alpha diversity indices: Nearest Relatedness Index, NRI, (equivalent to -1 standardized effect size of mean pairwise distances in communities, which estimates the average phylogenetic relatedness between all possible pairs of bacterial taxa within a tick) and Nearest Taxon Index, NTI, (− 1 times standardized effect size of mean nearest taxon distances in communities, which calculates the mean nearest phylogenetic neighbor among the bacterial taxa within a tick). Thus, while NTI reflects the phylogenetic structuring near the tips of the tree, NRI reflects structuring across the whole tree. The ratio between these two measures (i.e., NTI/NRI) provides a measure of phylogenetic clustering among OTUs: if NTI/NRI is positive, it suggests that there is phylogenetic clustering of OTUs (i.e., closely related OTUs are more likely to co-occur than by chance), whereas negative values indicate phylogenetic overdispersion (i.e., co-occurring OTUs are less related than expected by chance) [30]. We performed the analyses using the picante package [31]. To create null models, we randomized the bacterial community compositions obtained from the data using same distance matrices, but randomizing the bacterial OTU labels across taxa.
For each alpha diversity measures we used linear mixed models with the R package lme4 [32] to test for associations between bacterial alpha diversity and tick life stage/sex, sampling month, and linear and quadratic terms of elevation with full interactions between linear terms. Sampling location was included as a random effect in the model. We used a model selection approach based on Akaike's Information Criterion and model fit with conditional R 2 in the package piecewiseSEM [33] to test which combination of factors best describes variation in tick bacterial alpha diversity. The model selection is presented in the supplementary methods (S2), while the final model is described in the results.

Bacterial Beta Taxonomic Diversity
We analyzed tick bacterial beta diversity on pairwise matrices using five different indices. Two of the indices measure taxonomic beta diversity: Bray-Curtis dissimilarity, which takes into account the abundance of OTUs [34] and Jaccard index, which takes only presence-absence of OTUs into account [35], thus providing information on both aspects of beta diversity. The other three indices measure phylogenetic beta diversity: weighted UniFrac (wUF), which takes into account the unshared branch lengths for all OTUs and weights OTUs based on OTU counts [36], equivalent to the phylogenetic alpha diversity index; βNTI, which is a between-community equivalent of NTI (see above) [37], and Bray-Curtis-based Raup-Crick (BC-RC) abundance which measures the deviance of observed turnover while taking into account OTU relative abundances (i.e., between-community equivalent of NRI) [38]. Bray-Curtis and Jaccard indices were calculated with the R package vegan, UniFrac with mothur, βNTI with MicEco package [39], and BC-RC with code presented by Stegen et al. [38].
First, we performed permutational ANOVA with dissimilarity matrices using the package vegan to test for association between the first three measures of beta diversity (Jaccard, Bray-Curtis, UniFrac) and elevation, tick life stage/ sex, sampling location, and sampling month. Permutational ANOVA partitions distance matrices among sources of variation and fit linear models [40]. Initial models included all variables and interactions and, if non-significant, they were dropped during the model selection by removing the variables with the highest p-value, starting with the least significant interaction. Model selection was evaluated based on R 2 values of remaining variables. The results of the permutational ANOVA are visualized by performing non-metrical multidimensional scaling on Bray-Curtis dissimilarities and plotting the samples on the two first axes.
Second, we used an analysis of multivariate homogeneity of group dispersion using the package vegan to test whether elevation, sampling location, sampling month, or tick life stage/sex is associated with tick microbiota composition, again using the first three measures of beta diversity (Jaccard, Bray-Curtis, UniFrac). This analysis is a multivariate analogue of Levene's test for homogeneity of variances [41] and tests whether variation in community composition among groups is similar.

Influences of Ecological Processes on Community Turnover
We used two measures of beta diversity (βNTI and BC-RC) for the analysis of ecological processes on community turnover. To study which ecological processes shape within-tick bacterial community composition, we used the phylogenetic signal of organismal niches as described by Stegen et al. [38, Elevational Changes in Bacterial Microbiota Structure and Diversity in an Arthropod-Disease… 871 42]. By assuming that closely related taxa are ecologically more similar to each other and thus their niches are more similar, we can infer which processes govern community composition. Stochastic dynamics should lead to random community assembly, environmental filtering should lead to a community consisting of taxa that are more closely related than expected by chance, whereas strong competition should lead to a community consisting of less closely related taxa. Finally, environmental change should lead to increased phylogenetic turnover. Two cases of deterministic processes are possible: if βNTI < − 2, phylogenetic turnover is lower than expected by chance suggesting consistent selective pressures (homogenous selection), if βNTI > 2, phylogenetic turnover is higher than expected by chance, suggesting shifts in selective pressure due to environmental change (variable selection). If βNTI is between − 2 and 2, it suggests stochastic processes determine community composition [42]. If BC-RC < − 0.95, the compositional turnover between communities is low, thus suggesting a strong dispersal between two communities (homogenizing dispersal). If BC-RC > 0.95, turnover is high due to a low rate of dispersal leading to ecological drift (dispersal limitation). Finally, in situations of moderate dispersal and weak selection, it is possible that none of these four processes shape community composition (undominated) [42]. We analyzed processes within sites and compared these within-site results for tick sex/stages and elevations. Due to low sample sizes per site, nymphs and samples from high elevations were not included in community turnover analysis.
Additional analyses of phylosymbiosis between tick population genetic structure and microbiota structure and random forest classification are presented in the Supplementary material (S3 and S4). Furthermore, we present sensitivity analyses in the Supplementary material (S6), which confirm the robustness of our findings to unequal sample sizes across elevations and life stages.

Results
We sequenced the microbiota of 92 Ixodes ricinus ticks and five negative controls, resulting in 13 214 477 amplicons. No amplification was observed in the negative controls. After contig assembly and quality control, 1 802 719 sequences were retained for OTU analysis. There was a median of 1 661 quality-controlled amplicons per tick, with an interquartile range of 5 744. 79 samples with more than 500 amplicons per sample and a Good's coverage estimator ≥ 0.95 was included in the diversity analyses (Fig. 2).
In total, 5 181 bacterial OTUs were identified. The median number of OTUs per rarified sample was 83 OTUs, with a 95% confidence interval of 30-121 OTUs. After excluding OTUs that occurred in only one sample, 864 OTUs were used in subsequent analyses. Four OTUs were present in at least 90% of the samples and represented 38.9% of all amplicons: Candidatus Midichloria (Otu0001) (which was present in all samples), Pseudomonas (Otu0002), and Sphingomonas (Otu0006 and Otu0009).

Ixodes Ricinus Microbiota Alpha Diversity
There was a significant elevation and sampling month effect on tick microbiota alpha diversity based on Faith's phylogenetic index. Lower bacterial diversity was observed at higher elevations (2.0 index points per 1000 m) and diversity increased from June  Fig. 3). No other environmental variables were significantly associated with tick microbiota alpha diversity (Table S1)

Ixodes ricinus microbiota beta diversity
First, the analysis of tick microbiota beta diversity based on Jaccard index revealed significant differences in microbiota compositions along elevational clines. In addition, tick stage/ sex was a significant predictor of beta diversity across all beta diversity indices (Table 3). No other variable was significantly associated with microbial beta diversity (Table S2).  Second, a significantly larger group dispersion was observed at lower elevations (Bray-Curtis dissimilarity: F 8,70 = 5.9, adj. p < 0.001; Jaccard distance: F 8,70 = 11.4, adj. p < 0.001; wUF: F 8,70 = 3.1, adj. p = 0.02, Fig. 4), suggesting that among-tick variation in bacterial community composition is higher at lower elevations. Additionally, significant group dispersion was observed across sampling locations (Bray-Curtis index F 2,76 = 5.3, adj. p = 0.02), while for other indices and variables, no significant heterogeneity in bacterial community composition was observed (Table S3). Sensitivity analyses demonstrated that the results of the diversity analyses were not biased by unequal sample sizes across elevations or life stages (Supplementary material S6).

Influences of Ecological Processes on Community Turnover
In general, Undominated ecological processes were the most common relationships among communities (Fig. 5, Supplementary Materials Table S5-6). It suggests a moderate rate of dispersal among communities and relatively  weak selection. There were no statistically significant differences in community turnover process between samples from the same vs. different sites (χ 2 4 = 2.45, p = 0.65), between low and medium elevation sites (χ 2 4 = 4.10, p = 0.39), and between females and males (χ 2 4 = 6.88, p = 0.14). As between-site and within-site processes had similar distributions, the scale of the processes affecting tick microbiota is likely either larger (i.e., geographical) or smaller (i.e., within-tick) than our study. Among the groups with larger sample sizes, it is notable that female ticks showed a very low proportion of variable selection (1.1%), whereas variable selection was of substantially higher importance in male ticks (6.3%). In contrast, homogenous selection was more pronounced in female ticks (7.3%) compared to male ticks (3.3%). The sum of deterministic processes (i.e., homogenous and variable selection combined) was similar in females and males (8.4% and 9.5%, respectively). It suggests that deterministic processes have a similar, although limited, effect in shaping the microbiota composition of female and male ticks, while the specific type of selection (homogenous vs. variable) differs between sexes. Stochastic processes (i.e., dispersal limitation, homogenizing dispersal, and especially undominated processes) were found to shape tick bacterial community composition in both sexes (females, 84.8%, males, 82.8%).

Discussion
We observed significant changes in bacterial alpha and beta diversity in Ixodes ricinus ticks along replicated elevational gradients in the Swiss Alps. Alpha diversity measured as Faith's phylogenetic distance significantly decreased with increasing elevation, mirroring elevational diversity patterns observed in many macroscopic taxa [3], as well as in human skin microbiota [12]. In addition, microbiota composition (beta diversity) measured as Jaccard's index differed along elevational clines with significantly lower variation in bacterial community composition at higher elevations. These differences in tick bacterial beta diversity along elevational clines contrast patterns observed in previous studies that found either no association with elevation (e.g., in soil bacteria [6]) or higher beta diversity at higher elevations (e.g., in mammalian skin and gut microbiota [11,12]). We also observed seasonal changes in tick microbiota composition, with higher diversity later in the season.
Elevational changes in diversity were observed for some of the diversity indices, but not for others. Faith's phylogenetic diversity measures the sum of the branch lengths of a phylogenetic tree. In contrast to taxonomic indices of alpha diversity (e.g., inverse Simpson's index), it does not consider OTU abundances. Similarly, Jaccard beta diversity index is calculated on only presence and absence of OTUs, whereas Bray-Curtis index considers OTU abundances. Thus, our results suggest that the elevational differences in tick bacterial diversity are mainly driven by the presence or absence of rarer species, rather than differences in relative OTU abundances. 40% of all amplicons belonged to only four OTUs (of 864 OTUs included in the analysis in total), suggesting a highly skewed abundance distribution. There is currently not a good general understanding of the functional importance of tick bacterial microbiota or how the relative abundance of different OTUs affect its functionality. While we expect that the most abundant OTUs are also the most functionally relevant, relatively rare OTUs may have substantial effects on the hosts or within-host microbial interactions [43].
Currently we can only speculate about the factors that may mediate the observed changes in tick bacterial microbiota diversity along elevational clines. Elevation is strongly associated with temperature, soil moisture, tick host community structure, and land use [1], which might all directly or indirectly shape microbial colonization and thus microbiota diversity. A recent study on I. scapularis ticks in Canada found that ticks at the range expansion front had a different microbiota compared to ticks in the core range [44]. There is evidence that ticks have been undergoing a range expansion to higher elevations in recent years because of climate warming [15], so mountain tops represent a range expansion front [15]. Genetic diversity is typically reduced at range expansion fronts [45], which might contribute to differences in microbiota composition [10]. We directly tested for such effects in our study but found no indication that host genetic diversity or differentiation explains variation in microbial diversity (see Supplementary material S4, Table S4).
The analysis of ecological processes on community turnover suggested that stochastic processes have a strong effect on tick microbiota composition. Homogenizing dispersal, which could be facilitated by cofeeding, (i.e., two or several ticks sharing microbes by feeding in close proximity on the same host), however, played a minor role. Ticks generally only feed from one host individual per life stage, which could provide opportunities for deterministic processes, yet these were rarely observed. Furthermore, no evidence for higher turnover at the expansion edge has been found, although our ability to detect such effects might be limited by the low number of ticks sampled at high elevation sites. Understanding the timeframe which shapes within-host microbiota composition will be essential to better understand the factors that contribute to the elevational differences in microbial diversity observed in our study. However, unfortunately it is next to impossible to longitudinally follow changes in tick microbiota, making this a challenging task.
In addition to elevational and seasonal effects on the bacterial microbiota of ticks we observed substantial differences in tick microbiota composition across tick life stages and sexes. In line with our findings a lower microbiota diversity Elevational Changes in Bacterial Microbiota Structure and Diversity in an Arthropod-Disease… 875 in female ticks has been previously observed in I. scapularis and I. affinis ticks [46,47]. Sex differences in microbiota composition could be due to sex differences in physiology or behavior, including host preference. In many species males have larger home ranges than females, which might lead to exposure to a more diverse bacterial community and explain the higher microbial diversity in male ticks. Yet, to our knowledge movement patterns of male and female ticks have not been studied to date. No difference in seasonal activity is observed between adult male and female ticks [48], but during the nymph stage, female ticks become more engorged, i.e., take up more blood [49], which might partly explain the observed sex differences in microbiota composition [23]. In a previous study, using a joint OTU distribution modeling approach [21], we found that individual-level variation in OTU presence or absence masked the effects of largerscale ecological factors. In contrast, we found here clear signals of variation in microbiota community structure and composition due to elevation and tick life stage and sex. This suggests that the whole community might respond in a different way than individual microbial species in relation to environmental variation, emphasizing the need to consider host-associated microbiota both at the whole community as well as the individual species level.
In conclusion, we found that alpha diversity of tick-associated bacteria decreased with increasing elevation and that variation in within-tick bacterial communities was much more pronounced at lower elevations. Both of these effects were mainly driven by the presence-absence of rarer species rather than differences in relative OTU abundances. Given that bacterial microbiota composition influences the vector competence of ticks [50], understanding the functional consequences of the observed elevational differences in microbiota composition for tick-borne disease dynamics will be an important next step.