# Temporal patterns in Ixodes ricinus microbial communities: an insight into tick-borne microbe interactions

E. Lejal, J. Chiquet, J. Aubert, S. Robin, A. Estrada-Peña, O. Rue, C.

Midoux, M. Mariadassou, Xavier Bailly, A. Cougoul, et al.

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

► To cite this version:

E. Lejal, J. Chiquet, J. Aubert, S. Robin, A. Estrada-Peña, et al. Temporal patterns in Ixodes ricinus microbial communities: an insight into tick-borne microbe interactions. Microbiome, 2021, 9 (1), pp.153. 10.1186/s40168-021-01051-8. hal-03287263

HAL Id: hal-03287263

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

Submitted on 19 Aug 2021

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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Temporal patterns in *Ixodes ricinus*
microbial communities: an insight into tick-
borne microbe interactions

E. Lejal¹, J. Chiquet², J. Aubert², S. Robin², A. Estrada-Peña³, O. Rue⁴,⁵, C. Midoux⁴,⁵,⁶, M. Mariadassou⁴,⁵, X. Bailly⁷,⁸, A. Cougoul⁷,⁸, P. Gasqui⁷,⁸, J. F. Cosson¹, K. Chalvet-Monfray⁷,⁸, M. Vayssier-Taussat⁹ and T. Pollet¹⁰*

## Abstract

**Background:** Ticks transmit pathogens of medical and veterinary importance and are an increasing threat to human and animal health. Assessing disease risk and developing new control strategies requires identifying members of the tick-borne microbiota as well as their temporal dynamics and interactions.

**Methods:** Using high-throughput sequencing, we studied the *Ixodes ricinus* microbiota and its temporal dynamics. 371 nymphs were monthly collected during three consecutive years in a peri-urban forest. After a Poisson lognormal model was adjusted to our data set, a principal component analysis, sparse network reconstruction, and differential analysis allowed us to assess seasonal and monthly variability of *I. ricinus* microbiota and interactions within this community.

**Results:** Around 75% of the detected sequences belonged to five genera known to be maternally inherited bacteria in arthropods and to potentially circulate in ticks: *Candidatus Midichloria, Rickettsia, Spiroplasma, Arsenophonus* and *Wolbachia*. The structure of the *I. ricinus* microbiota varied over time with interannual recurrence and seemed to be mainly driven by OTUs commonly found in the environment. Total network analysis revealed a majority of positive partial correlations. We identified strong relationships between OTUs belonging to *Wolbachia* and *Arsenophonus*, evidence for the presence of the parasitoid wasp *Ixodiphagus hookeri* in ticks. Other associations were observed between the tick symbiont *Candidatus Midichloria* and pathogens belonging to *Rickettsia*. Finally, more specific network analyses were performed on TBP-infected samples and suggested that the presence of pathogens belonging to the genera *Borrelia, Anaplasma* and *Rickettsia* may disrupt microbial interactions in *I. ricinus*.

**Conclusions:** We identified the *I. ricinus* microbiota and documented marked shifts in tick microbiota dynamics over time. Statistically, we showed strong relationships between the presence of specific pathogens and the structure of the *I. ricinus* microbiota. We detected close links between some tick symbionts and the potential presence of either pathogenic *Rickettsia* or a parasitoid in ticks. These new findings pave the way for the development of new strategies for the control of ticks and tick-borne diseases.

## Introduction

Ticks are vectors of many zoonotic pathogens and are an important and increasing threat to human and animal

health. While these arthropods are among the main vectors of pathogens that affect humans and animals worldwide, it is now well established that tick-borne pathogens (TBPs) coexist with many other microorganisms in ticks. These other members of the microbiota, commensal and/or symbionts, are likely to confer multiple detrimental, neutral, or beneficial effects to their

<footnote>* Correspondence: thomas.pollet@inrae.fr
¹⁰UMR ASTRE, CIRAD, INRAE, Campus de Baillarguet, Montpellier, France
Full list of author information is available at the end of the article</footnote>
tick hosts and can play various roles in nutritional adaptation, development, reproduction, defense against environmental stress and immunity [1, 2]. They may also interact with tick-borne pathogens thereby influencing the competence of the tick vector [3, 4]. Identifying and characterising the tick microbiota is thus crucial to better understand tick-microbe interactions. With the development of high-throughput sequencing technologies, the number of studies dealing with tick microbiota has considerably increased in the past decade and revealed unexpected microbial diversity in ticks [5–14]. The microbiota of several tick species belonging to the genera *Ixodes*, *Dermacentor*, *Haemaphysalis*, *Rhipicephalus* and *Amblyomma* has been studied [15], revealing key details on microbial communities in ticks and improving our knowledge of tick microbiota ecology in general. However, while *Ixodes ricinus* is the main tick species present in western Europe, able to transmit the widest range of pathogens (including the agent of Lyme disease, *Borrelia burgdorferi* s.l.), the ecological factors driving variations in its microbiota have been the subject of few studies [16]. Moreover, several of the previous studies of the *Ixodes ricinus* microbiota [6, 12, 16–22], may have overestimated microbiota diversity due to contamination during the extraction step and the absence of negative controls in their analysis [22]. Studies of *I. ricinus* microbiota pointed to complex microbial assemblages inhabiting ticks including pathogens, specific endosymbionts, commensal and environmental microorganisms. As already reported for pathogens [23–26], the other members of these biological assemblages are probably dynamic and likely to vary over time. While the observation scale affects our view of the dynamics of tick microbiota [27], little information is currently available on the temporal patterns of the *I. ricinus* microbiota. Does the *I. ricinus* microbiota vary with the season? Are these potential temporal patterns repeated annually? Answering these questions is a crucial first step to better understand both the tick and tick microbiota ecology. While it is now accepted that tick microbiota may also play a role in driving transmission or multiplication of tick-borne pathogens [28], little information is currently available on the interactions between *I. ricinus*-borne microbiota members and on the potential co-occurrence of pathogens and the *I. ricinus* microbiota [16]. This information is needed to identify potential strategies for the control of ticks and tick-borne diseases in the future. *I. ricinus* nymphs were collected monthly in three consecutive years in a peri-urban forest. High-throughput sequencing was used to identify the *I. ricinus* microbiota and its temporal dynamics. Using multivariate and partial correlation network analyses, the aim of this study was to identify direct statistical associations between members of the microbiota, including pathogenic

genera, and to assess the influence of the presence of TBPs on tick microbiota structure and interactions.

# Material and methods

## Tick collection

Questing *Ixodes ricinus* nymphs were collected for 3 years by dragging (from April 2014 to May 2017) in the Sénart forest in the south of Paris. More details on the sampling location and design, and tick collection, are available in Lejal et al. [26].

## Tick homogenisation and DNA extraction

In total, 998 nymphs have been collected over the 3 years [26]. As detailed in our previous studies [26, 29], ticks were first washed once in ethanol 70% for 5 min and rinsed twice in sterile MilliQ water. They were then individually homogenised in 375 µL of Dulbecco's modified Eagle's medium with decomplemented foetal calf serum (10%) and six steel beads using the homogeniser Precellys®24 Dual (Bertin, France) at 5500 rpm for 20 s. DNA extraction was performed on 100 µL of tick homogenate, using the NucleoSpin® Tissue DNA extraction kit (Macherey-Nagel, Germany).

## DNA amplification and multiplexing

Among the 998 nymphs, 557 have been chosen to characterise the temporal dynamics of *I. ricinus* microbiota. Thanks to our previous study using microfluidic PCR [26], we knew the infection rate of all collected ticks. For each collecting month, while all infected nymphs have been compulsory integrated in the analysis, non-infected nymphs (at least 15 per month) were added until reaching a minimum of 30 ticks (infected and non-infected) per month, when it was possible. When less than 30 ticks were collected in a month, all ticks were added in the analysis. In addition to tick samples, 45 negative controls have been performed to distinguish tick microbial OTUs from contaminants [22]. As detailed in Lejal et al. [22], DNA amplifications were performed on the V4 region of the 16S rRNA gene using the primer pair used by Galan et al. [30] (16S-V4F: 5'-GTGCCAGCMG CCGCGGTAA-3' and 16S-V4R: 5'-GGACTACHVGGG TWTCTAATCC-3'), producing a 251-bp amplicon. Different 8 bp indexes were added to primers allowing *in fine* the amplification and multiplexing of all samples. All the PCR amplifications were carried out using the Phusion® High-Fidelity DNA Polymerase amplification kit (Thermo Scientific, Lithuania). For each sample, 5 µL of DNA extract were amplified in a 50 µL final reaction volume, containing 1X Phusion HF buffer, 0.2 µM of dNTPs, 0.2 U/mL of Phusion DNA polymerase and 0.35 µM of forward and reverse primer. The following thermal cycling procedure was used: initial denaturation at 98 °C for 30 s, 35 cycles of denaturation at 98 °C for
10 s, annealing at 55 °C for 30 s, followed by extension at 72 °C for 30 s. The final extension was carried out at 72 °C for 10 min. PCR products were checked on 1.5% agarose gels, cleaned and quantified and finally pooled at equimolar concentrations and sent to the sequencing platform (GenoScreen, France) [22].

## Sequencing and data processing

The equimolar mix was concentrated and sequenced by GenoScreen (Lille, France) using MiSeq Illumina 2 × 250 bp chemistry. All the quality controls and different steps of sequence analyses have been performed (more details are available in [22]). Based on results obtained in the 45 negative controls, sequences considered as contaminants were removed from the dataset [22] (Additional file 1). Due to this large OTU filtration, we identified 186 samples with less than 500 sequences. We considered this number of sequences too low to be analysed and we thus removed these samples. Finally, the microbiota of 371 *I. ricinus* nymphs was analysed from a final dataset composed of 907,941 sequences.

## Data preprocessing

Three distinct datasets have been analysed during this study. First, the whole one, including all the nymphs analysed in 16S rRNA gene sequencing, except the only two nymphs collected in November (369 analysed samples). Second, the No *Wolbachia/Arsenophonus* dataset, corresponding to a reduced dataset, where samples harbouring OTUs belonging to *Arsenophonus* and/or *Wolbachia* genera (presenting a number of sequences significantly higher than the highest number of sequences detected for the same OTUs in negative controls, in a 95% confidence interval) were removed (295 samples analysed in this data set). Third, the TBP dataset, composed of samples undoubtedly identified as TBPs positive for only one genus of TBPs (co-infected samples were not included due to the too low number of samples and the difficulty to interpret the results) according to 16S rRNA gene sequencing results as well as microfluidic PCR detection [26]. We considered a sample as positive if a TBP species was previously detected in microfluidic PCR and if at least one OTU identified in 16S rRNA gene sequencing, corresponding to the same pathogenic genera, presented a number of sequences significantly higher than the highest number of sequences detected for this OTU in negative controls (using a 95% confidence interval). Before analysing the data, we applied standard filtering to all the three data sets to remove OTUs associated with weak counts that could hamper the statistical analysis. First, OTUs whose total number of sequences was lower than the total number of samples have been removed. Second, a filter related to the yearly prevalence of each OTU: those

consistently detected in less than 10% of the samples of each year were removed. The application of these filters on the three datasets, the whole, the No *Wolbachia/Arsenophonus* dataset and the TBPs dataset, led to the selection of 89, 82 and 74 OTUs to be included in the subsequent statistical analyses, respectively.

## Statistical analysis

The statistical framework used to describe our data sets is the multivariate Poisson lognormal (PLN) distribution [31]. This statistical distribution is adapted to multivariate count data and shows an expected over-dispersion property compared to the standard Poisson distribution. The main idea behind the PLN model is to represent all the dependency structure between the OTUs in a latent (hidden) multivariate Gaussian layer, while a Poisson distribution in the observation space of the data is used to model counts and noise. Moreover, the PLN framework can be easily interpreted as a Generalised Multi-variate Linear Model (GLM); thus it naturally allows one to include structuring covariates or offsets like in a standard GLM. All the statistical analyses performed in our paper rely on the PLN model and its variants implemented in the R package PLNmodels (version 0.11.0-9005) [32]. In particular, PLNmodels includes variants to perform standard multivariate analyses for count tables, such as PCA (principal component analysis), LDA (linear discriminant analysis) or sparse network reconstruction (aka sparse inference of (inverse) covariance matrix). Additional methodological details can be found in Chiquet et al. [33, 34].

In all the models fitted in this paper, we accounted for the sampling effort (that is, the sequencing depth of each sample) by adding an offset term corresponding to the (log) total sum of counts per sample obtained prior to any filtering. Moreover, in order to correct for any spurious effect induced by a given year of sampling that may hide other effects such as seasonality, we included a covariate to account for the year of sampling in the two first models fitted on both the whole and the No *Wolbachia/Arsenophonus* dataset. While investigating the effect of TBPs on tick microbiota, we also systematically included covariates to account for the year and the season of sampling while dealing with the TBP dataset. A covariate accounting for the presence of TBPs was also considered in the establishment of TBP network in order to evaluate the importance of this variable on tick microbiota network establishment.

Principal component analysis was performed with the PCA variant of PLN and the PLNPCA R function [34], which performs probabilistic Poisson PCA. Network analysis was performed with the PLNnetwork function, which adds a sparsity constraint on the inverse covariance matrix in the latent Gaussian layer [33] to
reconstruct direct associations. Working from the binary adjacency matrix generated from PLNnetwork, we performed a Bernoulli stochastic block model using the BM Bernoulli function of Stochastic Block Model package [35]. This model is a mixture model for binary graph which assumes that nodes belong to different groups according to their connectivity pattern. The number of groups is selected according an integrated completed likelihood criterion. Thanks to this approach, OTUs presenting similar connection profiles were identified and assigned to different clusters.

A differential analysis was also performed using the edgeR package (version 3.30.0) [36] on the TBPs dataset to compare OTU abundances between positive and negative samples. edgeR uses the negative binomial (NB) distribution to model the read counts for each OTU in each sample and computes an empirical Bayes estimate of the NB dispersion parameter for each OTU, with abundance levels specified by a log-linear model. For the same reasons as previously described, we introduced two covariates to control for the year and the season when testing the group main effect (TBP positive or negative). Data were normalised for differences between library sizes using the Trimmed Mean of M value (TMM) method [37]. Models are fitted with the glmFit function, which implements generalised linear methods developed by McCarthy et al. [38]. For each OTU, we tested the group effect using likelihood ratio statistics (glmLRT function). Differentially abundant OTUs were defined as those with p values < 0.05 after adjustment for multiple testing using the Bonferroni procedure.

Networks generated from the TBPs dataset were also compared to each other. For this purpose, a weighted version of Kendall's τ, which integrates the edge appearance rank within families of networks (negative, *Rickettsia*, *Borrelia*, *Anaplasma* and total corrected for TBP effects), was calculated and used to compare them with each other.

All the statistical analyses were performed on R 4.0.2 [39].

## Results

### *Ixodes ricinus* microbiota diversity and composition

Considering the microbiota of all the 371 nymphs, we detected 353 OTUs. Among them, 307 belonged to 109 identified genera, and 46 OTUs belonged to multi-affiliated or unknown genera spread over 15 families. The mean Shannon diversity index is estimated to 2.1 (SD = ± 0.8) and varies between 0.3 (nymph collected in April 2016) and 3.8 (nymph collected in October 2014). Bacterial genera with proportions higher than 0.5% of all sequences in the dataset belonged to *Arsenophonus*, *Candidatus* Midichloria, *Rickettsia*, *Wolbachia*, *Spiroplasma*, *Methylobacterium*, *Mycobacterium*, *Pseudomonas*, *Stenotrophomonas*, *Williamsia*,

*Rickettsiella*, *Chryseobacterium*, *Borrelia*, *Bacillus*, *Anaplasma*, *Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium* and two multi-affiliated OTUs belonging to Rhizobiaceae and Microbacteriaceae families. In total, these sequences represented 93% of all sequences in the dataset (Fig. 1).

Although these genera represent a large number of sequences identified in the dataset, it is important to note that their presence is not equally distributed in all the samples. To investigate this point, we determined, for the five first represented genera in terms of total number of sequences (*Arsenophonus*, *Ca*. Midichloria, *Rickettsia*, *Wolbachia* and *Spiroplasma*), the number of samples where they had a relative abundance of at least 10%. For *Arsenophonus*, while this genus corresponded to 30% of the total number of sequences in the dataset, we detected sequences of this genus in only 16.7% of samples. These 16.7% of samples contained 96.5% of the total number of sequences in the dataset for this genus. We found the endosymbiont *Ca*. Midichloria, in 53% of nymphs which contained 90% of the total number of *Ca*. Midichloria sequences. Concerning *Rickettsia*, only 7.8% of samples concentrated 97% of the total number of sequences. Similarly, 92.1% of sequences belonging to *Wolbachia* were detected in only 12.4% of the samples. Finally, 8.1% of samples harbouring *Spiroplasma* contained 67% of the total number of sequences corresponding to this genus in the dataset.

Because, based on the literature, we hypothesised that *Arsenophonus* and *Wolbachia* OTUs could not be members of the tick microbiota (this hypothesis is developed in the discussion part), we also interested to their distribution in our samples. Indeed, investigating the cumulated percentages of sequences belonging to *Wolbachia* and *Arsenophonus* per sample, we observed that these 2 genera, considered together, represent on average 66.4% of sequences in the 74 infected samples and represent more than 90% of the sequences for 26 of them (data not shown).

## Temporal dynamics of the *Ixodes ricinus* microbiota

### Principal component analysis performed on the whole dataset

The principal component analysis was performed to determine the temporal variation of the tick microbiota, based on the relative abundances of OTUs (Fig. 2). The first two principal components PC1 and PC2 explained 24.04% and 11.16% of the total variance, respectively (Fig. 2A). All samples were clustered into three groups (deduced by graphic reading). The first one was projected in the lower right quarter and was opposed to the rest of the analysed samples, mainly according to the first axis. The small “outsider” cluster was composed of
<figure><img src="image_4.png" /></figure>

nymphs collected during different months and that seems to be randomly projected regarding the month variable. The remaining samples were divided into two clusters which were represented by ticks collected in March/April and May/June/July/August/September respectively. Note that these two clusters partially overlap. Four main genera gathering 17 OTUs seem to drive the formation of the “outsider” cluster: *Wolbachia* (OTUs 1, 2, 3, 4, 11 and 12), *Arsenophonus* (OTUs 1, 2, 3, 5, 6 and 7), *Spiroplasma* (OTUs 1, 3, 4 and 5) and *Pseudomonas* (OTU 2) (Fig. 2B). The identification of the remaining OTUs driving the two other clusters is trickier as they are not distinctly separated. In any way, it seems that OTUs belonging to the Beijerinckiaceae (Beijerinckiaceae_alphal.cluster_1, Beijerinckiaceae_1174–901–12_1, *Methylobacterium*_1-2-3, *Methylorosula* _1, *Methylocella* _1, *Massilia* _1, *Aquabacterium* _1 and *Cupriavidus* _1), Xantobacteriaceae (Xanthobacteriaceae_Multi_1-2), Burkholderiaceae (Burkholderiaceae-UK_1, Multi_1, *Burkholderia-Caballeronia-Paraburkholderia* _1 and *Acidovorax* _1), Rhizobiaceae (Rhizobiaceae-Multi_1 and *Allo-Neo-Para-Rhizobium* _1) and Microbacteriaceae (Microbacteriaceae-Multi_1, *Amnibacterium* _1) families as well as OTUs *Actinomycetospora* _1, *Williamsia* _1, *Pseudomonas* _1-4-7, *Luteibacter* _4, *Sphingomonas* _1-3-4-6, *Mycobacterium* _1-3 and *Kineococcus* _1 strongly explained the variance along the negative part of the second axis.

Principal component analysis performed on the No
Wolbachia/Arsenophonus dataset

This second principal component analysis was per-
formed to determine the temporal variation of the tick
microbiota excluding the *Wolbachia*- and *Arsenopho-
nus*-positive samples (Additional file 2) as we hypothe-
sised that these two genera are not real members of the
tick microbiota (this hypothesis is developed in the dis-
cussion part). Here, the first two principal components
PC1 and PC2 explained 19.47% and 11.05% of the total
variance, respectively (Fig. 3A). The clustering according
to months is still visible with samples clustered into
three groups (deducted by graphic reading): the cluster 1
mainly composed with samples collected in February/
March, the cluster 2 with samples collected in April and
the third cluster regrouping mainly samples collected in
May/June/July/August/September. Samples collected in
October seemed to be distributed in both clusters 1 and
3 (mainly on the left part of the plot, separated from the
rest of the community *via* the first axis), while those col-
lected in January do not seem to follow any particular
distribution. Ticks collected in April (cluster 2) are dis-
tributed all along the second axis, but only on the right
part of the plot, and therefore seemed to be separated
from the rest of the community mainly *via* the first axis.
By contrast, ticks in cluster 1, mainly collected in Febru-
ary/March, seemed to be distributed all along the first
axis, but mainly on the bottom of the plot, and are
<figure><img src="image_5.png" /></figure>

consequently separated from the rest of the community
mainly through the second axis. For the cluster 3, all
samples seemed to be distributed on the top left corner
of the plot, distinct from the two other clusters 1 and 2
mainly through the first and second axis, respectively.

Looking at the correlation circles (Fig. 3B), we can see
that the positive part of the first axis (where are mainly dis-
tributed nymphs sampled in April), seemed to be mainly ex-
plained by *Rickettsia_3* and *Pseudomonas_5*. To note that
their effect is not as strong as those observed on the
<figure><img src="image_6.png" /></figure>

negative part of the same axis and implicating mainly OTUs belonging to the families Microbacteriaceae (Microbacteriaceae-Multi_1 and *Amibacterium_1*), Beijerinckiaceae

(Beijerinckiaceae_1174-901-12_1-2-3, Beijerinckiaceae_1-2, *Methylobacterium_1-2-3-4-5-6-7*, *Methylorosula_1-2, *Methylocella_1*), Burkholderiaceae (Burkholderiaceae_
UK_1) and Xhantobacteriaceae (Xhantobacteriaceae-Multi_2), as well as the OTUs *Sphingomonas_1-4*, *Williamsia_1*, *Actinomycetospora_1*, *Jatrophihabitans_2*, *Acidiphilium-1* and *Kineococcus_1*. The negative part of the second axis (where samples of cluster 1 are mainly distributed), seemed to be mainly explained by *Mycobacterium_1-5-6*, *Isosphaeraceae-UK_1*, *Halomonas_1*, *Rhodanobacter_1* and *Luteibacter_1*.

### Ixodes ricinus microbiota correlations

Thanks to the network analysis performed on the whole dataset, we observed a total of 224 significant partial correlations between 89 OTUs. Interestingly, 97.8% of these

partial correlations were positive (Fig. 4). Positive partial correlations frequently occurred between OTUs belonging to the same genus or family, as it is the case for OTUs belonging to: *Wolbachia*, *Arsenophonus*, *Spiroplasma*, *Ca*. Midichloria, *Rickettsia*, *Mycobacterium*, *Sphingomonas* and Beijerinckiaceae (including Beijerinckiaceae_1174-901-12, Beijerinckiaceae_alphal.cluster, *Methylobacterium*, *Methylorosula* and *Methylocella*). By contrast, this pattern was not observed between any of the *Pseudomonas* OTUs.

In the network, five main clusters of OTUs presenting similar connection profiles were detected (Additional file 3). Clusters 1 and 5 seem to be composed of

<figure><img src="image_7.png" /></figure>
environmental OTUs. Cluster 1 contained ten OTUs (Rhizobiaceae-Multi_1; Allo-Neo-Para-Rhizobium_1; *Pseudomonas* 7; *Stenotrophomonas* 1, *Luteibacter* 1; *Acidovorax* 1; *Burkholderiaceae*-Multi_1; *Aquabacterium* 1; *Cupriavidus* 1; *Pedobacter* 1), including five belonging to the Betaproteobacteriales while other OTUs belong to four other orders: Sphingobacterales Xanthomonadales Rhizobiales and Pseudomonadales. Cluster 5 is composed of 17 OTUs (*Acidiphilium* 1; *Methylobacterium* 6; *Methylobacterium* 1; *Methylobacterium* 7; *Methylobacterium* 3; *Methylobacterium* 5; *Methylobacterium* 2; *Methylobacterium* 4; *Methyloresula* 1; Beijerinckiaceae-alpha1.cluster_1; Beijerinckiaceae_1174-901-12_2; Beijerinckiaceae_1174-901-12_1; Beijerinckiaceae_1174-901-12_3; *Methylocella* 1; Burkholderiaceae-UK_1; *Amnibacterium* 1; *Kineococcus* 1), mainly belonging to the Rhizobiales order (13 OTUs). Other OTUs of this cluster belong to Betaproteobacteriales, Micrococcales, Kineosporiales and Acetobacterales orders. Two other clusters are mainly composed of OTUs belonging to symbiotic or pathogenic genera (clusters 2 and 4). Cluster 2 is composed of 6 OTUs (*Ca*. Midichloria_4; *Ca*. Midichloria_3; *Rickettsia* 5; *Rickettsia* 1; *Rickettsia* 3; *Pseudomonas* 5) mainly belonging to the Rickettsiales order, and more precisely to the genera *Ca*. Midichloria and *Rickettsia*, while one of them belong to the Pseudomonadales order, more precisely to the *Pseudomonas* genus. Cluster 4 is composed of 17 OTUs (*Wolbachia* 11; *Wolbachia* 12; *Wolbachia* 2; *Wolbachia* 4; *Wolbachia* 3; *Wolbachia* 1; *Pseudomonas* 2; *Arsenophonus* 7; *Arsenophonus* 6; *Arsenophonus* 3; *Arsenophonus* 5; *Arsenophonus* 1; *Arsenophonus* 2; *Spiroplasma* 1; *Spiroplasma* 3; *Spiroplasma* 5; *Spiroplasma* 4), including six belonging to the *Wolbachia* genera of the Rickettsiales order, six to the *Arsenophonus* genus of the Enterobacterales, four to the *Spiroplasma* genus of the Entomoplasmatales order and one to the *Pseudomonas* genus of the Pseudomonadales order. Cluster 3 is the largest with 39 OTUs that have in common that they have little or no correlations with other OTUs. Finally, although the following OTUs do not belong to the same clusters, negative partial correlations were observed between *Pseudomonas* 1 and both *Arsenophonus* 5 and *Wolbachia* 3. *Wolbachia* 3 was also negatively correlated to *Pseudomonas* 4, *Mycobacterium* 1 and *Halomonas* 1.

## Links between the presence of tick-borne pathogens and the *Ixodes ricinus* microbiota

### Microbiota structure comparison between TBP-positive and TBP-negative samples

As expected, the comparison of OTU abundances between each group of TBP-positive samples (*Rickettsia*, *Borrelia* and *Anaplasma*) and TBP-negative samples allowed to detect a significantly higher abundance of

OTUs belonging to the genus of the tested TBP (Additional files 4, 5, 6). Concerning *Rickettsia*-positive samples, seven other OTUs were also harbouring significantly different abundances compared to negative samples: *Ca*. Midichloria_3, *Pseudomonas* 3-5-8, *Bacillus* 1 and Rhizobiaceae-Multi_1 were significantly more abundant in *Rickettsia*-positive samples while *Spiroplasma* 1 was significantly less abundant. In *Borrelia*- and *Anaplasma*-positive samples, no other OTUs than those belonging to the corresponding tested genera were significantly over- or underrepresented.

### Microbial network comparison between TBP-positive and TBP-negative samples

The network analysis was performed on 5 different datasets. First on samples positive for *Rickettsia*, *Borrelia* or *Anaplasma*; then on samples negative for all these genera and finally on all the samples included in the above-mentioned datasets, but considering a covariate correcting for the effect of TBP presence. As performed for the total network, OTUs presenting similar connection profiles were investigated *via* SBM. However, due to the high parsimony level of these networks, this analysis could only distinguish at best two clusters: connected samples from those presenting no correlations. Nonetheless, the determination of Kendall's τ, integrating the edge appearance rank for all the calculated networks, allowed us to compare them and observe significant differences between all the obtained networks (Additional file 7).

In the network analysis performed on positive ticks for *Rickettsia*, we observed 17 significant partial correlations, five were negative (Fig. 5A). Several members of the *Rickettsia* genus were positively correlated to OTUs belonging to *Ca*. Midichloria (*Rickettsia* 5/*Ca*. Midichloria_3) and *Pseudomonas* (*Rickettsia* 1/*Pseudomonas* 3) genera, as already observed in the total dataset network. *Rickettsia* 1, *Pseudomonas* 3 and *Ca*. Midichloria_3 were negatively correlated to *Bacillus* 1. This latter OTU was also negatively correlated with 2 OTUs belonging to environmental genera, *Burkholderia-Caballeronia-Paraburkholderia* 1 and *Chryseobacterium* 1 (that was also positively correlated to *Ca*. Midichloria_3), and positively correlated with *Anaplasma* 1. *Bacillus* 1 therefore appeared as a key member of this network exhibiting several partial correlations with environmental bacteria but particularly with some pathogenic and symbiotic genera. Several positive partial correlations were also observed between OTUs belonging to environmental genera, linking several OTUs within their own genus (*Methylobacterium* and *Mycobacterium*) or between several different genera or families (Beijerinckiaceae and *Aquabacterium*; *Stenotrophomonas* and *Acinetobacter* as well as
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Network analysis. Representation of the significant partial correlations detected between OTUs of the TBP dataset. A Considering only samples positive for *Rickettsia*, B considering only samples positive for *Borrelia*, C considering only samples positive for *Anaplasma*, D considering only negative samples and E considering positive and negative samples as well as a covariate accounting for the presence of TBPs. OTU circles are colored by taxonomic order. These circles represent nodes of the networks. Their size is proportional to the sum of the incoming edge weights. Thickness of the edge is proportional to the strength of the observed partial correlation. Positive partial correlations are represented by red edges, negative partial correlations are represented by turquoise edges. For representation reasons, the OTU *Burkholderia-Caballeronia-Paraburkholderia_1* was abbreviated to B-C-P_1</figcaption></figure>

Microbacteriaceae and *Amnibacterium*). We then performed a network analysis considering only ticks positive for *Borrelia* (Fig. 5B). Only 8 significant partial correlations were observed, including 2 negatives. Negative partial correlations were observed between first *Borrelia*-

RF_1 (responsible of Relapsing Fever) and *Borrelia*-LB_1 (responsible of Lyme Borreliosis) and then between *Borrelia*-RF_1 and *Rickettsiella_2* (Fig. 5B). A positive partial correlation implying *Rickettsiella_2* and *Bacillus_1*, and another one implying *Spiroplasma_1* and 2 were
also observed as well as several positive partial correla-
tions implying environmental genera (Bacillus_1, Mas-
silia_1, Aquabacterium_1 and Mycobacterium_2, as well
as Mycobacterium_5 and Luteibacter_1). In the network
analysis performed on ticks only positive for Anaplasma,
we observed 23 significant partial correlations with 2
identified as negative (Fig. 5C). A positive partial correl-
ation was observed between Rickettsiella_2 and Parra-
coccus_1. Spiroplasma_2 was positively correlated with
Jatrophihabitans_2 and several members of the Rhizo-
biales order (Methylobacterium_3, Methylocella_1 and
Xhantobacteraceae_Multi_1). Other partial correlations
concerned environmental bacteria, with members of the
Rhizobiales order, correlated with each other or with
OTUs belonging to environmental genera corresponding
to other orders (Acetobacterales, Bacillales, Betaproteo-
bacterales, Corynebacterales, Micrococcales and
Xanthomonadales). No correlations were observed be-
tween Anaplasma and the other members of the bacter-
ial community. The negative network (Fig. 5D),
comprising only free pathogen samples, was composed
of 33 positive partial correlations and no negative correl-
ation. In this network, a strong partial correlation was
observed between Spiroplasma_1 and 2. Most of the
other correlations were observed between members of
the Rhizobiales order, mainly between each other, but
also with environmental OTUs belonging to other or-
ders, such as Acetobacterales, Betaproteobacterales,
Frankiales, Micrococcales, Pseudonocardiales and Sphin-
gomonadales. Although it contains less partial correla-
tions (25), the total network corrected for TBP effects
(Fig. 5E) is very similar to the negative network. Indeed,
partial correlations observed between Spiroplasma_1
and 2, as well as most partial correlations involving Rhi-
zobiales order between each other and with OTUs be-
longing to other orders are also observed in this
network.

**Discussion**

Ixodes ricinus microbiota diversity and composition

Investigating the microbiota of 371 *I. ricinus* nymphs, we identified 109 bacterial genera, i.e. within the range of previous observations in this tick species [6, 12]. Similarly, the mean Shannon diversity index (=2.1) [22] is in the range reported in the literature for *Ixodes* ticks [8, 20, 40, 41]. However, these values are known to fluctuate, mainly according to the tick instars, species or localization, and are therefore difficult to compare [5, 8, 10, 12, 17, 41–45]. Furthermore, not all these studies used negative controls to identify potential contaminating OTUs and remove them from their datasets, thus calling for caution when attempting to draw conclusions concerning these differences, particularly knowing that such OTUs can represent more than 50% of the

sequences detected in tick samples [22]. In our dataset,
three quarters of the sequences belonged to five genera:
*Arsenophonus*, *Ca*. Midichloria, *Rickettsia*, *Wolbachia*,
*Spiroplasma*, all of which corresponded to well-known
maternally inherited bacteria in arthropods [46]. Their
high prevalence in the dataset suggest they are widely
distributed in a large proportion of tick samples and was
clearly the case in *Ca*. Midichloria, in which most of the
detected sequences were found in half the tick samples.
This result is in line with previous reports that almost
100% of *I. ricinus* females and larvae, and almost 50% of
nymphs and males were infected by this endosymbiont
[47–50]. The four other genera (*Rickettsia*, *Spiroplasma*,
*Wolbachia* and *Arsenophonus*) were mainly detected in
less than 17% of samples (7.8%, 8.1%, 12.4% and 16.7%,
respectively). The contrast between the high proportions
of sequences corresponding to these genera and the
small number of samples infected is due to their high
relative abundance in infected samples. Focusing on
*Rickettsia*, the sequences corresponded to *R. helvetica*,
which have already been detected in these samples [26]
and known to exhibit a high bacterial load in ticks [51,
52]. Concerning *Arsenophonus* and *Wolbachia*, while the
high number of sequences detected in the dataset as well
as their high cumulated percentage in infected samples
suggest they are dominant members of the *I. ricinus*
microbiota, their detection could be due to the presence
of the parasitoid wasp, *Ixodiphagus hookeri*, in ticks [53,
54]. This high prevalence, almost excluding endogenous
tick bacteria, prompted us to wonder about the potential
role of this symbiont in the colonisation of *I. ricinus* by
*I. hookeri*. Indeed, as previously hypothesised by Plantard
et al. [53], the presence of *Wolbachia* could induce an
immune reaction in the tick that could have a negative
impact on other microorganisms hosted by the tick. This
kind of strategy, dealing with the host immune system,
could allow the symbiont proliferation and so the egg
maintenance, as it has already been observed in *Dros-
ophila stimulans*, in which the presence of *Wolbachia*
reduced the ability of the host immune system to encap-
sulate the eggs of the parasitoid wasp *Leptopilina hetero-
toma* [55]. Finally, among the most frequently detected
genera in the *I. ricinus* microbiota, *Spiroplasma* is a ma-
ternally inherited symbiont particularly well known in
arthropods, including ticks [56]. With most of their se-
quences detected in only 8% of tick samples, this bacter-
ial genus probably represents secondary symbionts [56],
which are not essential for the survival of their host but
may confer conditional adaptive advantages, such as pro-
tecting their host against pathogens and natural enemies
[57, 58]. This is probably the case for *I. ricinus* ticks po-
tentially infected with the parasitoid wasp *Ixodiphagus
hookeri*, as we discuss later in this section. Among the
remaining genera, several are known to circulate in *I*.
ricinus ticks, including *Rickettsiella* (a maternally inherited bacterium in arthropods), *Borrelia* (Lyme Borreliosis and Relapsing Fever) and *Anaplasma*. Some others, including *Methylobacterium*, *Mycobacterium*, *Pseudomonas*, *Stenotrophomonas*, *Williamsia*, *Bacillus*, *Chryseobacterium*, *Acinetobacter* and *Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium*, as well as multi-affiliated OTUs belonging to Rhizobiaceae and Microbacteriaceae, which have already been detected in *Ixodes ricinus* ticks [6, 12, 16, 17, 19], and have also been identified in the environment [59–64]. While very few of the abovementioned references reported performing negative controls while studying tick microbiota [12, 16], it is important to keep in mind that some of these genera could also be contaminants that arise during the extraction or amplification steps [22, 30, 65–67]. Our dataset was subject to a thorough contaminant filtering, based on negative control composition [22], thereby reducing the risk that these OTUs correspond to contaminants. Nonetheless, the question remains as to whether these environmental taxa belong to the internal tick microbiota or to the cuticular microbiome? Probably to both. Indeed, Binetruy et al. [68] demonstrated that surface sterilisation of the tick using ethanol (which we used in the present study) was not as efficient as bleach at removing bacterial DNA from the tick cuticle. These authors also showed that environmental taxa, only detected in ethanol cleaned whole ticks, belonged to the same family as those directly detected on the tick cuticle by swabbing. However, they also found some of these taxa in the gut of ticks. While studying the microbiota in organs of *I. ricinus* ticks previously cleaned with bleach, Hernández-Jarguín [19] also detected bacterial genera commonly found in the environment.

## Temporal variability of *Ixodes ricinus* microbiota

We first performed a principal component analysis of the whole dataset to characterise the temporal dynamics of the *I. ricinus* microbiota. Our results provide evidence that the *I. ricinus* microbial communities vary markedly, with two main clusters, the first grouping the microbial communities identified in March–April and the second grouping the microbial communities identified from May to September. Interestingly, this temporal pattern was observed for samples belonging to each of the 3 years, suggesting predictable temporal differences in the *I. ricinus* microbial community structure. A third cluster was also identified and corresponded to nymphs collected across the 3 years. This “outsider” cluster was mainly driven by members of the genera *Wolbachia*, *Arsenophonus* and *Spiroplasma*. As mentioned above, the presence of *Wolbachia* and *Arsenophonus* species in *I. ricinus* ticks is probably due to the presence of the parasitoid wasps, *Ixodiphagus hookeri* [53, 54], which

depends entirely on the tick for its development. In their study, Plantard et al. [53] reported that almost 100% of *I. hookeri* were infested by *Wolbachia pipientis*. These authors also showed that all unfed tick nymphs parasitised by *I. hookeri* also harboured *Wolbachia*, while (with only one exception) non-infected ticks were *Wolbachia*-free. *Wolbachia* was shown to be vertically transmitted from the female *I. hookeri* to its eggs and the subsequent generation harboured *Wolbachia*. Similarly, Bohacsova et al. [54] detected the symbiont *Arsenophonus nasoniae* only in nymphs infected by the wasp *I. hookeri*, and almost 30% of the wasp population was infested by *A. nasoniae*. In this context, through the study of tick microbial communities, we suggest that identifying the temporal dynamics of both *Wolbachia* and *Arsenophonus* could serve as a proxy to characterise both the infection rate and temporal dynamics of *I. hookeri* in *I. ricinus*, even if it would probably be easier and more efficient to use primers or probes that specifically match the parasitoid DNA. While the presence of *Wolbachia* and *Arsenophonus* is thus probably linked to the presence of *I. hookeri* in ticks, we hypothesised that their presence in our dataset could alter our characterisation of the temporal dynamics of the *I. ricinus* microbial communities. We therefore chose to perform another PCA analysis after removing all samples harbouring one or both genera *Wolbachia* and *Arsenophonus*. As expected, the new PCA performed on the “reduced” dataset improved separation of the clusters. We finally distinguished three main clusters. Whatever the sampling year, this analysis mainly distinguished ticks collected in February/March (cluster 1), from those sampled in April (cluster 2) or from those sampled from May to September (cluster 3), with ticks sampled in October distributed in both clusters 1 and 3, suggesting that the factors which drove the tick microbial communities were probably the same in the three consecutive years. Interestingly, we also observed that the OTUs that best explained the dynamics of the tick microbial community structure were neither tick symbionts nor pathogens but were genera or families that are characteristic of the environment (i.e. *Sphingomonas*, *Williamsia*, *Actinomycetospora*, *Jatrophihabitans*, *Acidiphilium*, *Kineococcus*, *Mycobacterium*, *Halomonas*, *Microbacteriaceae*, *Beijerinckiaceae*, *Xanthobacteriaceae*, *Burkholderiaceae*, *Rhodanobacteraceae* and *Isosphareaceae*). Considering that these OTUs belong to the cuticular microbiome (as discussed above), the fluctuations we observed could therefore correspond to variations in the environmental microbiota [69]. However, these OTUs could also correspond to bacteria acquired by ticks from the environment during water uptake for example, or through the tick’s other orifices such as spiracles and the anal port. This was suggested by Narasimhan and Fikrig [70] who observed that ticks
hatched in a sterile environment harboured a signifi-
cantly different gut microbiota than those hatched in
"normal" conditions [71]. We hypothesise that the tem-
poral variations we observed in our analysis could thus
also be due to differences in the composition of the en-
vironmental microbiota. Otherwise, abiotic factors such
as temperature, which are likely to follow the same
yearly pattern, have been shown to influence the diver-
sity of the microbiota of *I. scapularis* ticks [45]. Further-
more, blood meal and host identity have also been
shown to influence the diversity of tick microbiota [11,
40, 71, 72]. Finally, considering the life cycle of ticks,
one other hypothesis which could explain this pattern is
that tick microbiota structure and composition could be
linked to the feeding status of ticks. While quantifying
the lipid content of *Ixodes ricinus* ticks sampled monthly
over a period of several years, Randolph et al. [73] and
Abdullah et al. [74] managed to discriminate tick feeding
cohorts questing at different periods of the year. Both
studies reported that lipid reserves were higher in ticks
collected at the end of the year than in those sampled
from the end of spring and throughout summer. Ac-
cording to Abdullah et al. [74], ticks sampled at the end
of the year (high lipid reserves), correspond to those that
succeeded in obtaining their blood meal in the previous
spring and emerged in autumn. The ticks sampled after
the end of the year and at the beginning of the following
spring, presenting a relatively high lipidic reserves that
trend to increase towards March, were hypothesised to
correspond to a mix of tick populations: (1) those just
emerging, with high lipid reserves and (2) those that
already emerged in autumn, with lower lipid levels. Ticks
still questing at the end of spring or during summer
were hypothesised to correspond to those that failed to
find a host and had almost exhausted their lipid reserves.
In this case, one could hypothesise that the ticks in clus-
ter 3 (mainly ticks sampled from May to September)
correspond to those with very few energy reserves, while
ticks in cluster 1 (sampled in February/March) corres-
pond to those with higher energy reserves. Ticks sam-
pled in October, distributed in clusters 1 and 3, could
correspond to a mix of two populations: (1) those with
high lipid reserves (i.e. which had fed the previous
spring and had just emerged) and (2) those whose energy re-
serves were almost exhausted (which had fed the previ-
ous year), as observed for nymphs sampled in October
2015 and 2016 in Abdullah et al. [74]. However, this hy-
pothesis does not elucidate the clustering of the ticks
sampled in April, thus suggesting that other factors (re-
lated to the life cycle of ticks or their environment),
could be at the origin of this pattern. One can further
hypothesise that the community patterns of bacteria are
also related to physiological stress of the tick, the micro-
biota "supporting" the tick's physiological requirements,

perhaps in agreement with differences in nutritional
stress, as suggested above, or perhaps to abiotic stress
linked to temperature or saturation deficit. In any case,
further studies will be necessary to elucidate this point
and confirm or invalidate these hypotheses.

*Ixodes ricinus* microbial community interactions

Identifying and understanding tick-borne microbe inter-
actions is a precondition for the development of new
strategies to control ticks and tick-borne diseases using
the tick microbiome. Using partial correlation networks,
we assessed which tick bacterial OTUs were correlated
in order to identify possible associations between tick
microbial community members.

First, used on the whole dataset, network analysis re-
vealed that more than 97% of detected interactions be-
tween members of the *I. ricinus* microbial community
were positive. Taxa with positive associations have usu-
ally been interpreted as functional guilds of organisms
performing similar or complementary functions [75, 76]
or featuring interactions shaped by interspecies cross-
feeding [77], although sometimes they may mainly re-
flect shared habitat preferences [78]. Similarly, negative
associations may reflect interactions including competi-
tion and niche partitioning. The vast majority of positive
correlations observed therefore suggest that tick micro-
bial communities favour mutualism and perform similar
or complementary functions. But let us keep in mind
that the correlations observed in our study were ob-
tained by examining entire individuals and may impair
the observation of associations at a finer scale (e.g. or-
gans) [27].

Due to the low infection rate observed in the nymphs
studied [26], the low proportion of ticks in our dataset
that tested positive for pathogens could bias the analyses
and conceal important information such as crucial inter-
actions between TBPs and members of the tick micro-
biota. To overcome this bias, we decided to perform
supplemental analysis to compare TBP-positive and TBP-
negative samples. Further, to get rid of effects that could
hamper or skew the observation of the impact of TBPs,
samples that appeared to be parasitised by *I. hookerii* (i.e.
with numerous sequences belonging to *Arsenophonus
and/or Wolbachia*) were excluded from this analysis, and
a covariate correcting for the effect of the sampling season
was added, compared to the analysis performed on the
whole dataset. The resulting negative network was mainly
composed of partial correlations linking environmental
OTUs, and very few correlations remained compared to
the network obtained using the whole dataset. These re-
sults suggest that most correlations previously observed in
the total network are linked to the presence of TBPs, the
response to the parasitoid or the seasonal effect. This fur-
ther strengthens the importance of these variables in the
tick microbial community structure and underlines the
adaptability of the tick microbiota to variable conditions.
Comparing negative samples with those infected by patho-
gens enabled us to demonstrate that with TBPs, the struc-
ture and correlations of tick microbial communities were
considerably modified in presence of TBPs and according
to the concerned pathogen. Indeed, in ticks infected by
*Rickettsia*, the proportions of several OTUs (i.e. *Spiro-
plasma* symbiont and environmental OTUs) were signifi-
cantly affected. Moreover, more negative correlations
were detected in TBP samples, suggesting much more
competition between members of the tick microbial com-
munities in presence of TBPs. In addition, in the “free
pathogen” network, we observed a lot of correlations in-
volving environmental OTUs (i.e. belonging to Rhizo-
biales, Acetobacterales, Betaproteobacterales, Frankiales,
Micrococcales, Pseudonocardiales and Sphingomonadales
orders), whereas these correlations decreased or even dis-
appeared in infected ticks, suggesting that the presence of
TBPs affects these microbial interactions. Two main hy-
potheses emerged from these results: the tick microbiota
are initially disturbed and thus favour infection by patho-
gens, or the presence of pathogens affects the structure of
the tick microbiota. Because the structure and interactions
of tick microbial communities are completely different de-
pending on the pathogen concerned, we first suggest that
the presence of pathogens in ticks is likely to affect the
other members of the tick microbiota. However, disturb-
ance of the microbiota could facilitate the installation of
pathogens in ticks. To give an example, the presence of
*Anaplasma phagocytophilum* in *I. scapularis* could modify
the gut microbiota whose modification could, in turn, dis-
rupt the integrity of the gut membrane and facilitate entry
by the pathogen [79]. Based on our data, we are currently
unable to conclude and further experiments are required
to this end.

Interactions between non-pathogenic OTUs

First, we detected in the a cluster of OTUs (cluster 4)
grouping several OTUs belonging to *Wolbachia*, *Arsen-
phonus* and *Spiroplasma* OTUs and one *Pseudomonas*
OTU (*Pseudomonas_2*), positively linked to each other
in the whole network. As previously mentioned, detect-
ing *Wolbachia* and *Arsenophonus* in *I. ricinus* is prob-
ably evidence for the presence of the parasitoid *I.
hookeri* in these ticks [53, 54]. While the presence of the
deer-associated *Anaplasma phagocytophilum* in ticks has
been previously shown to be positively correlated with
the presence of *I. hookeri*, due to the way of life/mode of
hunting of this parasitoid [80], we found no correlation
between this pathogen and *Wolbachia* or *Arsenophonus*
OTUs in the present study. Conversely, *Spiroplasma*,
identified as arthropod symbionts including ticks [46, 56,
81–84], were closely linked to the dynamics of these two

genera. Identified as a “male killer” in many other ar-
thropods [85–90], our results on nymphs alone do not
allow us to conclude on this potential role. Nonetheless,
assessing the proportions of contaminant sequences
when characterising both nymph and adult tick micro-
bial communities with high-throughput sequencing, we
detected a high proportion of Spiroplasma sequences in
males [22], suggesting that these Spiroplasma species do
not cause specific male mortality in I. ricinus ticks. Fur-
thermore, it should be noted that no evidence of sex-
ratio distortion was found by Binetruy et al. in the popu-
lation of Rhipicephalus decoloratus ticks infected with
Spiroplasma ixodetis [56]. Assuming that the detection
of both Wolbachia and Arsenophonus is probably due to
the infection of ticks by I. hookeri, we hypothesise that
Spiroplasma is upregulated by the presence of the para-
sitoid in ticks and may be a defensive response mechan-
ism against I. hookeri, as previously reported in
Drosophila melanogaster infested by two species of para-
sitoid wasps [91]. Nevertheless, we cannot rule out the
hypothesis that, while identified in a large variety of tick
species [56], Spiroplasma may correspond to a symbiont
of I. hookeri. Besides, while Pseudomonas_2 was identi-
fied in this cluster and was positively linked to its other
members, it should be noted that two other Pseudo-
monas OTUs (Pseudomonas_1 and 4) were negatively
correlated with several Wolbachia and Arsenophonus
OTUs. In our dataset, very few OTUs belonging to the
same genus displayed contrasting patterns of correlation,
but Pseudomonas is an exceptionally versatile genus in-
cluding plant pathogens, human and animal opportunis-
tic pathogens, and saprophytic bacteria found in water
and soil exhibiting great adaptation to their environment
[92]. Notably, some members of this genus have been re-
ported to influence arthropod survival. For example, P.
fluorescens strains have been reported to confer better
survival to the Varroa destructor mite and the Galleria
mellonella waxworm in presence of a fungal pathogen of
Beauveria bassiana arthropods [93, 94]. Conversely, P.
entomophila has been reported to be entomopathogenic
for Drosophila [95, 96]. Such contrasted behaviour in ar-
thropods could explain the contrasting correlations ob-
served with Pseudomonas OTUs, and suggest that, like
Spiroplasma OTUs, some of them may be involved in
defense mechanisms against I. hookerii.

Two other groups of OTUs, comprising non-
pathogenic OTUs, were identified in the total network,
and correlations involving several members of these
clusters were still observed in the negative and covariate
corrected networks. These were the OTUs belonging to
clusters 1 and 5 that belong to the Betaproteobacterales,
Sphingobacterales, Xanthomonadales, Rhizobiales and
Pseudomonadales orders in the case of cluster 1, and of
OTUs belonging to Rhizobiales, Betaproteobacterales,
Micrococcales, Kineosporiales and Acetobacterales orders in the case of cluster 5. All these OTUs belong to genera and families of microorganisms commonly found in the environment [62, 97–102], and known to potentially infect ticks. The high connectivity displayed by these OTUs may reflect similar functions between OTUs adapted to contrasting environmental conditions (functional redundancy). However, such a scenario would have probably resulted in many more negative correlations between these OTUs. On the other hand, these findings may indicate reduced functional capabilities of each OTU and hence a complementary functional role for these taxa in ticks. All these hypotheses need to be investigated in future works on the functional ecology of tick microbiota.

Finally, two correlations linking OTUs belonging to the symbiotic genus *Rickettsiella* to OTUs belonging to the environmental genera *Paracoccus* and *Bacillus* were observed in the *Anaplasma* and *Borrelia* networks, respectively. While the *Rickettsiella* genus has been detected in several species of ticks [46], notably in *I. ricinus* ticks [12, 17–19, 46], little is known about the importance of this maternally inherited bacterium in tick species. Its correlations with environmental OTUs in the presence of TBPs could mean it plays a particular role in the pathogenic context in *Ixodes ricinus* ticks.

## Interactions between tick-borne microbiota and pathogens

Interestingly, we found correlations between *Ca. Midichloria* and *Rickettsia* OTUs that were strong enough to be observed not only in the *Rickettsia* network, but also in the whole dataset network where these OTUs were grouped in cluster 2. While our sequencing approach does not allow us to identify *Rickettsia* at the species level, we hypothesise that they correspond to pathogenic agents, as *Rickettsia helvetica* was specifically detected in these same samples [26]. In addition, *Ca. Midichloria* abundance was significantly higher in *Rickettsia*-positive samples than in the negative ones. Considering the strong prevalence of *Ca. Midichloria* in *I. ricinus* ticks, the positive relationship between members of this genus and members of *Rickettsia* genus suggests a facilitating role for *Ca. Midichloria* in *I. ricinus* colonisation by *Rickettsia*. This hypothesis is in accordance with the results obtained by Budachetri et al. [28] who observed a positive correlation between *Ca. Midichloria mitochondrii* load and *Rickettsia parkeri* presence in the tick *Amblyomma maculatum*. The fact that *Rickettsia* OTUs were positively correlated with maternally inherited bacteria is particularly surprising, especially because members of this genus have been frequently reported to be involved in antagonist relationships with symbiotic or pathogenic genera [9, 10, 16, 103–107]. The potential

complementarity of these two genera should be examined in more detail in the future by characterising the bacterial transcriptome or metabolome of ticks infected with these bacteria alone or in association. Interestingly, the OTU of *Borrelia Relapsing fever* (*Borrelia RF*) was negatively correlated with *Rickettsiella* and *Borrelia Lyme Borreliosis* OTUs (*Borrelia LB*). The negative correlations between both *Borrelia* groups (relapsing fever vs Lyme borreliosis) we identified using the 16S rDNA are in full agreement with previous results we obtained using the same tick samples and using a high-throughput microfluidic real-time PCR with specific primers and probes to detect specifically tick-borne pathogens [26]. Even though this association was not significantly underrepresented compared to a random distribution, our previous results highlighted higher prevalences of *Borrelia burgdorferi s.l.* when *B. miyamotoi* prevalences were low, and vice versa. Competition for the same niche may explain the negative correlations observed between these two groups of *Borrelia*. However, these findings contrast with those reported by Aivelo et al. [16] as these authors found positive correlations between the relapsing fever spirochete *B. miyamotoi* with *Rickettsiella* and two *Borrelia* species belonging to the group *Borrelia* Lyme borreliosis (*B. garinii*, *B. afzelii*). Because clinical co-infections with several TBPs are commonly reported [108–110] and are known to affect both disease symptoms and severity [111, 112], it is now crucial to identify the conditions in which *Borrelia* RF and *Borrelia* LB are in competition, or on the contrary, in which they could collaborate and thus co-infect ticks. Furthermore, and contrary to what it has been observed on *Ixodes scapularis* in the literature [113, 114], we did not detect significant differences in the structure of the tick microbiota in *Borrelia*-positive samples compared to the negative ones, with no OTUs under- or overrepresented in these samples. The difference of tick species investigated could explain such contrasting results, especially since most of the genera overrepresented in the Brinkerhoff study [114] were not detected in our study. Furthermore, Chauhan et al. [113] focused their analysis on *Borrelia burgdorferi* s.l. investigated in females, while we investigated *Borrelia* (including Lyme borreliosis and relapsing fever *Borrelia*)-positive nymphs that could be another explanation for such differences of observation in our study. Otherwise, we also observed several partial correlations between pathogens and OTUs usually known to be “environmental” bacteria. This was notably the case of a *Pseudomonas* OTU positively correlated with *Rickettsia*. In addition, this OTU and two others belonging to the same genus were more abundant in *Rickettsia*-positive samples. While several *Pseudomonas* OTUs were previously identified as contaminants and removed from the
dataset analysed [22], those remaining are involved in
several interactions with different members of the tick
microbial community, such as *Wolbachia* and *Arsenophonus*, as already discussed above, but also with TBPs,
demonstrating the versatility of members of this genus
and their importance in the tick microbiota structure.
*Bacillus* also appears to be a key member of the *Ixodes ricinus* microbial community linked to the presence of
TBPs. While Adegoke et al. [115] observed higher abundances of this bacterial genus when the parasite *Theileria* was present in the tick *Rhipicephalus microplus*, our findings demonstrate a positive correlation between *Bacillus* and *Anaplasma* and a negative one with *Rickettsia*. The latter, implying *Bacillus* and *Rickettsia*, suggest potential competition between these bacteria and may be a first step to develop a future tool to control tick infections by *Rickettsia*. Furthermore, because the dynamics of environmental bacteria found in ticks varies over the year, probably due to contrasting environmental conditions, the vegetation and the tick hosts, their positive or negative interactions with pathogens suggest that their presence or absence is an important factor to take into account to better understand the temporal dynamics of TBPs. Finally, all these correlations involving “environmental” OTUs suggest that their detection in tick microbiota is probably not only the result of accidental ingestion, but more likely reflects their true adaptation within the tick microbial community.

Conclusion

Here we reported the identification of the *Ixodes ricinus*
microbiota in nymphs collected monthly in three consecutive
years. These results allowed us to show that (1) the
*Ixodes ricinus* microbiota is not stable over time but displays a recurrent temporal pattern that is mainly explained
by the dynamics of environmental taxa; (2) the presence of
TBPs is likely to disturb tick microbial community structure
and hence tick/microbe interactions; (3) some specific
symbionts and “environmental” bacteria may play a key role
in the presence and the dynamics of *I. ricinus*-borne pathogens
and in the defense against parasitoid species. While
the microbial correlations identified in this ecosystem study
need to be confirmed in the near future using experimental
approaches, our new findings suggest that a large part of
the tick microbiome, including environmental taxa, could
play a role in the infectious risk associated with *Ixodes ricinus*, either through bacteria-bacteria interactions or interactions
with the tick and its natural enemies. The tick
microbiome diversity would thus be considered as a promising
resource for the development of new controls strategies
against tick and tick-borne diseases.

Abbreviations

*I. ricinus*: *Ixodes ricinus*; PCA: Principal component analysis; OTU: Operational taxonomic unit; TBP: Tick-borne pathogens

**Supplementary Information**

The online version contains supplementary material available at https://doi.org/10.1186/s40168-021-01051-8.

**Additional file 1.** Distinction of contaminants OTUs, removed from the dataset, from genuine members of tick microbiota.

**Additional file 2** Graphic representation of samples presenting a significant number of sequences belonging to *Arsenophonus* and *Wolbachia* genera. This figure is adapted from the principal component analysis performed on the whole dataset, presented according to axes 1 (24.04%) and 2 (11.16%). (A) Sample projection of the PCA. Samples are colored according to the month of tick sampling. Plotted samples are named as following: ID_Month.Year. Samples presented in larger format and darker colour are those that have been removed from the dataset for the subsequent analyses. (B) Correlation circle of the PCA. OTUs are colored by taxonomic order. OTUs presented in larger format and darker colour are those belonging to *Arsenophonus* and *Wolbachia* genera.

**Additional file 3.** OTUs clusters identified from the total network. OTUs presenting similar connection profiles were identified and assigned to different clusters according to a stochastic block model obtained from the corresponding binary adjacency matrix, obtained from the network inferred by PLN network, which places an edge for each non-null partial correlation (whether positive or negative).

**Additional file 4** Microbiota structure comparison between *Rickettsia* positive and negative samples. Differential analysis comparing OTUs abundances between *Rickettsia* positive and negative samples. Differentially abundant OTUs were defined as those with a p-values < 0.05 after adjustment for multiple testing using the Bonferroni procedure.

**Additional file 5** Microbiota structure comparison between *Borrelia* positive and negative samples. Differential analysis comparing OTUs abundances between *Borrelia* positive and negative samples. Differentially abundant OTUs were defined as those with a p-values < 0.05 after adjustment for multiple testing using the Bonferroni procedure.

**Additional file 6** Microbiota structure comparison between *Anaplasma* positive and negative samples. Differential analysis comparing OTUs abundances between *Anaplasma* positive and negative samples. Differentially abundant OTUs were defined as those with a p-values < 0.05 after adjustment for multiple testing using the Bonferroni procedure.

**Additional file 7** TBPs positive and negative network comparison. The different families of network generated in this study (negative, *Rickettsia*, *Borrelia*, *Anaplasma* and total corrected for TBPs effect) were compared by calculating a weighted version of the Kendall's τ, integrating the edge appearance rank within families of networks (negative, *Rickettsia*, *Borrelia*, *Anaplasma* and total corrected for TBPs effects) as well as associated Pvalues.

Acknowledgements

This work was supported by the métaprogramme "Metaomics and microbial
ecosystems" (MEM) and the Métaprogramme "Adaptation of Agriculture and
Forests to Climate Change" (ACCAF) granted by the French National Institute
for Agriculture, Alimentation and Environment Research (France). The authors
would like to thank Sabine Delannoy and the IdentyPath genomic platform
that allowed us to perform part of the experiments in their laboratory. We also
thank Maxime Galan for providing index sequences and stimulating discussions.
We are grateful to the INRAE MIGALE bioinformatics facility (MIGALE, INRAE,
2020. Migale bioinformatics Facility, doi: 10.15454/1.5572390655343293E12) for
providing help and computing and storage resources. Finally, the authors
would like to thank Olivier Plantard for enriching discussions particularly on the
links between tick microbiota and energetic reserves.

Declarations

Authors' contributions

Conceived and designed the experiments: TP, MVT, JFC, KCM, EL. Performed
the experiments: EL. Analysed the data: EL, AE, JC, AL, SR, TP, MM, CM, OR.
Wrote the paper: EL, AE, JC, AL, SR, OR, MM, CM, KCM, XB, AC, PG, JFC, MVT,
TP. All the authors read and approved the final manuscript.
## Funding

This work was supported by the métaprogramme "Metaomics and microbial ecosystems" (MEM) and the métaprogramme "Adaptation of Agriculture and Forests to Climate Change" (ACCAF) granted by the French National Institute for Agricultural Research (France). The salary of Emilie Lejal, the PhD student working on this project was funded by the Ile-de-France region. Julie Aubert, Julien Chiquet and Stéphane Robin were partially supported by the French ANR grant ANR-17-CE32-0011 Next Generation Biomonitoring (NGB).

## Availability of data and materials

The datasets used for this study can be found in the European Nucleotide Archive. Project accession number: PRJEB36903 (ERP120162) - Sample accession numbers: ERS4353953-ERS4354625 - Run accession numbers: ERR3956669-ERR3957340.

## Ethics approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Competing interests

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Author details

¹UMR BIPAR, Animal Health Laboratory, INRAE, ANSES, Ecole Nationale Vétérinaire d'Alfort, Université Paris-Est, Maisons-Alfort, France. ²Université Paris-Saclay, AgroParisTech, INRAE, UMR MIA-Paris, 75005 Paris, France.
³Faculty of Veterinary Medicine, University of Zaragoza, Zaragoza, Spain.
⁴INRAE, MalAGE, Université Paris-Saclay, Jouy-en-Josas, France. ⁵INRAE, BioinfOmics, MIGALE Bioinformatics Facility, Université Paris-Saclay, Jouy-en-Josas, France.
⁶INRAE, PROSE, Université Paris-Saclay, Antony, France.
⁷Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA, 63122 Saint Genes Champanelle, France. ⁸Université de Lyon, INRAE, VetAgro Sup, UMR EPIA, 69280 Marcy l'Etoile, France. ⁹INRAE, Animal Health Department, Nouzilly, France. ¹⁰UMR ASTRE, CIRAD, INRAE, Campus de Baillarguet, Montpellier, France.

Received: 12 November 2020 Accepted: 17 March 2021

Published online: 03 July 2021

## References

1. Duron O, Noël V, McCoy KD, Bonazzi M, Sidi-Boumedine K, Morel O, et al. The recent evolution of a maternally-inherited endosymbiont of ticks led to the emergence of the Q fever pathogen, Coxiella burnetii. PLoS Pathog. 2015;11(5):e1004892. https://doi.org/10.1371/journal.ppat.1004892.

2. Duron O, Morel O, Noël V, Buysse M, Binetruy F, Lancelot R, et al. Tick-bacteria mutualism depends on B vitamin synthesis pathways. Curr Biol. 2018;28:1896-1902.e5.

3. Vayssier-Taussat M, Kazimirova M, Hubalek Z, Hornok S, Farkas R, Cosson J-F, et al. Emerging horizons for tick-borne pathogens: from the 'one pathogen-one disease' vision to the pathobiome paradigm. Future Microbiol. 2015;10(12):2033-43. https://doi.org/10.2217/fmb.15.114.

4. Bonnet SI, Binetruy F, Hernández-Jarguín AM, Duron O. The tick microbiome: why non-pathogenic microorganisms matter in tick biology and pathogen transmission. Front Cell Infect Microbiol. 2017;7 [cited 2019 Aug 30] Available from: https://www.frontiersin.org/articles/10.3389/fcimb.2017.00236/full.

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

6. Nakao R, Abe T, Nijhof AM, Yamamoto S, Jongejan F, Ikemura T, et al. A novel approach, based on BLSOMs (Batch Learning Self-Organizing Maps), to the microbiome analysis of ticks. ISME J. 2013;7(5):1003-15. https://doi.org/10.1038/ismej.2012.171.

7. Zhang X-C, Yang Z-N, Lu B, Ma X-F, Zhang C-X, Xu H-J. The composition and transmission of microbiome in hard tick, Ixodes persulcatus, during blood meal. Ticks Tick Borne Dis. 2014;5(6):864-70. https://doi.org/10.1016/j.ttbdis.2014.07.009.

8. Van Treuren W, Ponnusamy L, Brinkerhoff RJ, Gonzalez A, Parobek CM, Juliano JJ, et al. Variation in the microbiota of Ixodes ticks with regard to geography, species, and sex. Appl Environ Microbiol. 2015;81(18):6200-9. https://doi.org/10.1128/AEM.01562-15.

9. Gall CA, Reif KE, Scoles GA, Mason KL, Mousel M, Noh SM, et al. The bacterial microbiome of Dermacentor andersoni ticks influences pathogen susceptibility. ISME J. 2016;10(8):1846-55. https://doi.org/10.1038/ismej.2015.266.

10. René-Martellet M, Minard G, Massot R, Tran Van V, Valiente Moro C, Chabanne L, et al. Bacterial microbiota associated with Rhipicephalus sanguineus (s.l.) ticks from France, Senegal and Arizona. Parasit Vectors. 2017;10:416.

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

12. Estrada-Pena A, Cabezas-Cruz A, Pollet T, Vayssier-Taussat M, Cosson J-F. High throughput sequencing and network analysis disentangle the microbial communities of ticks and hosts within and between ecosystems. Front Cell Infect Microbiol. 2018;8:236. https://doi.org/10.3389/fcimb.2018.00236.

13. Varela-Stokes AS, Park SH, Stokes JV, Gavron NA, Lee SI, Moraru GM, et al. Tick microbial communities within enriched extracts of Amblyomma maculatum. Ticks Tick Borne Dis. 2018;9(4):798-805. https://doi.org/10.1016/j.ttbdis.2018.02.022.

14. Tokarz R, Tagliafierro T, Sameroff S, Cucura DM, Oleynik A, Che X, et al. Microbiome analysis of Ixodes scapularis ticks from New York and Connecticut. Ticks Tick Borne Dis. 2019;10(4):894-900. https://doi.org/10.1016/j.ttbdis.2019.04.011.

15. Cabezas-Cruz A, Pollet T, Estrada-Pena A, Allain E, Bonnet SI, Moutailler S. Handling the microbial complexity associated to ticks. Ticks Tick Borne Pathogen. 2018; [cited 2019 Sep 1]; Available from: https://www.intechopen.com/books/ticks-and-tick-borne-pathogens/handling-the-microbial-complexity-associated-to-ticks.

16. Aivelo T, Norberg A, Tschirren B. Bacterial microbiota composition of Ixodes ricinus ticks: the role of environmental variation, tick characteristics and microbial interactions. PeerJ. 2019;7:e8217.

17. Carpi G, Gagnacci F, Wittekindt NE, Zhao F, Qi J, Tomsho LP, et al. Metagenomic profile of the bacterial communities associated with Ixodes ricinus Ticks. PLoS ONE. 2011;6(10):e25604. https://doi.org/10.1371/journal.pone.0025604.

18. Gofton AW, Oskam CL, Lo N, Beninati T, Wei H, McCarl V, et al. Inhibition of the endosymbiont "Candidatus Midichloria mitochondrii" during 16S rRNA gene profiling reveals potential pathogens in Ixodes ticks from Australia. Parasit Vectors. 2015;8(1):345. https://doi.org/10.1186/s13071-015-0958-3.

19. Hernández-Jarguín A, Díaz-Sánchez S, Villar M, de la Fuente J. Integrated metatranscriptomics and metaproteomics for the characterization of bacterial microbiota in unfed Ixodes ricinus. Ticks Tick Borne Dis. 2018;9(5):1241-51. https://doi.org/10.1016/j.ttbdis.2018.04.020.

20. Portillo A, Palomar AM, de Toro M, Santibáñez S, Santibáñez P, Oteo JA. Exploring the bacteriome in anthropophilic ticks: To investigate the vectors for diagnosis. PLoS ONE. 2019;14:e0213384.

21. Guizzo MG, Neupane S, Cucera M, Perner J, Frantová H, da Silva Vaz IJ, et al. Poor unstable midgut microbiome of hard ticks contrasts with abundant and stable monospecific microbiome in ovaries. Front Cell Infect Microbiol, Frontiers. 2020;10 [cited 2020 Sep 21] Available from: https://www.frontiersin.org/articles/10.3389/fcimb.2020.00211/full?report=reader.

22. Lejal E, Estrada-Pena A, Marsot M, Cosson J-F, Rué O, Mariadassou M, et al. Taxon appearance from extraction and amplification steps demonstrates the value of multiple controls in tick microbiota analysis. Front Microbiol, Frontiers. 2020;11 [cited 2020 Sep 14] Available from: https://www.frontiersin.org/articles/10.3389/fmicb.2020.01093/full.

23. Gassner F, van Vliet AJH, Burgers SLGE, Jacobs F, Verbaarschot P, Hovius EKE, et al. Geographic and temporal variations in population dynamics of Ixodes ricinus and associated Borrelia infections in The Netherlands. Vector Borne Zoonotic Dis. 2010;11:523-32.

24. Coipan EC, Jahfari S, Fonville M, Maassen C, van der Giessen J, Takken W, et al. Spatiotemporal dynamics of emerging pathogens in questing Ixodes ricinus. Front Cell Infect Microbiol. 2013;3:36.

25. Takken W, van Vliet AJH, Verhulst NO, Jacobs FHH, Gassner F, Hartemink N, et al. Acarological risk of Borrelia burgdorferi Sensu Lato infections across blood meal. Ticks Tick Borne Dis. 2014;5(6):864-70. https://doi.org/10.1016/j.ttbdis.2014.07.009.
space and time in The Netherlands. Vector Borne Zoonotic Dis. 2016;17:99–107.

26. Lejal E, Marsot M, Chalvet-Monfray K, Cosson J-F, Moutailler S, Vayssier-Taussat M, et al. A three-years assessment of Ixodes ricinus-borne pathogens in a French peri-urban forest. Parasit Vectors. 2019;12(1):551. https://doi.org/10.1186/s13071-019-3799-7.

27. Pollet T, Sprong H, Lejal E, Krawczyk AI, Moutailler S, Cosson J-F, et al. The scale affects our view on the identification and distribution of microbial communities in ticks. Parasit Vectors. 2020;13(1):36. https://doi.org/10.1186/s13071-020-3908-7.

28. Budachetri K, Kumar D, Crispell G, Beck C, Dasch G, Karim S. The tick endosymbiont Candidatus Midichloria mitochondrii and selenoproteins are essential for the growth of Rickettsia parkeri in the Gulf Coast tick vector. Microbiome. 2018;6(1):141. https://doi.org/10.1186/s40168-018-0524-2.

29. Lejal E, Moutailler S, Šimo L, Vayssier-Taussat M, Pollet T. Tick-borne pathogen detection in midgut and salivary glands of adult Ixodes ricinus. Parasit Vectors. 2019;12(1):152. https://doi.org/10.1186/s13071-019-3418-7.

30. Galan M, Razzauti M, Bard E, Bernard M, Brouat C, Charbonnel N, et al. 16S rRNA amplicon sequencing for epidemiological surveys of bacteria in wildlife. mSystems. 2016;1:e00032–16.

31. Aitchison J, Ho CH. The multivariate Poisson-log normal distribution. Biometrika. 1989;76(4):643–53. https://doi.org/10.1093/biomet/76.4.643.

32. Chiquet J, Mariadassou M, Robin S. PLNmodels: Poisson lognormal models. 2019. [cited 2019 Oct 19]. Available from: https://CRAN.R-project.org/packagge=PLNmodels

33. Chiquet J, Mariadassou M, Robin S. Variational inference for sparse network reconstruction from count data. arXiv: 2018:180603120 [stat] [cited 2019 Oct 18]; Available from: http://arxiv.org/abs/1806.03120.

34. Chiquet J, Mariadassou M, Robin S. Variational inference for probabilistic Poisson PCA. Ann Appl Stat. 2018;12:2674–98.

35. Leger J-B. Blockmodels: A R-package for estimating in Latent Block Model and Stochastic Block Model, with various probability functions, with or without covariates. arXiv: 2016:160207587 [stat] [cited 2020 Nov 9]; Available from: http://arxiv.org/abs/201602.07587.

36. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

37. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25. https://doi.org/10.1186/gb-2010-11-3-r25.

38. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40(10):4288–97. https://doi.org/10.1093/nar/gks042.

39. R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2020. Available from: https://www.R-project.org/

40. Heise SR, Elshahed MS, Little SE. Bacterial diversity in Amblyomma americanum (Acari: Ixodidae) with a focus on members of the Genus Rickettsia. J Med Entomol. 2010;47(2):258–68. https://doi.org/10.1093/jmedent/47.2.258.

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

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

43. Trout Fryxell RT, DeBruyn JM. The microbiome of Ehrlichia-infected and uninfected lone star ticks (Amblyomma americanum). Stevenson B, editor. PLoS ONE. 2016;11:e0146651.

44. Gurfield N, Grewal S, Cua LS, Torres PJ, Kelley ST. Endosymbiont interference and microbial diversity of the Pacific coast tick, Dermacentor occidentalis, in San Diego County, California. PeerJ. 2017;5:e3202. https://doi.org/10.7717/peerj.3202.

45. Thapa S, Zhang Y, Allen MS. Effects of temperature on bacterial microbiome composition in Ixodes scapularis ticks. Microbiol Open. 2019;8(5):e00719. https://doi.org/10.1002/mbo.719.

46. Duron O, Binetruy F, Noël V, Cremaschi J, McCoy KD, Arnathau C, et al. Evolutionary changes in symbiont community structure in ticks. Mol Ecol. 2017;26(11):2905–21. https://doi.org/10.1111/mec.14094.

47. Beninati T, Lo N, Sacchi L, Genchi C, Noda H, Bandi C. A novel alpha-proteobacterium resides in the mitochondria of ovarian cells of the tick Ixodes ricinus. Appl Environ Microbiol. 2004;70(5):2596–602. https://doi.org/10.1128/AEM.70.5.2596-2602.2004.

48. Lo N, Beninati T, Sassera D, Bouman EA, Santagati S, Gern L, et al. Widespread distribution and high prevalence of an alpha-proteobacterial symbiont in the tick Ixodes ricinus. Environ Microbiol. 2006;8:1280–7.

49. Sassera D, Beninati T, Bandi C, Bouman EAP, Sacchi L, Fabbi M, et al. 'Candidatus Midichloria mitochondrii', an endosymbiont of the tick Ixodes ricinus with a unique intramitochondrial lifestyle. Int J Syst Evol Microbiol. 2006;56(11):2535–40. https://doi.org/10.1099/ijsem.0.64386-0.

50. Olivieri E, Epis S, Castelli M, Varotto Boccazzi I, Romeo C, Desirò A, et al. Tissue tropism and metabolic pathways of Midichloria mitochondrii suggest tissue-specific functions in the symbiosis with Ixodes ricinus. Ticks Tick Borne Dis. 2019;10(5):1070–7. https://doi.org/10.1016/j.ttbdis.2019.05.019.

51. Schicht S, Schnieder T, Strube C. Rickettsia spp. and coinfections with other pathogenic microorganisms in hard ticks from Northern Germany. J Med Entomol. 2012;49(3):766–71. https://doi.org/10.1603/ME11204.

52. Raulf M-K, Jordan D, Fingerle V, Strube C. Association of Borrelia and Rickettsia spp. and bacterial loads in Ixodes ricinus ticks. Ticks Tick Borne Dis. 2018;9(1):18–24. https://doi.org/10.1016/j.ttbdis.2017.10.014.

53. Plantard O, Bouju-Albert A, Malard M-A, Hermouet A, Capron G, Verheyden H. Detection of Wolbachia in the Tick Ixodes ricinus is due to the presence of the Hymenoptera Endoparasitoid Ixodiphagus hookeri. PLoS ONE. 2012; 7(1):e30692. https://doi.org/10.1371/journal.pone.0030692.

54. Bohacsova M, Mediannikov O, Kazimirova M, Raoult D, Sekeyova Z. Arsenophonus nasoniae and Rickettsiae infection of Ixodes ricinus due to parasitic wasp Ixodiphagus hookeri. PLoS ONE. 2016;11(2):e0149950. https://doi.org/10.1371/journal.pone.0149950.

55. Fytrou A, Schofield PG, Kraaijieveld AR, Hubbard SF. Wolbachia infection suppresses both host defence and parasitoid counter-defence. Proc R Soc B Biol Sci. 2006;273:791–6.

56. Binetruy F, Bailly X, Chevillon C, Martin OY, Bernasconi MV, Duron O. Phylogenetics of the Spiroplasma ixodetis endosymbiont reveals past transfers between ticks and other arthropods. Ticks Tick Borne Dis. 2019; 10(3):575–84. https://doi.org/10.1016/j.ttbdis.2019.02.001.

57. Oliver KM, Russell JA, Moran NA, Hunter MS. Facultative bacterial symbionts in aphids confer resistance to parasitic wasps. PNAS. 2003;100(4):1803–7. https://doi.org/10.1073/pnas.0335320100.

58. Guay J-F, Boudreault S, Michaud D, Cloutier C. Impact of environmental stress on aphid clonal resistance to parasitoids: Role of Hamiltonella defensa bacterial symbiosis in association with a new facultative symbiont of the pea aphid. J Insect Physiol. 2009;55(10):919–26. https://doi.org/10.1016/j.jinsphys.2009.06.006.

59. Towner K. The Genus Acinetobacter. In: Dworkin M, Falkow S, Rosenberg E, Schleifer KH, Stackerbrandt E, editors. The Prokaryotes: Volume 6: Proteobacteria: Gamma Subclass. New York: Springer New York; 2006. p. 746–58 [cited 2019 Oct 11]. Available from: https://doi.org/10.1007/0-387-30746-X_25.

60. Falkinham JO. Environmental sources of nontuberculous Mycobacteria. Clin Chest Med. 2015;36(1):35–41. https://doi.org/10.1016/j.ccm.2014.10.003.

61. Kämpfer P, McInroy JA, Glaeser SP. Chryseobacterium rhizoplanae sp. nov., isolated from the rhizoplane environment. Antonie Van Leeuwenhoek. 2015; 107:533–8.

62. Horn H, Keller A, Hildebrandt U, Kämpfer P, Riederer M, Hentschel U. Draft genome of the Arabidopsis thaliana phyllosphere bacterium, Williamsia sp. ARP1. Stand Genomic Sci. 2016;11:8.

63. An S, Berg G. Stenotrophomonas maltophilia. Trends Microbiol. 2018;26(7): 637–8. https://doi.org/10.1016/j.tim.2018.04.006.

64. Green PN, Ardley JK. Review of the genus Methylobacterium and closely related organisms: a proposal that some Methylobacterium species be reclassified into a new genus, Methylorubrum gen. nov. Int J Syst Evol Microbiol. 2018;68(9):2727–48. https://doi.org/10.1099/ijsem.0.002856.

65. Salter SJ, Cox MJ, Turek EM, Calus ST, Cookson WO, Moffatt MF, et al. Reagent and laboratory contamination can critically impact sequence-based microbiome analyses. BMC Biol. 2014;12(1):87. https://doi.org/10.1186/s12915-014-0087-z.

66. Glassing A, Dowd SE, Galandiuk S, Davis B, Chiodini RJ. Inherent bacterial DNA contamination of extraction and sequencing reagents may affect interpretation of microbiota in low bacterial biomass samples. Gut Pathogens. 2016;8(1):24. https://doi.org/10.1186/s13099-016-0103-7.
67. Weyrich LS, Farrer AG, Eisenhofer R, Arriola LA, Young J, Selway CA, et al. Laboratory contamination over time during low-biomass sample analysis. Mol Ecol Resour. 2019;19(4):982–96. https://doi.org/10.1111/1755-0998.13 011.

68. Binetruy F, Dupraz M, Buyse M, Duron O. Surface sterilization methods impact measures of internal microbial diversity in ticks. Parasit Vectors. 2019; 12(1):268. https://doi.org/10.1186/s13071-019-3517-5.

69. Baldrian P. Microbial activity and the dynamics of ecosystem processes in forest soils. Curr Opin Microbiol. 2017;37:128–34. https://doi.org/10.1016/j. mib.2017.06.008.

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

71. Narasimhan S, Rajeevan N, Liu L, Zhao YO, Heisig J, Pan J, et al. Gut microbiota of the tick vector Ixodes scapularis modulate colonization of the Lyme disease spirochete. Cell Host Microbe. 2014;15(1):58–71. https://doi. org/10.1016/j.chom.2013.12.001.

72. Landesman WJ, Mulder K, Allan BF, Bashor LA, Keesing F, LoGiudice K, et al. Potential effects of blood meal host on bacterial community composition in Ixodes scapularis nymphs. Ticks Tick Borne Dis. 2019;10(3):523–7. https://doi. org/10.1016/j.ttbdis.2019.01.002.

73. Randolph SE, Green RM, Hoodless AN, Peacey MF. An empirical quantitative framework for the seasonal population dynamics of the tick Ixodes ricinus. Int J Parasitol. 2002;32(8):979–89. https://doi.org/10.1016/S0020-7519(02) 00030-9.

74. Abdullah S, Davies S, Wall R. Spectrophotometric analysis of lipid used to examine the phenology of the tick Ixodes ricinus. Parasit Vectors. 2018;11(1): 523. https://doi.org/10.1186/s13071-018-3102-3.

75. Eiler A, Heinrich F, Bertilsson S. Coherent dynamics and association networks among lake bacterioplankton taxa. ISME J. 2012;6(2):330–42. https://doi.org/10.1038/ismej.2011.113.

76. Chow C-ET, Kim DY, Sachdeva R, Caron DA, Fuhrman JA. Top-down controls on bacterial community structure: microbial network analysis of bacteria, T4-like viruses and protists. ISME J. 2014;8(4):816–29. https://doi.org/10.1038/ ismej.2013.199.

77. Fuhrman JA, Steele JA. Community structure of marine bacterioplankton: patterns, networks, and relationships to function. Aquat Microb Ecol. 2008; 53:69–81. https://doi.org/10.3354/ame01222.

78. Berdjeb L, Parada A, Needham DM, Fuhrman JA. Short-term dynamics and interactions of marine protist communities during the spring–summer transition. ISME J. 2018;12(8):1907–17. https://doi.org/10.1038/s41396-018- 0097-x.

79. Abraham NM, Liu L, Jutras BL, Yadav AK, Narasimhan S, Gopalakrishnan V, et al. Pathogen-mediated manipulation of arthropod microbiota to promote infection. PNAS. 2017;114(5):E781–90. https://doi.org/10.1073/pnas.1 613422114.

80. Krawczyk AI, Bakker JW, Koenraadt CJM, Fonville M, Takumi K, Sprong H, et al. Tripartite interactions among Ixodiphagus hookeri, Ixodes ricinus and deer: differential interference with transmission cycles of tick-borne pathogens. Pathogens. 2020;9:339.

81. Tully JG, Rose DL, Yunker CE, Carle P, Bové JM, Williamson DL, et al. Spiroplasma ixodetis sp. nov., a new species from Ixodes pacificus ticks collected in Oregon. Int J Syst Evol Microbiol. 1995;45:23–8.

82. Henning K, Greiner-Fischer S, Hotzel H, Ebsen M, Theegarten D. Isolation of Spiroplasma sp. from an Ixodes tick. Int J Med Microbiol. 2006;296:157–61.

83. Weinert LA, Tinsley MC, Temperley M, Jiggins FM. Are we underestimating the diversity and incidence of insect bacterial symbionts? A case study in ladybird beetles. Biol Lett. 2007;3(6):678–81. https://doi.org/10.1098/rsbl. 2 007.0373.

84. Duron O, Bouchon D, Boutin S, Bellamy L, Zhou L, Engelstädter J, et al. The diversity of reproductive parasites among arthropods: Wolbachia do not walk alone. BMC Biol. 2008;6(1):27. https://doi.org/10.1186/1741-7007-6-27.

85. Hurst GDD, von der Schulenburg JHG, Majerus TMO, Bertrand D, Zakharov IA, Baungarda J, et al. Invasion of one insect species, Adalia bipunctata, by two different male-killing bacteria. Insect Mol Biol. 1999;8(1):133–9. https:// doi.org/10.1046/j.1365-2583.1999.810133.x.

86. Majerus TMO, Schulenburg JHGVD, Majerus MEN, Hurst GDD. Molecular identification of a male-killing agent in the ladybird Harmonia axyridis (Pallas) (Coleoptera: Coccinellidae). Insect Mol Biol. 1999;8(4):551–5. https:// doi.org/10.1046/j.1365-2583.1999.00151.x.

87. Majerus MEN, Hinrich J, Schulenburg GVD, Zakharov IA. Multiple causes of male-killing in a single sample of the two-spot ladybird, Adalia bipunctata

(Coleoptera: Coccinellidae) from Moscow. Heredity. 2000;84(5):605–9. https://doi.org/10.1046/j.1365-2540.2000.00710.x.

88. Jiggins FM, Hurst GDD, Jiggins CD, v d Schulenburg JH, Majerus MEN. The butterfly Danaus chrysippus is infected by a male-killing Spiroplasma bacterium. Parasitology. 2000;120:439–46.

89. Tinsley MC, Majerus MEN. A new male-killing parasitism: Spiroplasma bacteria infect the ladybird beetle Anisosticta novemdecimpunctata (Coleoptera: Coccinellidae). Parasitology. 2006;132(6):757–65. https://doi. org/10.1017/S0031182005007989.

90. Sanada-Morimura S, Matsumura M, Noda H. Male killing caused by a Spiroplasma symbiont in the small brown planthopper, Laodelphax striatellus. J Hered. 2013;104(6):821–9. https://doi.org/10.1093/jhered/est052.

91. Xie J, Butler S, Sanchez G, Mateos M. Male killing Spiroplasma protects Drosophila melanogaster against two parasitoid wasps. Heredity. 2014; 112(4):399–408. https://doi.org/10.1038/hdy.2013.118.

92. Peix A, Ramírez-Bahena M-H, Velázquez E. Historical evolution and current status of the taxonomy of genus Pseudomonas. Infect Genet Evol. 2009;9(6): 1132–47. https://doi.org/10.1016/j.meegid.2009.08.001.

93. Meikle WG, Mercadier G, Guermache F, Bon M-C. Pseudomonas contamination of a fungus-based biopesticide: implications for honey bee (Hymenoptera: Apidae) health and Varroa mite (Acari: Varroidae) control. Biol Control. 2012;60(3):312–20. https://doi.org/10.1016/j.biocontrol.2011.12. 004.

94. Meikle WG, Bon M-C, Cook SC, Gracia C, Jaronski ST. Two strains of Pseudomonas fluorescens bacteria differentially affect survivorship of waxworm (Galleria mellonella) larvae exposed to an arthropod fungal pathogen, Beauveria bassiana. Biocontrol Sci Tech. 2013;23:220–33.

95. Vodovar N, Vinals M, Liehl P, Basset A, Degrouard J, Spellman P, et al. Drosophila host defense after oral infection by an entomopathogenic Pseudomonas species. PNAS. 2005;102(32):11414–9. https://doi.org/10.1073/ pnas.0502240102.

96. Vodovar N, Vallenet D, Cruveiller S, Rouy Z, Barbe V, Acosta C, et al. Complete genome sequence of the entomopathogenic and metabolically versatile soil bacterium Pseudomonas entomophila. Nat Biotechnol. 2006; 24(6):673–9. https://doi.org/10.1038/nbt1212.

97. Evtushenko LI, Takeuchi M. The Family Microbacteriaceae. In: Dworkin M, Falkow S, Rosenberg E, Schleifer KH, Stackerbrandt E, editors. The Prokaryotes: Volume 3: Archaea Bacteria: Firmicutes, Actinomycetes. New York: Springer New York; 2006. p. 1020–98. [cited 2019 Oct 19]. Available from: https://doi.org/10.1007/0-387-30743-5_43.

98. Jiang Y, Wiese J, Tang SK, Xu LH, Imhoff JF, Jiang CL. Actinomycetospora chiangmaiensis gen. nov., sp. nov., a new member of the family Pseudonocardiaceae. Int J Syst Evol Microbiol. 2008;58:408–13.

99. Lee SD. Kineococcus rhizosphaerae sp. nov., isolated from rhizosphere soil. Int J Syst Evol Microbiol. 2009;59(9):2204–7. https://doi.org/10.1099/ijs.0. 008599-0.

100. Yamamura H, Ashizawa H, Nakagawa Y, Hamada M, Ishida Y, Otoguro M, et al. Actinomycetospora rishiriensis sp. nov., isolated from a lichen. Int J Syst Evol Microbiol. 2011;61(11):2621–5. https://doi.org/10.1099/ijs.0.028753- 0.

101. Eberl L, Vandamme P. Members of the genus Burkholderia: good and bad guys. F1000Res. 2016:5 [cited 2019 Oct 19]; Available from: https://www. ncbi.nlm.nih.gov/pmc/articles/PMC4882756/.

102. Viana AT, Caetano T, Covas C, Santos T, Mendo S. Environmental superbugs: the case study of Pedobacter spp. Environ Pollut. 2018;241:1048–55. https:// doi.org/10.1016/j.envpol.2018.06.047.

103. Burgdorfer W, Hayes SF, Mavros AJ. Nonpathogenic rickettsiae in Dermacentor andersoni: a limiting factor for the distribution of Rickettsia rickettsii. Rickettsiae Rickettsial Dis. 1981;58:585–94.

104. Macaluso KR, Sonenshine DE, Ceraul SM, Azad AF. Rickettsial infection in Dermacentor variabilis (Acari: Ixodidae) inhibits transovarial transmission of a second Rickettsia. J Med Entomol. 2002;39(6):809–13. https://doi.org/10.1 603/0022-2585-39.6.809.

105. Steiner FE, Pinger RR, Vann CN, Grindle N, Civitello D, Clay K, et al. Infection and co-infection rates of Anaplasma phagocytophilum variants, Babesia spp. , Borrelia burgdorferi, and the Rickettsial Endosymbiont in Ixodes scapularis (Acari: Ixodidae) from Sites in Indiana, Maine, Pennsylvania, and Wisconsin. J Med Entomol. 2008;45(2):289–97. https://doi.org/10.1093/jmedent/45.2.289.

106. Telford SR. Status of the "East side hypothesis" (transovarial interference) twenty five years later. Ann N Y Acad Sci. 2009;1166(1):144–50. https://doi. org/10.1111/j.1749-6632.2009.04522.x.
107. Williams-Newkirk AJ, Rowe LA, Mixson-Hayden TR, Dasch GA. Characterization of the bacterial communities of life stages of free living lone star ticks (Amblyomma americanum). Munderloh UG, editor. PLoS ONE. 2014;9:e102130.

108. Tijsse-Klasen E, Sprong H, Pandak N. Co-infection of Borrelia burgdorferi sensu lato and Rickettsia species in ticks and in an erythema migrans patient. Parasit Vectors. 2013;6(1):347. https://doi.org/10.1186/1756-3305-6-347.

109. Moniuszko A, Dunaj J, Święcicka I, Zambrowski G, Chmielewska-Badora J, Żukiewicz-Sobczak W, et al. Co-infections with Borrelia species, Anaplasma phagocytophilum and Babesia spp. in patients with tick-borne encephalitis. Eur J Clin Microbiol Infect Dis. 2014;33(10):1835-41. https://doi.org/10.1007/s10096-014-2134-7.

110. Hoversten K, Bartlett MA. Diagnosis of a tick-borne coinfection in a patient with persistent symptoms following treatment for Lyme disease. BMJ Case Rep. 2018; [cited 2018 Nov 19]; Available from: http://europepmc.org/abstracct/med/30262525.

111. Krause PJ, Telford SR, Spielman A, Sikand V, Ryan R, Christianson D, et al. Concurrent Lyme disease and Babesiosis: evidence for increased severity and duration of illness. JAMA. 1996;275(21):1657-60. https://doi.org/10.1001/jama.1996.03530450047031.

112. Diuk-Wasser MA, Vannier E, Krause PJ. Coinfection by Ixodes tick-borne pathogens: ecological, epidemiological, and clinical consequences. Trends Parasitol. 2016;32(1):30-42. https://doi.org/10.1016/j.pt.2015.09.008.

113. Chauhan G, McClure J, Hekman J, Marsh PW, Bailey JA, Daniels RF, et al. Combining citizen science and genomics to investigate tick, pathogen, and commensal microbiome at single-tick resolution. Front Genet. 2020:10 [cited 2021 Jan 14] Available from: https://www.frontiersin.org/articles/10.3389/fgene.2019.01322/full?utm_source=S-TWT&utm_medium=SNET&utm_campaign=ECO_FGENE_XXXXXXX_auto-dlvrit.

114. Brinkerhoff RJ, Clark C, Ocasio K, Gauthier DT, Hynes WL. Factors affecting the microbiome of Ixodes scapularis and Amblyomma americanum. PLoS ONE. 2020;15:e0232398.

115. Adegoke A, Kumar D, Bobo C, Rashid MI, Durrani AZ, Sajid MS, et al. Tick-borne pathogens shape the native microbiome within tick vectors. Microorganisms. 2020;8:1299.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions