# Thermostable Keystone Bacteria Maintain the Functional Diversity of the Ixodes scapularis Microbiome Under Heat Stress

Alejandra Wu-Chuang, Dasiel Obregon, Agustín Estrada-Peña, Alejandro Cabezas-Cruz

► To cite this version:

Alejandra Wu-Chuang, Dasiel Obregon, Agustín Estrada-Peña, Alejandro Cabezas-Cruz. Thermostable Keystone Bacteria Maintain the Functional Diversity of the Ixodes scapularis Microbiome Under Heat Stress. Microbial ecology, 2022, 84 (4), pp.1224-1235. 10.1007/s00248-021-01929-y. hal-03988873

HAL Id: hal-03988873

https://hal.inrae.fr/hal-03988873

Submitted on 14 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Thermostable Keystone Bacteria Maintain the Functional Diversity of the *Ixodes scapularis* Microbiome Under Heat Stress

Alejandra Wu-Chuang<fnref n="1" /> · Dasiel Obregon<fnref n="2" /> · Agustín Estrada-Peña<fnref n="3" /> · Alejandro Cabezas-Cruz<fnref n="1" />

Received: 25 August 2021 / Accepted: 10 November 2021 / Published online: 24 November 2021
© The Author(s), under exclusive licence to Springer Science+Business Media, LLC, part of Springer Nature 2021

## Abstract

Variations in the composition and diversity of tick microbiome due to high temperatures may influence the hierarchy of community members as a response to environmental change. Modifications in the community structure are hypothesized to drive alterations in the presence and/or abundance of functional pathways in the bacterial metagenome. In this study, this hypothesis was tested by using published 16S rRNA datasets of *Ixodes scapularis* males incubated at different temperatures (i.e., 4, 20, 30, and 37 °C) in a laboratory setting. Changes in community structure and functional profiles in response to temperature shifts were measured using co-occurrence networks and metagenome inference. Results from laboratory-reared ticks were then compared with those of field-collected ticks. The results from laboratory-reared ticks showed that high temperature altered the structure of the microbial community and decreased the number of keystone taxa. Notably, four taxa were identified as keystone in all the temperatures, and the functional diversity of the tick microbiome was contained in the four thermostable keystone their associated bacterial taxa. Three of the thermostable keystone taxa were also found in free-living ticks collected in Massachusetts. Moreover, the comparison of functional profiles of laboratory-reared and field-collected ticks revealed the existence of an important set of metabolic pathways that were common among the different datasets. Similar to the laboratory-reared ticks, the keystone taxa identified in field-collected ticks alongside their consortia (co-occurring taxa) were sufficient to retain the majority of the metabolic pathways in the functional profile. These results suggest that keystone taxa are essential in the stability and the functional resiliency of the tick microbiome under heat stress.

**Keywords** *Ixodes scapularis* · Microbiome · Heat stress · Microbial networks · Keystone bacteria · Functional redundancy

## Introduction

The tick microbiome is an important player in tick physiology and vector competence [39, 52]. Efforts have been concentrated on the characterization of the microbial diversity as well as the understanding of the factors that modulate the tick microbiota. Diverse factors such as tick species, gender, developmental stage, blood meal, and host species among others modify the composition of the tick microbiome [12,

37, 45, 49, 51, 52]. Besides these biotic factors, the tick microbiome is also sensible to abiotic factors such as soil type, humidity, season and temperature variation [19, 28, 35, 45, 48, 50, 55]. Among these, the temperature was shown to be an important environmental variable that shapes microbial composition and diversity in *Ixodes scapularis* in a laboratory setting [48]. Specifically, the study of Thapa et al. [48] showed that the incubation of *I. scapularis* in a gradient of temperature, including 4 °C, 20 °C, 30 °C, and 37°C at a constant humidity of 80–95%, have an impact on tick survival, where only 15% of ticks survived at 30 and 37 °C after 10 days of incubation. Moreover, they reported that incubation of ticks at high temperatures (30 and 37 °C) induces large changes in taxonomic composition and decreases bacterial diversity compared to ticks incubated at lower temperatures (4 and 20 °C).

Despite advances in the characterization of the taxonomic composition of the tick microbiome under abiotic stressors such as high temperature, little is known about the impact

<figure><img src="image_2.png" /></figure>

<footnote>Alejandro Cabezas-Cruz
alejandro.cabezas@vet-alfort.fr</footnote>

<footnote n="1">¹ Anses, INRAE, Ecole Nationale Vétérinaire d'Alfort,
UMR BIPAR, Laboratoire de Santé Animale,
Maisons-Alfort F-94700, France</footnote>

<footnote n="2">² School of Environmental Sciences, University of Guelph,
Guelph, Ontario N1G 2W1, Canada</footnote>

<footnote n="3">³ Faculty of Veterinary Medicine, University of Zaragoza,
Zaragoza, Spain</footnote>
of abiotic disturbance on the functional profiles of the tick microbiome. Previous studies showed that the functional profiles of tick microbiomes are highly redundant [15, 16, 40]. For instance, Estrada-Peña et al. [15] showed that up to 198 bacterial genera could contribute to a single pathway in *I. scapularis* microbiome. Furthermore, analysis of the resistance of *I. scapularis* microbiome to disturbing factors such as *Anaplasma phagocytophilum* infection and exposure to an antimicrobial peptide showed that while these factors affected the taxonomic composition, they had limited impact on the functional traits of tick microbiome [16]. These results suggest that a moderate reduction in the abundance of any taxa or random loss of few taxa would have only a marginal impact on the overall function of the tick microbiome since other bacteria can also perform an identical function. Thus, redundancy can contribute to the functional resistance of the microbiome in stressful ecological conditions.

How the functional stability of the tick microbiome is maintained despite taxonomic disturbance is currently unknown. Microbial communities found in the soil, plants, and animals harbor keystone taxa, which play important roles in stabilizing species diversity and ecosystem function [2, 24, 54]. Keystone bacteria were recently identified in tick microbiota using co-occurrence networks [33, 34, 53]. Keystone taxa are importance for the tick microbiome as their depletion using anti-microbiota vaccines was associated to major disturbance in the tick microbial community [33, 34, 53].

We hypothesized that the impact of heat stress on tick microbiota may drive changes in the structure and hierarchical organization of the community as well as and it functionality, as central microbial taxa may lose their keystone-ness (i.e., the quality of being keystone [9, 36] at different temperatures and thus, impact the structure and functional traits of tick microbiome. In this study, 16S rRNA sequence data reported by Thapa et al. [48] was used to test whether keystone bacteria contribute to the functional stability of the microbiome of ticks incubated at different temperatures. Then, the 16S rRNA sequence data reported by Thapa et al. [49] was used to test whether keystone bacteria play a similar role in from wild-caught ticks.

# Methods

## Data Source

16S rRNA gene sequencing datasets published by Thapa et al. [48] and Thapa et al. [49] were used for the analysis. The original studies described the changes in taxonomic composition and bacterial diversity of the microbiome of *Ixodes scapularis* reared in the laboratory or collected from vegetation (detailed below). Datasets were generated

by amplification of hypervariable region four (V4) of the bacterial 16S rRNA gene using barcoded universal primers (515F/806R) and sequenced in an Illumina MiSeq system.

(i) The study of Thapa et al. [48] investigated the impact of environmental temperature on bacterial microbiome by incubation of unfed, female and male *I. scapularis* ticks under different temperatures (4, 20, 30, and 37 °C) at a stable humidity of 80-95% for 10 days or until tick death.

(ii) The study of Thapa et al. [49] examined the influence of sex and location on the bacterial microbiome of free-living *I. scapularis* ticks collected from vegetation in Texas (during autumn) and Massachusetts (during late spring), USA.

As female ticks from the mentioned studies [48, 49] were reported to have a low taxonomic richness, we concentrated our analysis on the microbiome of male ticks. The raw sequences data are publicly available at the sequence read archive (SRA) of the National Center for Biotechnology Information (NCBI) under the project accession numbers PRJNA471905 [48] and PRJNA464062 [49].

## Original Raw Sequences Processing

The 16S data sets were downloaded from SRA repository [29], and deinterlaced using the SRA-Toolkit (http://ncbi.github.io/sra-tools/). The raw sequences (paired-end) were processed and analyzed with the microbiome analysis package QIIME2 (v. 2021.4) [7]. Briefly, fastq files were denoised and merged using DADA2 software package [10]. The resulting amplicon sequence variants (ASVs) were aligned with MAFFT [26]. Taxonomy annotation of ASVs was done using the q2-feature-classifier [6] classify-sklearn naive Bayes taxonomy classifier against the 16S SILVA database (release 132). The resulting taxonomic data tables were collapsed at genus level, then filtered excluding rare taxa (i.e., with less than 10 total reads and presents in less than 30% of samples of each dataset).

## Microbial Diversity and Differential Taxonomic Composition Analysis

Alpha diversity metric (i.e., ASV richness) was estimated using the q2-diversity plugin in Qiime2 environment, and analyzed using Kruskal-Wallis statistical test. For assessment of the differential abundant taxa, the tables with taxonomic profiles were collapsed at genus level, and normalized using centered log ratio (clr) transformation. The comparisons were performed using Kruskal-Wallis statistical test through the ALDEx2 R package [17]. Furthermore, we identified the differential abundant taxa that directly responded to the continuous variable temperature, for which correlational
analyses were performed using aldex.corr function from the ALDEX2 R package. This tool returned the average Pearson coefficient and its significance for each feature.

## Microbial Co-occurrence Networks

Co-occurrence networks were built for each temperature condition using the taxonomic profiles. Co-occurrence networks are graphical visualization of microbial community assemblies where nodes represent bacterial taxa and the connecting edges represent significant correlations between taxa. The Sparse Correlations for Compositional data (SparCC) method [18] implemented in R studio [43] was used to detect significant ($p < 0.05$) correlations (SparCC = 0.5) among bacterial taxa. The visualization and assessment of the networks was performed using the software Gephi 0.9.2 [4]. Network topological parameters were calculated: number of nodes and edges, modularity (the strength of division of a network into modules), network diameter (the shortest path between the two most separated nodes), average degree (the average number of links per node), weighted degree (the sum of the weight of all the edges connected to a node), and clustering coefficient (the degree to which nodes in a network tend to form clusters). For the module (i.e., cluster of taxa that co-occur more often among them than with other bacterial taxa in the network) detection, we used the Louvain method [5] available in Gephi 0.9.2. This algorithm assigns a modularity class value for each node; thus, nodes with the same value of modularity class belong to the same module.

## Identification of Keystone Taxa in Microbial Networks

For the identification of keystone taxa, we used three different criteria: (i) eigenvector centrality, which measures the importance of a node in a co-occurrence network while considering the relevance of their neighbors [44], (ii) ubiquitousness (i.e., bacterial taxa present across all the samples at one condition); and (iii) abundance. Then, cut-off values of 5 and 0.75 were selected for the mean of the abundance and the eigenvector centrality, respectively. Scatter plot was prepared using the software GraphPad 8 Prism (GraphPad Software Inc., San Diego, CA, USA).

## Prediction of Functional Traits in Tick Microbiome

For the metabolic profiling of each sample, PICRUSt2 software (v. 2.4.1) [14] was used for the prediction of functional gene abundances based on 16S rRNA amplicon sequences. Briefly, the ASVs were aligned and placed into a reference tree (NSTI cut-off value of 2) that contains 20000 full 16S rRNA sequences from prokaryotic genomes, which is used to infer gene family copy numbers of each ASVs and finally determine

gene family abundance per sample. Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (KO) [25], Enzyme Classification numbers (EC), and Cluster of Orthologous Genes (COGs) [47] were used as gene family catalogues for the predictions. Pathway profiles were inferred from structured pathway mapping based on MetaCyc database [11].

## Functional Diversity and Contribution of Taxa to the Functional Profile of Tick Microbiome

Functional diversity among the different temperatures was assessed using the open-source software PAST (v. 4.07) [22]. Non-metric multidimensional scaling (NMDs), based on a distance matrix computed with Bray-Curtis similarity index and a PERMANOVA test, was used to compare the functional diversity among the four temperature conditions. A pairwise PERMANOVA test was applied later to determine differences between specific groups. A permutation test was also applied to test differences in the homogeneity of dispersions of the group using the function betadisper from the Vegan R-package.

Linkages between ASVs (collapsed at genus level) and predicted functions (pathways) were assessed using the function "Taxa contribution" from PICRUSt2 metagenome predictions.

The number of shared functions between the whole bacterial community and the bacterial taxa of interest as well as the number of shared keystone taxa amongst the tick microbiome at different temperature conditions was calculated using Venn diagrams through the online tool http://bioinformatics.psb.ugent.be/webtools/Venn/. The number of shared functions between laboratory and field datasets was done in R studio [43] using the package "Venn".

## Definition of Environmental Conditions of Field-Collected Ticks

The environmental variables included average maximum and minimum temperatures and precipitation of the months in which ticks were collected in Texas and Massachusetts, USA. The data were obtained from the TerraClimate repository (http://www.climatologylab.org/terraclimate.html), accessed on July 2021.

## Results

### Decreased Bacterial Diversity at High Temperature is Associated with Shifts in Bacterial Community Assembly

Analysis of bacterial diversity using single nucleotide exact amplicon sequence variants (ASV) showed that the

<figure><img src="image_3.png" /></figure>
microbiome of *I. scapularis* ticks incubated at 37 °C presented a significant decreased microbial richness compared to those incubated at 4 °C, 20 °C, and 30 °C (pairwise Kruskal-Wallis, *p* < 0.05, Supplementary Figure S1), which is consistent with the results of Thapa et al. [48]. Differential analysis showed that the abundance of 44 taxa changed significantly among the four temperatures (Kruskal-Wallis, *p* < 0.05, Supplementary Figure S2). Differences in taxa abundance among the different conditions came predominantly from microbiome of ticks incubated at 37°C, where most taxa presented decreased abundance compared to the microbiome of ticks incubated at 4, 20, and 30 °C. Moreover, correlation analysis showed that most taxa that changed significantly in abundance among the four conditions were strongly correlated to the temperature gradient (Spearman rank correlation, *p* < 0.05; Supplementary Figure S2). Specifically, all these taxa presented a negative correlation with the increase of the environmental temperature.

We then asked whether the changes in diversity and taxonomic composition induced by the temperature shift impacted the bacterial community assembly in *I. scapularis* microbiome. To answer this question, bacterial co-occurrence networks were inferred and used for the visualization and quantification of underlying interactions among bacteria of the tick microbiome. Network analysis showed distinct co-occurrence of bacterial genera across the different temperatures tested. Visual inspection of networks revealed that changes in environmental temperature reshape the bacterial interaction patterns with a decrease in the number of co-occurring taxa as temperature increases (Table 1; Fig. 1). For instance, the number of nodes decreased towards the increase of the temperature, while the number of edges decreased not only at 37 °C but also at 20 °C and 30 °C compared to the network at 4 °C. In addition, the modularity was decreased and the network diameter was increased at 37 °C. Moreover, the number of connections per node was significantly lower at 37 °C. These results suggest that

the decrease in bacterial abundance and diversity induced at 37 °C is also associated with the reduction of microbiome interactions.

## A Fraction of Highly Connected Keystone Taxa is Maintained Along the Temperature Gradient

We next asked whether the microbial networks at different temperatures support keystone taxa. The keystone taxa of each network were identified based on ubiquitousness, abundance, and high eigenvector centrality in the networks. Various ubiquitous bacteria, with high eigenvector centrality (> 0.75) and abundance (clr > 5) were defined as keystone taxa (Fig. 2). However, the results show that the number of keystone taxa decreased with the increase in temperature, since 15, 14, 7, and 4 keystone taxa were identified in ticks incubated at 4, 20, 30, and 37 °C, respectively (Fig. 2). Notably, four taxa (i.e., *Pseudomonas*, *Ralstonia*, *Acinetobacter*, and *Bradyrhizobium*) were identified as keystone in the different temperatures. These taxa are referred hereafter as thermostable keystone taxa. The network analysis further showed that the thermostable keystone taxa were positively connected between them in all the temperatures (Fig. 3).

The 4 thermostable keystone taxa were directly connected to 66, 39, 38, and 15 nodes at 4, 20, 30, and 37 °C, respectively. From these, 41, 32, 31, and 8 taxa were connected to all the thermostable keystone taxa at 4, 20, 30, and 37 °C, respectively (Fig. 4). A detailed analysis of nodes connected to *Pseudomonas*, *Ralstonia*, *Acinetobacter*, and *Bradyrhizobium* in the four conditions (hereafter referred to as “subnetworks”) revealed the existence of a group of 8–10 bacteria that persist as neighbors of each thermostable keystone taxa across the gradient of temperature (Fig. 5). Moreover, the type of connections between the thermostable keystone taxa and their respective persistent neighbors were all positive at 4 °C, except the negative connection between *Pseudomonas* or *Acinetobacter* and *Sediminibacterium*. Interestingly, the positive connection was kept at 20, 30, and 37 °C (Supplementary Table 1). These results suggest that thermostable keystone taxa maintain stable ecological associations with selected members of the microbial communities. It is noteworthy that the keystone taxa and their positively associated neighbor represents only 27.3% (39, total = 143), 24.8% (32, total = 129), 24.4% (31, total = 127), and 27.1% (16, total = 59) of the total of bacteria presents in the co-occurrence networks at 4, 20, 30, and 37 °C, respectively.

## Thermostable Keystone Bacteria Keep the Functional Stability of the Tick Microbiome

To study the impact of different temperatures on the metabolic potential of tick microbiome, we explored the functional diversity from predicted pathway profiles. The

<table><caption>**Table 1** Topological features of the microbial co-occurrence networks from microbiome of male *I. scapularis* incubated at different temperatures</caption><thead><tr><th>Topological features</th><th>4 °C</th><th>20 °C</th><th>30 °C</th><th>37 °C</th></tr></thead><tbody><tr><td>Nodes</td><td>143</td><td>129</td><td>127</td><td>59</td></tr><tr><td>Edges</td><td>1860</td><td>902</td><td>1009</td><td>108</td></tr><tr><td>Positive (%)</td><td>940 (50.5)</td><td>539 (59.7)</td><td>647 (64.1)</td><td>80 (74.1)</td></tr><tr><td>Negative (%)</td><td>920 (49.5)</td><td>363 (40.2)</td><td>362 (35.9)</td><td>28 (25.9)</td></tr><tr><td>Modularity</td><td>16.422</td><td>0.944</td><td>0.624</td><td>0.273</td></tr><tr><td>Network diameter</td><td>5</td><td>6</td><td>5</td><td>8</td></tr><tr><td>Average degree</td><td>26.014</td><td>13.984</td><td>15.89</td><td>3.661</td></tr><tr><td>Weighted degree</td><td>0.486</td><td>2.169</td><td>4.136</td><td>1.36</td></tr><tr><td>Clustering coefficient</td><td>0.611</td><td>0.52</td><td>0.516</td><td>0.609</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Co-occurrence networks in *I. scapularis* under heat stress. Co-occurrence bacterial networks were inferred from the microbiota of ticks incubated at 4, 20, 30, and 37 °C. Nodes represent bacterial taxa and connecting edges stand for a co-occurence correlations (SparCC > 0.5). Node sizes are proportional to the eigenvector centrality value and node colors are based on modularity class metric. Only nodes with at least one connection are displayed. Modules in each network were colored</figcaption></figure>

non-metric multidimensional scaling (NMDS) plot and a
paired PERMANOVA test showed that the functional traits
of ticks incubated at 37 °C were significantly different from
ticks incubated at 4, 20, and 30 °C (F = 5.716, p < 0.05,
Fig. 6). A permutation test for the evaluation of the homo-
geneity of dispersions in the four groups indicated that the
functional traits of ticks incubated at 37 °C were signifi-
cantly more disperse than ticks incubated at 20 and 30 °C
(F = 3.91, p < 0.05).

We next tested the contribution of the thermostable key-
stone taxa and their direct neighbor (Fig. 3) to the functional
profiles of tick microbiome. We found that the thermosta-
ble keystone taxa and their direct neighbors (i.e., including
those with positive and negative interactions) contributed
to more than 99% of the predicted pathways regardless the
incubation temperature (Fig. 7, Supplementary Table 2).
Interestingly, when only the taxa that co-occur positively
with the thermostable keystone taxa in the subnetworks were
considered for the analysis, the contribution to the functional

profile remains identical to the precedent result (Supplemen-
tary Figure S3, Supplementary Table 3). Hence, the thermo-
stable keystone taxa along with the bacteria that co-occur
positively with them were capable to provide most of the
metabolic functions regardless of the temperature gradient.

Field-Collected Ticks Contain Identical Keystone
Taxa Found in Laboratory-Reared Ticks

In order to determine if the keystone taxa found in a labora-
tory setting are the same in field-collected ticks, we used
published 16S rRNA datasets from *I. scapularis* collected
from vegetation in Texas and Massachusetts, USA [49]. The
authors reported that ticks were collected in different sea-
sons to follow the tick activity levels in different locations.
Thus, collection of ticks in Massachusetts was done during
late spring (May 2017) while ticks from Texas were col-
lected during autumn (December 2016, 2017, and November
2017). The average for maximum temperature in Texas was
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 2** Keystone taxa at different temperatures. Scatter plot showing the mean relative abundance, expressed as center log ratio (clr) value, and the eigenvector centrality of each bacterial taxon (dots or triangle) of the co-occurrence networks at different temperatures. The triangles represent bacterial taxa that are present across all the samples of *I. scapularis* incubated at a certain temperature. The vertical dotted line represents the eigencentrality cutoff value of 0.75 and the</figcaption></figure>

of 15 °C, 15 °C, and 22 °C while the minimum temperature was of 7 °C, 6 °C, and 13 °C in December 2016, December 2017, and November 2017, respectively. The monthly rainfall during these months was between 18 and 31 mm. On the other hand, in May 2017 at Massachusetts, the averages for maximum and minimum temperatures were of 16.8 °C and 8.9 °C, respectively while the monthly rainfall during this period was 151 mm. Our analysis revealed that the microbial community of ticks collected in different locations harbour different keystone bacteria (Fig. 8). Surprisingly, these keystone taxa (i.e., *Acinetobacter*, *Ralstonia* and *Bradyrhizobium*) identified in *I. scapularis* of Massachusetts were the same as three of the thermostable keystone taxa found in ticks derived from a laboratory setting. These results suggest that the identification of keystone taxa from laboratory-reared ticks can accurately reflect the microbial community hierarchy of some populations of free-living ticks.

As in the laboratory-reared tick dataset, we analyzed the contribution of the keystone taxa and their direct neighbor to the functional profile of field-caught tick microbiome. We found that the keystone taxa and their subnetwork contribute to 99.3% and 98.1% of predicted pathways in the

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

horizontal dotted line represents the clr mean cutoff value of 5. Red triangles that are above and to the right of the horizontal and vertical dotted line, respectively, represent the identified keystone bacteria at each temperature. The list of keystone taxa is displayed next to the scatter plots. Note that the first four keystone taxa remain the same in all temperatures

microbiome of ticks collected in Massachusetts (Fig. 9A, Supplementary Table 4) and Texas (Fig. 9B, Supplementary Table 4), respectively. Furthermore, we found that the majority of the metabolic pathways (i.e., 373, total 430; 86.7%) were common among laboratory-reared and field-collected ticks (Fig. 9C, Supplementary Table 5).

## Discussion

Several studies have demonstrated that the taxonomic composition and diversity of tick microbiome is highly influenced by diverse biotic [12, 37, 52] and abiotic factors [19, 45, 48, 50]. Yet, the functional traits of tick microbiome and their possible changes associated to environmental factors remains largely unknown. In this study, the impact of heat stress on the structure and functional profiles of the microbial communities of laboratory-reared and field-collected *I. scapularis* was assessed. To this aim, we used published and publicly accessible 16S rRNA datasets [48, 49]. As the study of Thapa et al. [48] reported that the microbiome of laboratory-reared female ticks was dominated by *Rickettsia* and
<figure><img src="image_9.png" /><figcaption>**Fig. 4** Number of neighbor taxa shared by the thermostable keystone taxa in the different temperatures. Venn diagram showing the number of bacteria that are common or unique among the taxa connected to the four thermostable keystone taxa *Pseudomonas*, *Acinetobacter*, *Bradyrhizobium*, and *Ralstonia* in the different temperatures</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 5** Number of taxa that persist as neighbors for each thermostable keystone taxa across the gradient of temperature. The number of bacteria genera shared by each thermostable keystone taxa, *Pseudomonas*, *Ralstonia*, *Acinetobacter*, and *Bradyrhizobium*, in the four temperatures is displayed as a Venn diagram</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 6** Impact of temperature on the diversity of the functional profiles. Non-metric Multidimensional Scaling (NMDS) based on Bray-Curtis similarity matrix of predicted pathways for each sample at 4 (dots), 20 (triangle), 30 (square), and 37 °C (diamond). Groups are surrounded by ellipses (95% Confidence interval around centroid). Permutational analysis of variance (PERMANOVA) *F* and *p* values are displayed at the top right of NMDS plot</figcaption></figure>

hence, it had a lower microbial diversity compared to male ticks, we decided to focus only on male dataset. The higher bacterial diversity of the *I. scapularis* male can spur more information of changes in microbial interactions and their contribution to the functional profile. Our results support

the stability of the functional microbiome despite changes in taxonomic composition due to heat stress. Notably, a fraction of microbes composed by four thermostable keystone taxa and their invariably associated microbial partners carry almost the totality of functions of the metagenome underlying the functional redundancy of tick microbiome.

Consistent with the results of Thapa et al. [48], we found a decrease in the diversity and relative abundance of bacterial taxa in the tick microbiome under heat stress. Community modifications were also reflected in the co-occurrence networks of ticks incubated at the different temperatures (i.e., 4, 20, 30, and 37 °C). The greater loss of bacterial taxa and decreased interconnection between them was observed in the networks inferred from microbial communities of ticks incubated at 37 °C. This microbial network presented the highest diameter and the lowest modularity, compared to microbial networks of ticks incubated at 4, 20, and 30 °C. Estrada-Peña et al. [16] showed that disturbing factors such as anti-tick immunity, *Anaplasma phagocytophilum* infection and antimicrobial peptide treatment increased the diameter of tick microbiome networks, which was associated with lower network robustness compared to the networks without disturbance. Furthermore, high network modularity has been proposed to be a proxy of microbiome stability [27]. Network parameters thus suggest that heat stress reduces the robustness and stability of tick microbial communities. Shifts on microbial interaction patterns as result of temperature variation has been previously described in soil bacterial communities [3, 8, 13, 20, 54].

Co-occurrence networks can also be used to identify keystone taxa which have high influence and importance in the microbial community [1, 2, 23]. We used three parameters (i.e., eigenvector centrality, ubiquitousness, and abundance) to infer the keystoneness of bacterial taxa in the tick microbiome [33]. As discussed by Wu-Chuang et al. [53], these criteria differ from those used by other authors to infer bacterial keystoneness in microbial communities. The impact of keystone taxa, predicted with the criteria used in this study, on tick physiology [33] and tick microbiome structure and function [33, 34] has been validated using anti-tick microbiota vaccines [33, 52, 53]. In this study, the number of predicted keystone taxa in tick microbiome decreased with high temperatures suggesting that environmental temperature variation can change the influence the importance of some bacteria in the community.

Despite a decrease in the number of keystone taxa, we found a core of four keystone bacterial genera that was systematically identified across the gradient of temperature. These thermostable keystone taxa were associated to a set of bacteria that together contributed to more than 99% of the predicted pathways regardless the incubation temperature. These findings suggest that thermostable keystone taxa may provide functional stability to the tick microbiome under heat stress. Similar result were found in soil bacterial
<figure><img src="image_14.png" /><figcaption>**Fig. 7** Functional contribution of thermostable keystone taxa and their direct neighbor. Venn diagram showing the number of metabolic pathways contributed by the thermostable keystone taxa and their positively and negatively associated direct neighbors (subnetwork) and compared to the total of predicted functions for the whole bacterial community at 4, 20, 30, and 37 °C</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 8** Keystone taxa in field-collected ticks. Scatter plot showing the mean relative abundance (expressed as clr) and the eigencentrality of each bacterial taxon (dots or triangle) of the microbial community from *I. scapularis* collected from Massachusetts (A) and Texas (B). The red dots or triangle represent bacterial taxa that are present across all the samples of *I. scapularis* collected from different locali-</figcaption></figure>

tions. The vertical dotted line represents the eigencentrality cutoff value of 0.75 and the horizontal dotted line represents the clr mean cutoff value of 5. Red triangles that are above and at the right of the horizontal and vertical dotted line, respectively, represent the identified keystone bacteria at each location. The list of keystone taxa is displayed next to the scatter plots

communities of two sites of the Atacama Desert, west of the Andes Mountains [32]. These authors showed that beside a persistent group of bacteria, the most abundant functional roles were conserved despite changes in environmental variables [32]. It is noteworthy that the keystone taxa with their positively correlated neighbors represent only a small fraction of the whole bacterial community suggesting that several bacteria contribute to the same set of functions in the tick microbiome, as previously described by Obregón et al. [40] and [15, 16]. A limitation here is that the accuracy

of functional inferences based in PICRUSt2 varies across sample types and functional categories [46]. Although the accuracy of PICRUSt2 predictions for tick microbiome has not been assessed, a recent study reported a sharp decrease in PICRUSt2 performance for inference in non-human (i.e., gorilla, mouse, and chicken) and environmental samples (i.e., soil), compared with human samples [46]. Having said that, it is worth mentioning that functional inference of tick microbiome using PICRUSt2 have been validated by PCR and qPCR in two independent studies [33, 34]. In both

<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 9** Field-collected and laboratory-reared ticks share a functional core. Venn diagram showing the number of metabolic pathways contributed by the keystone taxa and their direct neighbors compared to the total number of predicted pathways for the whole bacterial community of ticks collected from vegetation of Texas (A) and Massachusetts (B). Venn diagram showing the number of functional pathways shared among tick microbiomes of laboratory-reared and field collected ticks (C)</figcaption></figure>

studies, validation by PCR/qPCR concurred with PICRUSt2 predictions of pathway presence and abundance [33, 34].

Notably, three of the thermostable keystone taxa identified in laboratory-reared ticks were also found in field-collected *I. scapularis* from Massachusetts. These results suggest that the microbiome of ticks studied in a controlled environment can reflect accurately the hierarchical organization of microbial communities in some populations of free-living ticks. Similarly, Narasimhan et al. [38] found that field-collected and lab-reared nymphs presented similar gut microbial composition at phylum and genus levels. Although not in the microbiome context, other studies have also shown that laboratory data can reflect physiological or behavioral features of field-collected ticks [41, 42]. In contrast, *I. scapularis* ticks collected from Texas were found to harbour *Rickettsia* as the only keystone taxon. Different abiotic or biotic factors influencing tick populations and their microbial communities in Texas (collected in late spring) and Massachusetts (collected in autumn) may modify microbiota composition and/or influence the hierarchical organization of microbial communities.

## Conclusions

In this study, we found that heat stress induces changes in the tick microbial community including a reduction in the number of co-occurring bacteria and keystone taxa. Despite drastic alterations in the bacterial community structure, a core of keystone taxa persisted across the gradient of temperature. These thermostable keystone taxa and their shared set of positively correlated neighbors represented approximately 99%

of the functional diversity of the tick microbiome. The results obtained under laboratory conditions were partially validated in the field, as one of the tested populations of field-collected ticks harbored identical keystone taxa as the laboratory-reared ticks. The results suggest that keystone taxa keep the functional stability of the tick microbiome under environmental variations in temperature. Furthermore, we found the existence of a functional core shared by laboratory-reared and field-collected ticks, underlying that redundancy can contribute to the functional resistance of the microbiome in stressful ecological conditions. In summary, our study provides an insight into the role of keystone taxa in keeping the functional stability of the tick microbiome under abiotic stressors.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00248-021-01929-y.

**Acknowledgements** UMR BIPAR is supported by the French Government's Investissement d'Avenir program, Laboratoire d'Excellence "Integrative Biology of Emerging Infectious Diseases" (grant no. ANR-10-LABX-62-IBEID). Alejandra Wu-Chuang is supported by Programa Nacional de Becas de Postgrado en el Exterior "Don Carlos Antonio López" (grant no. 205/2018).

**Authors' Contributions** A.C.C. and D.O. conceived the study. A.W.C., D.O., and A.E.P. performed the analyses. A.W.C. prepared the figures and drafted the manuscript. All authors revised and accepted the last version of the manuscript. All authors read and approved the final manuscript.

**Funding** This research received no specific funding.

**Data Availability** All material relevant to this publication are available in the manuscript.

**Code Availability** Not applicable
## Declarations

**Ethics Approval** Not applicable

**Consent to Participate** Not applicable

**Consent for Publication** Not applicable

**Conflicts of Interest** The authors declare that they have no conflicts of interests.

## References

1. Agler, MT, Ruhe J, Kroll S, Morhenn C, Kim ST, Weigel D, Kemen EM (2016) Microbial hub taxa link host and abiotic factors to plant microbiome variation. PLoS Biology 14(1). https://doi.org/10.1371/journal.pbio.1002352

2. Banerjee S, Schlaeppi K, van der Heijden MGA (2018) Keystone taxa as drivers of microbiome structure and functioning. Nat Rev Microbiol 16(9):567–576. https://doi.org/10.1038/s41579-018-0024-1

3. Barreiro A, Lombao A, Martín A, Cancelo-González J, Carbalas T, Díaz-Raviña M (2020) Soil heating at high temperatures and different water content: Effects on the soil microorganisms. Geosci (Switzerland) 10(9):1–17. https://doi.org/10.3390/geosciences10090355

4. Bastian M, Heymann S, Jacomy M (2009) Gephi: an open source software for exploring and manipulating networks. In Third international AAAI conference on weblogs and social media

5. Blondel V, Guillaume J, Lambiotte R, Mech E (2008) Fast unfolding of communities in large networks. J Stat Mech: Theory Exp 2008: P10008. http://findcommunities.googlepages.com

6. Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, Huttley GA, Gregory Caporaso J (2018) Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2's q2-feature-classifier plugin. Microbiome 6(1):1–17. https://doi.org/10.1186/s40168-018-0470-z

7. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, Alexander H, Alm EJ, Arumugam M, Asnicar F, Bai Y, Bisanz JE, Bittinger K, Brejnrod A, Brislawn CJ, Brown CT, Callahan BJ, Caraballo-Rodríguez AM, Chase J, Cope EK, Da Silva R, Diener C, Dorrestein PC, Douglas GM, Durall DM, Duvallet C, Edwardson CF, Ernst M, Estaki M, Fouquier J, Gauglitz JM, Gibbons SM, Gibson DL, Gonzalez A, Gorlick K, Guo J, Hillmann B, Holmes S, Holste H, Huttenhower C, Huttley GA, Janssen S, Jarmusch AK, Jiang L, Kaehler BD, Kang KB, Keefe CR, Keim P, Kelley ST, Knights D, Koester I, Kosciolek T, Kreps J, Langille MGI, Lee J, Ley R, Liu YX, Loftfield E, Lozupone C, Maher M, Marotz C, Martin BD, McDonald D, McIver LJ, Melnik AV, Metcalf JL, Morgan SC, Morton JT, Naimey AT, Navas-Molina JA, Nothias LF, Orchanian SB, Pearson T, Peoples SL, Petras D, Preuss ML, Pruesse E, Rasmussen LB, Rivers A, Robeson MS, Rosenthal P, Segata N, Shaffer M, Shiffer A, Sinha R, Song SJ, Spear JR, Swafford AD, Thompson LR, Torres PJ, Trinh P, Tripathi A, Turnbaugh PJ, Ul-Hasan S, van der Hooft JJJ, Vargas F, Vázquez-Baeza Y, Vogtmann E, von Hippel M, Walters W, Wan Y, Wang M, Warren J, Weber KC, Williamson CHD, Willis AD, Xu ZZ, Zaneveld JR, Zhang Y, Zhu Q, Knight R, Caporaso JG (2019) Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol 37(8):852–857. https://doi.org/10.1038/s41587-019-0209-9

8. Burman E, Bengtsson-Palme J (2021) Microbial community interactions are sensitive to small changes in temperature. Front Microbiol 12:672910. https://doi.org/10.3389/fmicb.2021.672910

9. Cagua EF, Wootton KL, Stouffer DB (2019) Keysteness, centrality, and the structural controllability of ecological networks. J Ecol 107(4):1779–1790. https://doi.org/10.1111/1365-2745.13147

10. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP (2016) DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods 13(7):581–583. https://doi.org/10.1038/nmeth.3869

11. Caspi R, Billington R, Fulcher CA, Keseler IM, Kothari A, Krummenacker M, Latendresse M, Midford PE, Ong Q, Ong WK, Paley S, Subhraveti P, Karp PD (2018) The MetaCyc database of metabolic pathways and enzymes. Nucleic Acids Res 46(D1):D633–D639. https://doi.org/10.1093/nar/gkx935

12. Chicana B, Couper LI, Kwan JY, Tahiraj E, Swei A (2019) Comparative microbiome profiles of sympatric tick species from the far-western United States. Insects 10(10):1–12. https://doi.org/10.3390/insects10100353

13. Donhauser J, Niklaus PA, Rousk J, Larose C, Frey B (2020) Temperatures beyond the community optimum promote the dominance of heat-adapted, fast growing and stress resistant bacteria in alpine soils. Soil Biol Biochem 148:107873. https://doi.org/10.1016/j.soilbio.2020.107873

14. Douglas GM, Maffei VJ, Zaneveld JR, Yurgel SN, Brown JR, Taylor CM, Huttenhower C, Langille MGI (2020) PICRUSt2 for prediction of metagenome functions. Nat Biotechnol 38(6). https://doi.org/10.1038/s41587-020-0548-6

15. Estrada-Peña A, Cabezas-Cruz A, Obregón D (2020) Behind taxonomic variability: the functional redundancy in the tick microbiome. Microorganisms 8(11):1–16. https://doi.org/10.3390/microorganisms8111829

16. Estrada-Peña A, Cabezas-Cruz A, Obregón D (2020) Resistance of tick gut microbiome to anti-tick vaccines, pathogen infection and antimicrobial peptides. Pathogens 9(4):1–17. https://doi.org/10.3390/pathogens9040309

17. Fernandes DA, Reid J, Macklaim MJ, McMurrough TA, Edgell DR, Gloor BG (2014) Unifying the analysis of high-throughput sequencing datasets: characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. Microbiome 2:15. https://doi.org/10.1186/2049-2618-2-15

18. Friedman J, Alm EJ (2012) Inferring correlation networks from genomic survey data. PLoS Comput Biol 8(9). https://doi.org/10.1371/journal.pcbi.1002687

19. Gall CA, Scoles GA, Magori K, Mason KL, Brayton KA (2017) Laboratory colonization stabilizes the naturally dynamic microbiome composition of field collected *Dermacentor andersoni* ticks. Microbiome 5(1):133. https://doi.org/10.1186/s40168-017-0352-9

20. García FC, Bestion E, Warfield R, Yvon-Durochera G (2018) Changes in temperature alter the relationship between biodiversity and ecosystem functioning. Proc Natl Acad Sci USA 115(43):10989–10994. https://doi.org/10.1073/pnas.1805518115

21. Greay TL, Gofton AW, Paparini A, Ryan UM, Oskam CL, Irwin PJ (2018) Recent insights into the tick microbiome gained through next-generation sequencing. Parasites and Vectors 11(1). https://doi.org/10.1186/s13071-017-2550-5

22. Hammer Ø, Harper DAT, Ryan PD (2001) PAST: Paleontological Statistics Software Package for Education and Data Analysis. Palaeontol Electron 4(1): 9. http://palaeo-electronica.org/2001_1/past/issue1_01.htm

23. Herren CM, McMahon KD (2018) Keystone taxa predict compositional change in microbial communities. Environ Microbiol 20(6):2207–2217. https://doi.org/10.1111/1462-2920.14257
24. Jones P, Garcia BJ, Furches A, Tuskan GA, Jacobson D (2019) Plant host-associated mechanisms for microbial selection. In Frontiers in Plant Science, 10, 862. Frontiers Media S.A. https://doi.org/10.3389/fpls.2019.00862

25. Kanehisa M, Goto S (2000) KEGG: Kyoto Encyclopedia of Genes and Genomes. In Nucleic Acids Research, 28(1), pp. 27–30. Oxford University Press. https://doi.org/10.1093/nar/28.1.27

26. Katoh K, Misawa K, Kuma KI, Miyata T (2002) MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res 30(14):3059–3066. https://doi.org/10.1093/nar/gkf436

27. Krause AE, Frank KA, Mason DM, Ulanowicz RE, Taylor WW (2003) Compartments revealed in food-web structure. Nature 426:282–285. https://doi.org/10.1038/nature02115

28. Lalzar I, Harrus S, Mumcuoglu KY, Gottlieb Y (2012) Composition and seasonal variation of *Rhipicephalus turanicus* and *Rhipicephalus sanguineus* bacterial communities. Appl Environ Microbiol 78(12):4110–4116. https://doi.org/10.1128/AEM.00323-12

29. Leinonen R, Sugawara H, Shumway M (2011) The sequence read archive. Nucleic Acids Res 39(SUPPL. 1). https://doi.org/10.1093/nar/gkq1019

30. Liu Z, Wei H, Zhang J, Saleem M, He Y, Zhong J, Ma R (2021) Higher sensitivity of soil microbial network than community structure under acid rain. Microorganisms 9:118. https://doi.org/10.3390/microorganisms9010118

31. Ma B, Wang Y, Ye S, Liu S, Stirling E, Gilbert JA, Faust K, Knight R, Jansson JK, Cardona C, Röttjers L, Xu J (2020) Earth microbial co-occurrence network reveals interconnection pattern across microbiomes. Microbiome 8:82. https://doi.org/10.1186/s40168-020-00857-2

32. Mandakovic D, Rojas C, Maldonado J, Latorre M, Travisany D, Delage E, Bihouée A, Jean G, Díaz FP, Fernández-Gómez B, Cabrera P, Gaete A, Latorre C, Gutiérrez RA, Maass A, Cambiazo V, Navarrete SA, Eveillard D, González M (2018) Structure and co-occurrence patterns in microbial communities under acute environmental stress reveal ecological factors fostering resilience. Sci Rep 8(1):1–12. https://doi.org/10.1038/s41598-018-23931-0

33. Mateos-Hernández L, Obreg D, Maye J, Borneres J, Versille N, de la Fuente JL, Estrada-Peña A, Hodžić A, Šimo L, Cabezas-Cruz A (2020) Anti-tick microbiota vaccine impacts *Ixodes ricinus* performance during feeding. Vaccines 8(4):1–21. https://doi.org/10.3390/vaccines8040702

34. Mateos-Hernández L, Obregon D, Wu-Chuang A, Bornères J, Versillé N, de la Fuente J, Diaz-Sanchez S, Bermúdez-Humarán LG, Torres-Maravilla E, Estrada-Peña A, Hodžić A, Šimo L, Cabezas-Cruz A (2021) Anti-microbiota vaccines modulate the tick microbiome in a taxon-specific manner. Front Immunol. https://doi.org/10.3389/fimmu.2021.704621

35. Menchaca AC, Visi DK, Strey OF, Teel PD, Kalinowski K, Allen MS, Williamson PC (2013) Preliminary assessment of microbiome changes following blood-feeding and survivorship in the *Amblyomma americanum* nymph-to-adult transition using semiconductor sequencing. PLoS ONE 8(6):e67129. https://doi.org/10.1371/journal.pone.0067129

36. Mills LS, Soulé ME, Doak DF (1993) The keystone-species concept in ecology and conservation. Bioscience 43(4):219–224. https://doi.org/10.2307/1312122

37. Moore TC, Pulscher LA, Caddell L, von Fricken ME, Anderson BD, Gonchigoo B, Gray GC (2018) Evidence for transovarial transmission of tick-borne rickettsiae circulating in Northern Mongolia. PLoS Negl Trop Dis 12(8):e0006696. https://doi.org/10.1371/journal.pntd.0006696

38. Narasimhan S, Rajeevan N, Liu L, Zhao YO, Heisig J, Pan J, Eppler-Epstein R, Deponte K, Fish D, Fikrig E (2014) Gut microbiota of the tick vector *Ixodes scapularis* modulate colonization of the Lyme

disease spirochete. Cell Host Microbe 15(1):58–71. https://doi.org/10.1016/j.chom.2013.12.001

39. Narasimhan S, Fikrig E (2015) Tick microbiome: the force within. Trends Parasitol 31(7):315–323. https://doi.org/10.1016/j.pt.2015.03.010

40. Obregón D, Bard E, Abrial D, Estrada-Peña A, Cabezas-Cruz A (2019) Sex-Specific Linkages Between Taxonomic and Functional Profiles of Tick Gut Microbiomes. Front Cell Infect Microbiol 9(August):1–16. https://doi.org/10.3389/fcimb.2019.00298

41. Ogden NH, Lindsay LR, Beauchamp G, Charron D, Maarouf A, O’Callaghan CJ, Waltner-Toews D, Barker IK (2004) Investigation of relationships between temperature and developmental rates of tick *Ixodes scapularis* (Acari: Ixodidae) in the laboratory and field. J Med Entomol 41(4):622–633. https://doi.org/10.1603/0022-2585-41.4.622

42. Peavey CA, Lane RS (1996) Field and laboratory studies on the timing of oviposition and hatching of the western black-legged tick, *Ixodes pacificus* (Acari: Ixodidae). In Experimental & Applied Acarology 20:695–711

43. RStudio Team (2020) RStudio: Integrated Development for R. RStudio, PBC, Boston, MA. http://www.rstudio.com/

44. Ruhnau B (2000) Eigenvector-centrality—a node-centrality? Soc Netw 22:357–365. https://doi.org/10.1016/S0378-8733(00)00031-9

45. Swei A, Kwan JY (2017) Tick microbiome and pathogen acquisition altered by host blood meal. ISME J 11(3):813–816. https://doi.org/10.1038/ismej.2016.152

46. Sun S, Jones RB, Fodor AA (2019) Inference based PICRUSt accuracy varies across sample types and functional categories. bioRxiv. https://doi.org/10.1101/655746

47. Tatusov RL, Galperin MY, Natale DA, Koonin EV (2000) The COG database: A tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res 28(1):33–36. https://doi.org/10.1093/nar/28.1.33

48. Thapa S, Zhang Y, Allen MS (2018) Effects of temperature on bacterial microbiome composition in *Ixodes scapularis* ticks. MicrobiologyOpen 8(5). https://doi.org/10.1002/mbo.719

49. Thapa S, Zhang Y, Allen MS (2019) Bacterial microbiomes of *Ixodes scapularis* ticks collected from Massachusetts and Texas, USA. BMC Microbiol 19(1). https://doi.org/10.1186/s12866-019-1514-7

50. Trout Fryxell RT, DeBruyn JM (2016) The microbiome of *Ehrlichia*-infected and uninfected lone star ticks (*Amblyomma americanum*). PLoS ONE 11(1):1–19. https://doi.org/10.1371/journal.pone.0146651

51. van Treuren W, Ponnusamy L, Brinkerhoff RJ, Gonzalez A, Parobek CM, Juliano JJ, Andreadis TG, Falco RC, Ziegler LB, Hathaway N, Keeler C, Emch M, Bailey JA, Roe RM, Apperson CS, Knight R, Meshnick SR (2015) Variation in the microbiota of *Ixodes* ticks with regard to geography, species, and sex. Appl Environ Microbiol 81(18):6200–6209. https://doi.org/10.1128/AEM.01562-15

52. Wu-Chuang A, Hodžić A, Mateos-Hernández L, Estrada-Peña A, Obregón D, Cabezas-Cruz A (2021a) Current debates and advances in tick microbiome research. Curr Res Parasitol Vector-Borne Dis 100036. https://doi.org/10.1016/j.crpvbd.2021.100036

53. Wu-Chuang A, Obregon D, Mateos-Hernández L, Cabezas-Cruz A (2021) Anti-tick microbiota vaccines: how can this actually work? Biologia. https://doi.org/10.1007/s11756-021-00818-6

54. Xun W, Liu Y, Li W, Ren Y, Xiong W, Xu Z, Zhang N, Miao Y, Shen Q, Zhang R (2021) Specialized metabolic functions of keystone taxa sustain soil microbiome stability. Microbiome 9(1):1–15. https://doi.org/10.1186/s40168-020-00985-9

55. Zolnik CP, Prill RJ, Falco RC, Daniels TJ, Kolokotronis SO (2016) Microbiome changes through ontogeny of a tick pathogen vector. Mol Ecol 25(19):4963–4977. https://doi.org/10.1111/mec.13832