# Rickettsial pathogens drive microbiota assembly in Hyalomma marginatum and Rhipicephalus bursa ticks

Apolline Maitre, Alejandra Wu-Chuang, Lourdes Mateos-Hernández, Elianne

Piloto-Sardiñas, Angélique Foucault-Simonin, Vincent Cicculli, Sara

Moutailler, Jean-christophe Paoli, Alessandra Falchi, Dasiel Obregón, et al.

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

► To cite this version:

Apolline Maitre, Alejandra Wu-Chuang, Lourdes Mateos-Hernández, Elianne Piloto-Sardiñas, Angélique Foucault-Simonin, et al. Rickettsial pathogens drive microbiota assembly in Hyalomma marginatum and Rhipicephalus bursa ticks. Molecular Ecology, 2023, 10.1111/mec.17058. hal-04150492

HAL Id: hal-04150492

https://hal.science/hal-04150492

Submitted on 7 Jul 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Rickettsial pathogens drive microbiota assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks

Apolline Maitre<fnref n="1" />,<fnref n="2" />,<fnref n="3" />

Alejandra Wu-Chuang<fnref n="1" />

Lourdes Mateos-Hernández<fnref n="1" />

Elianne Piloto-Sardiñas<fnref n="4" />

Angélique Foucault-Simonin<fnref n="1" />

Vincent Cicculli<fnref n="3" />

Sara Moutailler<fnref n="1" />

Jean-Christophe Paoli<fnref n="2" />

Alessandra Falchi<fnref n="3" />

Dasiel Obregón<fnref n="5" /> ID

Alejandro Cabezas-Cruz<fnref n="1" /> ID

**Correspondence**

Alejandro Cabezas-Cruz and Apolline Maitre, ANSES, INRAE, Ecole Nationale Vétérinaire d'Alfort, UMR BIPAR, Laboratoire de Santé Animale, Maisons-Alfort F-94700, France.

Email: alejandro.cabezas@vet-alfort.fr and apolline.maitre@anses.fr

**Funding information**

Agence Nationale de la Recherche, Grant/Award Number: ANR-10-LABX-62-IBEID

**Handling Editor:** Sebastien Calvignac-Spencer

**Abstract**

Most tick-borne pathogens (TBPs) are secondarily acquired by ticks during feeding on infected hosts, which imposes ‘priority effect’ constraints, as arrival order influences the establishment of new species in a microbial community. Here we tested whether once acquired, TBPs contribute to bacterial microbiota functioning by increasing community stability. For this, we used *Hyalomma marginatum* and *Rhipicephalus bursa* ticks collected from cattle in different locations of Corsica and combined 16S rRNA ampli-con sequencing and co-occurrence network analysis, with high-throughput pathogen detection, and in silico removal of nodes to test for impact of rickettsial pathogens on network properties. Despite its low centrality, *Rickettsia* showed preferential connections in the networks, notably with a keystone taxon in *H. marginatum*, suggesting facilitation of *Rickettsia* colonisation by the keystone taxon. In addition, conserved patterns of community assembly in both tick species were affected by *Rickettsia* removal, suggesting that privileged connections of *Rickettsia* in the networks make this taxon a driver of community assembly. However, *Rickettsia* removal had minor impact on the conserved ‘core bacterial microbiota’ of *H. marginatum* and *R. bursa*. Interestingly, networks of the two tick species with *Rickettsia* have similar node centrality distribution, a property that is lost after *Rickettsia* removal, suggesting that this taxon drives specific hierarchical interactions between bacterial microbes in the microbiota. The study indicates that tick-borne *Rickettsia* play a significant role in the tick bacterial microbiota, despite their low centrality. These bacteria are influential and contribute to the conservation of the ‘core bacterial microbiota’ while also promoting community stability.

**KEYWORDS**

*Hyalomma marginatum*, microbiota, networks, *Rhipicephalus bursa*, *Rickettsia*, ticks

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

This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2023 The Authors. *Molecular Ecology* published by John Wiley & Sons Ltd.

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

<footnote n="2">²INRAE, UR 0045 Laboratoire de Recherches Sur Le Développement de L'Elevage (SELMET-LRDE), Corte, France</footnote>

<footnote n="3">³EA 7310, Laboratoire de Virologie, Université de Corse, Corte, France</footnote>

<footnote n="4">⁴Direction of Animal Health, National Center for Animal and Plant Health, Carretera de Tapaste y Autopista Nacional, San José de las Lajas, Cuba</footnote>

<footnote n="5">⁵School of Environmental Sciences, University of Guelph, Guelph, Ontario, Canada</footnote>
# 1 | INTRODUCTION

Ticks have intricate ecological interactions with microorganisms present in their microbiota (Wu-Chuang et al., 2021; Wu-Chuang, Obregon, Estrada-Peña, et al., 2022). Some of these microorganisms are considered human pathogens, such as *Rickettsia sibirica mongolii*-timonae (Jakimovski et al., 2022). Ticks also harbour endosymbionts like *Coxiella*-like (Zhong, 2012) or *Francisella*-like endosymbionts (Duron et al., 2018), as well as commensal bacteria (Wu-Chuang et al., 2021; Wu-Chuang, Obregon, Estrada-Peña, et al., 2022). While the distinction between human pathogens and other non-pathogenic microorganisms carried by ticks is still being investigated, it is possible to differentiate between obligate tick endosymbionts, which are vertically transmitted from mother to offspring and mainly reside in vector reproductive organs, and commensal bacteria, which are horizontally acquired and heavily influenced by the ecological conditions in which ticks are found (Hussain et al., 2022; Wu-Chuang et al., 2021; Wu-Chuang, Obregon, Estrada-Peña, et al., 2022). The composition of commensal bacteria in internal (e.g. salivary glands) and external (e.g. cuticle) tick tissues, referred here as bacterial microbiota (Wu-Chuang et al., 2021; Wu-Chuang, Obregon, Estrada-Peña, et al., 2022), is a dynamic entity that can be directly or indirectly influenced by biotic factors like pathogen infections (Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022) or abiotic factors like temperature (Thapa et al., 2019; Wu-Chuang et al., 2021; Wu-Chuang, Obregon, Estrada-Peña, et al., 2022).

Microbiota dynamics are relevant for tick biology and disease ecology. On the one hand, manipulation of tick microbiota by antibiotics affected the proportions of bacterial taxa in the microbiota and significantly reduced the reproductive fitness of *Dermacentor andersoni* (Clayton et al., 2015) and *Amblyomma americanum* (Zhong et al., 2007). Also, *Ixodes ricinus* fed on hosts immunised with anti-microbiota vaccines had an increased engorgement weight and lower bacterial diversity in the bacterial microbiota (Mateos-Hernández et al., 2020, 2021). On the other hand, empirical studies under controlled laboratory conditions have shown that tick bacterial microbiota dysbiosis can cause lower *Borrelia burgdorferi* colonisation in *Ixodes scapularis* (Narasimhan et al., 2014) or higher transstadial transmission of *Babesia microti* in *Haemaphysalis longicornis* (Wei et al., 2021). Furthermore, strong relationships have been found between the presence of specific pathogens and the diversity (Sperling et al., 2020) and assembly (Lejal et al., 2021; Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022) of the bacterial microbiota of questing ticks (Lejal et al., 2021), and ticks feeding on domestic animals (Sperling et al., 2020) or humans (Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022).

Most tick-borne pathogens (TBPs) are secondarily acquired by ticks during blood feeding on an infected host, which poses a challenge for these pathogens in establishing themselves within an already established bacterial community. The order and timing of their arrival influence the establishment of new species, known as

a 'priority effect' (Debray et al., 2022). How do TBPs interact with an already established bacterial community within a tick upon their arrival? What is the significance of TBPs in the process of community assembly? Additionally, it is worth considering how the order of colonisation influences both the outcome of community assembly and the ecological success of individual colonisers, as point out by Martínez et al. (2018) in their study. To address these questions, we turn our attention to emergent properties of bacterial co-occurrence networks (Röttjers & Faust, 2018), and asked whether the nesting of TBPs in the tick bacterial microbiota networks could be such that their removal could impact community assembly and emergent properties. Emergent properties (e.g. robustness and connectivity) explain the behaviour of complex systems such as the bacterial microbiota and these would not be observed if parts of the network are investigated in isolation (Aderem, 2005; Röttjers & Faust, 2018). A previous report showed that despite changes in the relative bacterial microbiota composition of laboratory-reared *I. scapularis* infected with *Anaplasma phagocytophilum* (Abraham et al., 2017; Estrada-Peña et al., 2020), the network robustness decreased only marginally in *A. phagocytophilum*-infected ticks (Estrada-Peña et al., 2020). In contrast to pathogen infection, anti-tick immunity significantly reduced the network robustness (Estrada-Peña et al., 2020). An additional study found that the connectivity and centrality of the taxon *Rickettsia* was significantly reduced in networks of *R. helvetica*-infected ticks (Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022). This suggests that TBP infection induces minimal changes on the tick bacterial microbiota robustness (Estrada-Peña et al., 2020), and that once established; pathogens occupy marginal positions in the bacterial community (Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022), with yet untested consequences for the robustness.

Analysing and comparing tick bacterial microbiota requires investigating the effects of microbial members on their environment and the interactions between bacteria within the microbiota. Diversity metrics alone are insufficient for studying these mechanisms (Röttjers & Faust, 2018). Network analysis provides a powerful tool to explore the complexity of microbial communities, representing interactions between bacterial taxa as nodes and co-occurrence as edges. With network analysis, the abundance of bacteria is not the sole metric for assessing their importance in the community. Network studies can provide information on community hierarchical organisation, architecture, assembly, stability, and interactions, which are crucial for understanding community dynamics (Guseva et al., 2022; Röttjers & Faust, 2018). In addition, networks allow measuring the strength of the bacterial interactions which are fundamental to understand community assembly and stability (Coyte et al., 2021). Examining positive and negative co-occurrence interactions through microbial networks helps analyse the consequences of these interactions on microbial fitness, population dynamics, and metabolic functions (Berg et al., 2020). Network analysis is, therefore, an ideal tool to assess the position, role, and importance of TBPs in community assembly.
In this study, our aim was to investigate the relationship between secondarily acquired rickettsial pathogens and the overall assembly of the tick bacterial microbiota. Specifically, we sought to determine whether the presence of these pathogens in the microbial community could influence the stability of the tick bacterial network. To achieve this, *Rhipicephalus bursa* and *Hyalomma marginatum* ticks were collected from cattle in different locations of Corsica. The research approach involved analysing the bacterial microbiota using 16S rRNA amplicon sequencing, conducting co-occurrence network analysis, and employing high-throughput pathogen detection. To evaluate the impact of rickettsial pathogens on network properties, we performed in silico removal of nodes from the networks. The results of our study revealed that a high occurrence of single rickettsial pathogen infection was associated with lower alpha diversity in the bacterial microbiota of *H. marginatum* compared to *R. bursa*. Removing *Rickettsia* from the networks had a minor impact on the 'core bacterial microbiota' but affected network robustness and connectivity, leading to reshaping of bacterial community assembly in both tick species. We conclude that once acquired by ticks, rickettsial pathogens play a major role in bacterial community assembly in *H. marginatum* and *R. bursa* with minimal impact on the 'core bacterial microbiota'.

# 2 | MATERIALS AND METHODS

## 2.1 | Tick samples

Ticks were collected between 2019 and 2021 in Ponte-Leccia slaughterhouse (42°46'45" N, 9°20'83"E) on Corsican breed cattle. Ticks were collected manually on bovine skins after slaughter. The identification number on the ear loop was collected to know the location of the collected bovines and therefore the ticks. All the bovines sampled were located in Corsica. Ticks were morphologically identified to the species level under a dissecting microscope (Nikon SMZ445) following the dichotomous keys (Estrada-Peña et al., 2004). In total, samples of eight species (i.e. *Dermacentor marginatus* (n=8), *Ixodes ricinus* (n=44), *Haemaphysalis punctata* (n=46), *Hyalomma marginatum* (n=702), *Hyalomma scupense* (n=198), *Rhipicephalus bursa* (n=1473), *Rhipicephalus sanguineus* (n=69) and *Rhipicephalus annulatus* (n=147)) were collected including different sexes (i.e. male and female) and different stages (i.e. nymphs and adults engorged or not). Adult females (engorged or not) from the two most abundant ticks in Corsica (Grech Angelini et al., 2016), *Rhipicephalus bursa* (n=24) and *Hyalomma marginatum* (n=23), were selected randomly for this study. The ticks were collected from different bovines, of which location of origin were listed in Table S1. Before DNA extraction, ticks were washed two times in miliQ sterile water and one time in 70% ethanol. Different surface sterilisation methods (i.e. ethanol vs. bleach) did not significantly affect the bacterial diversity of tick midguts, but the whole-body microbiota of bleach-treated ticks have significantly lower bacterial diversity compared with ethanol-treated

ticks, as bleach removes external microbes (Binetruy et al., 2019). Accordingly, ethanol washing was used here to; explicitly include the internal and external tick microbiota in the analysis, as we consider tick surface microbes as part of the tick microbiota. After washing, ticks were conserved in 70% ethanol and stored at -80°C until further processing.

## 2.2 | DNA extraction and pre-amplification

Ticks were homogenised using TissueLyser II (Qiagen) in a phosphate-buffered saline solution twice for 3 min at a frequency of 30 Hz. DNA was extracted using the Nucleospin Tissue kit (Macherey Nagel), following manufacturer instructions. Each tick DNA sample was used for TBPs detection and bacterial microbiota characterisation using microfluidic PCR and 16S rRNA amplicon sequencing, respectively. Tick DNA was pre-amplified for the microfluidic PCR to have better detection of the pathogens. The PreAmp Master Mix (Fluidigm) was used according to the manufacturer's instructions. Briefly, primers were pooled combining equal volume of primers (200 nM final each). The pre-amplification was performed in a final volume of 5 µL containing 1 µL Perfecta Preamp 5x, 1.25 µL pooled primers mix, 1.5 µL distilled water and 1.25 µL DNA. The thermocycling program consisted of one cycle at 95°C for 2 min, 14 cycles at 95°C for 15 s and 4 min at 60°C. The reactions were then diluted 1:10 in Milli-Q ultrapure water. Pre-amplified DNAs were stored at -20°C.

## 2.3 | Microfluidic PCR detection of TBPs in tick DNA samples and sequencing confirmation

To detect major TBPs (32 bacterial species, five bacteria genera, two parasite genera, one parasite phylum; Table S2), the BioMark™ real-time PCR system (Fluidigm) was used for high-throughput microfluidic real-time PCR amplification using the 48.48 dynamic arrays (Fluidigm). Amplifications were performed using 6-carboxyfluorescein (FAM)- and black hole quencher (BHQ1)-labelled TaqMan probes with TaqMan Gene expression master mix following manufacturer's instructions (Applied Biosystems). PCR cycling comprised 2 min at 50°C, 10 min at 95°C, followed by 40 cycles of two-step amplification of 15 s at 95°C and 1 min at 60°C. One negative water control was included per chip. To determine if factors present in the sample could inhibit the PCR, *Escherichia coli* strain EDL933 DNA was added to each sample as an internal inhibition control, and primers and probes specific to the *E. coli eae* gene were used. For more details, please see Michelet et al. (2014). Confirmation PCR was realised for positive samples of *Rickettsia* with primers amplifying an 850-bp region of the *gltA* gene (Mediannikov et al., 2004), and *Ehrlichia* with the primers amplifying a 590-bp region of the protein (*groEL*) gene (Dahmani et al., 2017). Positive PCR products were purified and sequenced using the Sanger Sequencing method.
## 2.4 | 16S rRNA amplicon sequencing and processing of raw sequences

A single lane of Illumina MiSeq system was used to generate 251-base paired-end reads from the V4 variable region of the 16S rRNA gene using bar-coded universal primers (515F/806R) in ticks. Two extractions reagent controls were set, in which the different DNA extraction steps were performed using the same conditions as for the samples but using water as a template. DNA amplification was then performed on the extraction control in the same conditions as for any other sample. The raw 16S rRNA paired sequences obtained from tick samples were deposited at the SRA repository (Bioproject No. PRJNA865094). 16S rRNA sequences were analysed using QIIME2 pipeline (v.2019.7; Bolyen et al., 2019). The demultiplexed raw sequences (obtained in fastq files) were denoised, quality trimmed and merged using DADA2 software (Callahan et al., 2016) implemented in QIIME 2 (Bolyen et al., 2019). The amplicon sequence variants (ASVs) obtained were aligned with q2-alignment of MAFFT (Katoh et al., 2002) and used to construct a phylogeny with q2-phylogeny of FastTree 2 (Price et al., 2010). Taxonomy was assigned to ASVs using a classify-sklearn naïve Bayes taxonomic classifier based on SILVA database (release 138; Bokulich et al., 2018). Only the target sequence fragments were used in the classifier (i.e. classifier trained with the 515F/806R primers) (Ren & Wu, 2016; Werner et al., 2012). Taxa that persisted across serial fractions of the samples using QIIME 2 plugin feature-table (core-features) were considered ubiquitous. To determine the impact of *Rickettsia* taxon in the microbial community, the taxon was put to 0 in the ASVs table, and the following tests were realised with the same procedure.

## 2.5 | Microbial diversity and taxonomic differential abundance

To test for differences in bacterial diversity within and between *H. marginatum* and *R. bursa* tick samples, we conducted analyses of alpha and beta diversity. Alpha diversity measures the bacterial richness within groups, while bacterial beta diversity compares the bacterial diversity between groups. Alpha diversity and beta diversity metrics were calculated using the q2-diversity plugin in QIIME 2 (Bolyen et al., 2019). Richness and evenness were calculated using 'observed features' (DeSantis et al., 2006) and Pielou evenness index (Pielou, 1966), respectively. Differences in alpha diversity metric between groups were assessed with Kruskal-Wallis test (*p* < .05) using QIIME 2 (Bolyen et al., 2019). The beta diversity was assessed using the Bray Curtis dissimilarity index (Bray & Curtis, 1957) and compared between groups using the PERMANOVA test (*p* < .05) as implemented in QIIME 2 (Bolyen et al., 2019). Dispersion, which measures bacterial variability between samples within the population, was calculated using 'betadisper' function (script available,

File S1) and Vegan R package implemented in Rstudio (RStudio Team, 2020). The dispersion between groups was compared using the ANOVA test (*p* < .05).

The clustering analysis, which measures the similarity between the tick microbial samples, was calculated for microbial samples of the two tick species to test for a clustering pattern according to the species, as well as within each tick species individually to test for a clustering pattern according to the pathogen infection status. The cluster analysis was conducted with the Jaccard coefficient of similarity (script available, File S1) with the use of Vegan (Oksanen et al., 2021) implemented on R studio (RStudio Team, 2020). The differences in bacterial taxa abundance between the two species were tested using a t-test and performed with the ANOVA-Like Differential Expression (ALDEX2; script available, File S1) package (Fernandes et al., 2013). Only taxa with significant differences (*p* < .05) were used for representation of the differential taxa relative abundance. Relative abundance was measured as centred log ratio (clr) transformation which uses the geometric mean of the read counts in the sample. The advantage of the clr transformation is that it makes the quantification scale free and therefore comparable between conditions (Fernandes et al., 2014). The resulting data were used to construct the heatmap with the heatmap.2 function (script available, File S1), implemented in R studio environment (RStudio Team, 2020).

## 2.6 | Inference of bacterial co-occurrence networks

Here we used co-occurrence network to compare the architecture (i.e. network topology), node hierarchy (i.e. centrality values and keystone taxa), and robustness between networks of the two tick species with and without *Rickettsia*. Co-occurrence networks were built for each dataset using the bacterial taxonomic profiles at genera level. The networks allow the graphic visualisation of the bacterial community assemblies. Bacterial taxa are represented by nodes and the significant correlations between taxa are represented by edges. Analyses of significant positive (weight > 0.75) or negative (weight < -0.75) correlations were performed using the Sparse Correlations for Compositional data (SparCC) method (script available, File S1) (Friedman & Alm, 2012) implemented in R studio environment (RStudio Team, 2020). The visualisation and measures of topological features (i.e. number of nodes and edges, network diameter, modularity, average degree, weighted degree and clustering coefficient) of the networks were performed using the software Gephi 0.9.2 (Bastian et al., 2009). Topological features can provide information about the bacterial community stability and robustness. Core co-occurrence analysis was conducted to determine the strongest co-occurrence associations in the bacterial microbiota using the same technic as standard co-occurrence networks, but with a higher SparCC threshold for positive (weight > 0.9) or negative (weight < 0.9) correlations.
## 2.7 | Inference of eigenvector and betweenness centrality indexes

Eigenvector and betweenness centrality measures are used to measure the influence and importance of each node within a network. The aim of this analysis was to compare the number and position of the most central nodes between the two bacterial networks. The eigenvector centrality is a value (between 0 and 1) quantifying the influence of a node by measuring the importance of the connections of its neighbours (Ruhnau, 2000). A node with a high eigenvector value means that the node is connected with many other high-score nodes. We consider high eigenvector centrality nodes the ones that were above the median value. The betweenness centrality calculates the number of shortest paths that pass through each node in the network. A node with a high betweenness centrality value is that in the shortest path of a high number of nodes in the network. We consider high betweenness centrality nodes the ones that were above the median value. Eigenvector and betweenness centrality values of each node were exported with Gephi (Bastian et al., 2009).

## 2.8 | Comparative network analysis

The differential network represents the correlations that are different between the same taxa in two bacterial networks. The statistical network estimation was realised using Network Construction and Comparison for Microbiome Data (NetComi; script available, File S1) (Peschel et al., 2021) on Rstudio (RStudio Team, 2020). Each node represents a bacterial taxon shared by the two networks. Each edge represents a type of connection.

Jaccard index was calculated to test for dissimilarities in local centrality measures (degree, betweenness centrality, closeness centrality and eigenvector centrality) between nodes in networks of the two tick species with and without *Rickettsia*. The Jaccard index tests the similarity between sets of 'most central nodes' of networks, which are defined as those nodes with a centrality value above the empirical 75% quartile. This index expresses the similarity of the sets of most central nodes as well as the sets of hub taxa between the two networks. Jaccard index range from 0 (completely different sets) to 1 (sets equal). The two *p*-values *P* (*J* ≤ *j*) and *P* (*J* ≥ *j*) for each Jaccard's index are the probability that the observed value of Jaccard's index is 'less than or equal' or 'higher than or equal', respectively, to the Jaccard value expected at random which is calculated taking into account the present total number of taxa in both sets (Real & Vargas, 1996).

## 2.9 | Parameters to describe keystone taxa

Keystone taxa are central in microbiota functioning, having an important role for its structure and/or assembly. Here, we defined keystone taxa using three criteria, as previously reported

(Mateos-Hernández et al., 2020, 2021): (i) ubiquitousness, (ii) high eigenvector centrality (≥0.75), (iii) high relative abundance (clr value above the average).

## 2.10 | Local connectivity of *Rickettsia* in the microbial community

To understand the relation of rickettsial pathogens with the rest of the bacterial microbiota, the *Rickettsia* taxon was visualised in association with all the taxa it was positively or negatively connected to (*Rickettsia* sub-networks). The sub-networks were exported and assessed individually. The analyses were performed on Gephi (Bastian et al., 2009), the strength of the edges was presented with the SparCC weight.

## 2.11 | Network robustness in nodes removal

The network robustness analysis provides information on how resistant the network is to perturbations such as node removal. In this analysis, the proportion of removed nodes needed to reach a loss in connectivity of 0.80 was recorded for each network after directed, cascading or random removal of nodes. For comparisons between the networks with and without *Rickettsia*, a delta value was calculated (i.e. the proportion of nodes needed to reach a loss in connectivity of 0.80 in networks with *Rickettsia* minus that in networks without *Rickettsia*). Directed removal of nodes consists in removing first the nodes with higher betweenness centrality. Cascading effect consists in removing first the nodes with high betweenness centrality, but recalculated each time a node is removed. The last type is a random removal of nodes. The robustness of networks is calculated with Network Strengths and Weaknesses Analysis (NetSwan; script available, File S1; Lhomme, 2015) in Rstudio (RStudio Team, 2020). A statistical unpaired t-test was conducted between networks with and without *Rickettsia* for each species and each node removal condition with GraphPad Prism version 8.0.1 (GraphPad Software).

# 3 | RESULTS

## 3.1 | *Rickettsial pathogens infection and microbiota diversity in H. marginatum and R. bursa*

A high-throughput pathogen detection analysis was conducted on samples from *H. marginatum* and *R. bursa* ticks. The results showed the presence of at least one pathogen in 100% (23/23) of the *H. marginatum* samples, whereas only 70.8% (17/24) of *R. bursa* samples presented at least one pathogen (Table 1). All *H. marginatum* samples (n=23) were positive to *R. aeschlimannii*. For *R. bursa*, 58.3% (14/24) of the samples were positive to *Rickettsia* spp., while 54.2% (13/24) were positive to *R. aeschlimannii*, 4.2% (1/24) were positive to *R. slovaca* and 20.8% (5/24) were positive to *R. raoultii*. Other
<table><caption>TABLE 1 Pathogen occurrence in *Hyalomma marginatum* and *Rhipicephalus bursa* samples.</caption><thead><tr><th>Vector-borne pathogen(s)</th><th>H. marginatum</th><th>% H. marginatum</th><th>R. bursa</th><th>% R. bursa</th></tr></thead><tbody><tr><td>Total infected ticks (≥1 pathogen)</td><td>23</td><td>100</td><td>17</td><td>70.8</td></tr><tr><td>Rickettsia spp.</td><td>23</td><td>100</td><td>14</td><td>58.3</td></tr><tr><td>Rickettsia aesclimannii</td><td>23</td><td>100</td><td>13</td><td>54.2</td></tr><tr><td>Rickettsia slovaca</td><td>0</td><td>0</td><td>1</td><td>4.2</td></tr><tr><td>Rickettsia raoultii</td><td>0</td><td>0</td><td>5</td><td>20.8</td></tr><tr><td>Anaplasma spp.</td><td>0</td><td>0</td><td>3</td><td>12.5</td></tr><tr><td>Anaplasma marginale</td><td>0</td><td>0</td><td>3</td><td>12.5</td></tr><tr><td>Ehrlichia spp.</td><td>1</td><td>4.3</td><td>0</td><td>0</td></tr><tr><td>Ehrlichia minasensis</td><td>1</td><td>4.3</td><td>0</td><td>0</td></tr><tr><td>Theileria spp.</td><td>1</td><td>4.3</td><td>11</td><td>45.8</td></tr></tbody></table>

pathogens such as *A. marginale* (*R. bursa*, 12.5%, total 24), *Theileria* spp. (*R. bursa*, 45.8%, total 24 and *H. marginatum*, 4.3%, total 23) and *Ehrlichia minasensis* (*H. marginatum*, 4.3%, total 23) were found at lower frequencies.

In *H. marginatum* samples, single infections were more frequently detected (91.3%) than co-infections (8.7%; Table S3). In contrast, *R. bursa* samples had a lower occurrence of single infections (16.6%) and a higher occurrence of co-infections (54.2%; Table S3). These results indicate different patterns of pathogen infection between the two tick species. *Hyalomma marginatum* ticks have a high occurrence of single infection by rickettsial pathogens but low pathogen diversity (within and between samples), while *R. bursa* samples have a lower occurrence of rickettsial pathogens infections but high pathogen diversity (within and between samples).

To understand how pathogen occurrence and diversity affect the microbiota, the bacterial richness and evenness of *R. bursa* were compared to those of *H. marginatum*. The analysis showed that *R. bursa* had significantly higher bacterial richness (Kruskal-Wallis, *p* < .05, Figure 1a) and evenness (Kruskal-Wallis, *p* < .05, Figure 1b). Additionally, the bacterial community compositions of *R. bursa* and *H. marginatum* were significantly different (PERMANOVA, *p* ≤ .001, Figure 1c), while no significant differences were found in the beta dispersion (ANOVA test, *p* = .069, Figure 1c). A Jaccard distance analysis revealed that the microbial samples formed two distinct clusters: one composed mostly of *H. marginatum* samples (22/29, 75.9%) and a few *R. bursa* samples (7/29, 24.1%), and the other cluster composed mostly of *R. bursa* samples (17/18, 94.4%) and a single *H. marginatum* sample (1/18, 5.6%; Figure 1d). These clusters indicate a strong clustering pattern according to the tick species. However, within each tick species, the Jaccard distance analysis did not show any clustering pattern when considering the pathogen infection status in *H. marginatum* (Figure S1a) or *R. bursa* (Figure S1b). Overall, the composition of the bacterial microbiota was significantly different between *H. marginatum* and *R. bursa*, and the differences in pathogen infection within a species did not impact the microbial composition.

Compositional analysis revealed that the two tick species shared 857 taxa, while 81 and 242 unique taxa (genera level) were detected

in *H. marginatum* and *R. bursa*, respectively (Figure 1e). Differential relative abundance analysis identified 18 taxa with significant changes in relative abundance among the shared taxa (Figure 1f). Some taxa such as *Alcaligenes*, *Shewanella*, *Flavobacterium*, *Ralstonia*, *Cupriavidus*, *Chryseobacterium*, *Rickettsia*, *Candidatus* Midichloria and *Francisella* were more abundant in *H. marginatum* samples, while others including *Oscillospiracae* UCG-002, *Ruminococcus gauvreauii* group, *Rothia*, *Clostridia*, *Monoglobus*, *Oscillospiracae* UCG-005, *Clostridia* UCG-014, *Eubacterium coprostanoligenes* group and *Coxiella* were more abundant in *R. bursa* samples. The *Rickettsia* taxon was present in all samples with a higher mean relative abundance in *H. marginatum* (Average clr 14.55±0.49) compared to *R. bursa* (Average clr 13.61±0.62) (*p* < .05). The contrasting patterns of rickettsial pathogen occurrence (Table S3), relative abundance of the *Rickettsia* taxon, and of alpha diversity (Figure 1a) between the two tick species suggests an impact of *Rickettsia* on the tick bacterial microbiota.

## 3.2 | Centrality and connectivity of *Rickettsia* in the *H. marginatum* and *R. bursa* microbial networks

The importance of *Rickettsia* in the bacterial networks was assessed by examining its centrality within each network. Specifically, we checked whether *Rickettsia* was among the most central nodes or connected with keystone taxa. The *Rickettsia* taxon was not among the nodes with high eigenvector centrality values in *H. marginatum* (*n* = 11 nodes, Figure 2a), and *R. bursa* (*n* = 15 nodes, Figure 2b), or the nodes with high betweenness centrality in *H. marginatum* (*n* = 2 nodes, Figure 2c) and *R. bursa* (*n* = 1 node, Figure 2d) networks. Most of the nodes had eigenvector centrality and betweenness centrality values below the mean values (80.4% for *H. marginatum*, Figure 2e and 79.7% for *R. bursa*, Figure 2f). To further explore centrality distribution between the two networks, we conducted a Jaccard index comparison test on local centrality measures. Jaccard index values for degree, hub taxa, closeness and eigenvector centrality between the two species networks was 0.40, which was higher than expected by random (*P* (≥Jacc) < .05, Table 2). However, the Jaccard value for
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** Comparison of microbial diversity and taxa bacterial abundance between *Hyalomma marginatum* and *Rhipicephalus bursa*. (a) Comparison of alpha diversity with observed features index for *H. marginatum* (purple) and *R. bursa* (green) samples (*p* = .022*). (b) Comparison of alpha diversity with Pielou evenness index for *H. marginatum* (purple) and *R. bursa* (green) samples (*p* = .011*). (c) Comparison of beta diversity with Bray Curtis dissimilarity index for *H. marginatum* (purple) and *R. bursa* (green) samples. Small circles represent samples, and ellipses represent centroid position for each group (*F* = 11.39; *p* = .0001**; stress = 0.1244). This test use Principal Coordinate Analysis (PCoA), it is used to explore and to visualise variability in a microbial community. ANOVA test was performed and showed that beta dispersion of the two sets of samples is not significantly different (*p* = .0686). (d) Jaccard clusterisation of *H. marginatum* and *R. bursa* samples. The Jaccard distance is represented between 0 and 2, the lines are proportional to this distance. Each sample is represented by a violet dot for *H. marginatum* and a green one for *R. bursa*. (e) Venn diagram displaying the comparison of taxa composition for all taxa. The purple and green circles represent *H. marginatum* and *R. bursa* samples, respectively. Numbers represent the amount of taxa found in each dataset and those shared by the two groups. (f) Dendrogram heatmap resulting from the heatmap.2 function implemented on R studio (RStudio Team, 2020). Taxa were clustered based on relative abundance (calculated as clr transformed values). Each column represents the CLR values for bacterial taxa per sample and per group. Each line represents bacterial taxa with significant changes between the two datasets. clr value range from -15 (red) to 15 (green).</figcaption></figure>

betweenness centrality did not deviate from random distribution (Table 2), suggesting that the two networks were mostly equivalent in terms of node centrality distribution.

To gain a deeper understanding of the distribution of node centrality in the two networks, we examined the 10 taxa with the highest eigenvector, betweenness and degree centrality values.
<figure><img src="image_5.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(g)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(h)</figcaption></figure>

<figure><img src="image_12.png" /></figure>
FIGURE 2 Distribution of node centrality in *Hyalomma marginatum* and *Rhipicephalus bursa* co-occurrence networks. Co-occurrence networks of *H. marginatum* (a) and (c) and *R. bursa* (b) and (d). Network (a) and (b) represent eigenvector centrality. Node colour ranges from white (0) to purple (1). Node size is proportional with eigenvector centrality measure. Network (c) and network (d) represent betweenness centrality. Node colour ranges from white (0) to purple (1000 for (c) and 2500 for (d)). Node size is proportional with betweenness centrality measure. Eigenvector centrality and Betweenness centrality ratio is represented for *H. marginatum* (e) and *R. bursa* (f). Each black dot represents a taxon. The red dots represent *Rickettsia* taxon. Direct neighbours of *Rickettsia* for (g) *H. marginatum* and (h) *R. bursa* microbial networks. Blue and red edges represent positive and negative co-occurring relations respectively, the number is the co-occurring association strength (the SparCC value).

<table><caption>TABLE 2 Jaccard index for *Hyalomma marginatum* and *Rhipicephalus bursa* networks.</caption><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">H. marginatum vs. R. bursa</th></tr><tr><th>Jacc<sup>a</sup></th><th>P (≤Jacc)</th><th>P (≥Jacc)</th></tr></thead><tbody><tr><td>Degree</td><td>0.40</td><td>.97</td><td>.04*</td></tr><tr><td>Betweenness centrality</td><td>0.35</td><td>.64</td><td>.45</td></tr><tr><td>Closeness centrality</td><td>0.40</td><td>.97</td><td>.04*</td></tr><tr><td>Eigenvector centrality</td><td>0.40</td><td>.97</td><td>.04*</td></tr><tr><td>Hub taxa</td><td>0.40</td><td>.97</td><td>.04*</td></tr></tbody></table>

Interestingly, 50% of the most central nodes were similar between the two networks in terms of eigenvector and degree centralities (Table S4), while only 10% were similar for betweenness centrality (Table S4). Despite *Rickettsia* being among the low-ranked nodes in both networks (Figure 2e,f), it is worth noting that its eigenvector centrality value of 0.11 and betweenness centrality value of 3.7 in the *H. marginatum* network were higher than those in *R. bursa* network, where *Rickettsia* had eigenvector centrality=0.03, and betweenness centrality=0.

Furthermore, our analysis of keystone taxa identification and local connectivity analysis revealed that *Rickettsia* showed a positive connection with the only keystone taxon (*Cutibacterium*) identified in *H. marginatum* (Figure 2g). However, no connection was found between the keystone taxa in *R. bursa* (*Staphylococcus* and *Jeotgalicoccus*) and *Rickettsia* in *R. bursa* network (Figure 2h). The results indicated similar distribution of node centrality in the two networks and slightly higher importance of *Rickettsia* in the bacterial microbiota of *H. marginatum*.

To further investigate the importance of *Rickettsia* in the networks, we conducted an in silico experiment where *Rickettsia* was removed from the networks to observe its effect on node centrality distribution. The comparison of centrality measures between networks with and without *Rickettsia* in *H. marginatum* or *R. bursa* showed that the removal of *Rickettsia* resulted in Jaccard indexes very close to 1 for all local centrality measures (Table 3). This suggests a high similarity in the proportion of central nodes in the networks with and without *Rickettsia* within each tick species. However, when comparing Jaccard values between *H. marginatum* and *R. bursa* networks without *Rickettsia*, we observed a random distribution of all centrality measures (Table 4).

To further analyse the taxonomic composition of the two tick species with and without *Rickettsia*, we identified the top 10 taxa with the highest eigenvector, betweenness and degree centrality values in networks without *Rickettsia*. Interestingly, the top taxa differed significantly between networks with and without *Rickettsia* for each node centrality metric (Table S4). To determine whether the loss of similarity in node centrality distribution was associated to the specific removal of *Rickettsia* taxon, we also removed one taxon in each network with similar centrality values as *Rickettsia* (referred as 'equivalent taxa') and compared the centrality distribution between networks. The removal of equivalent taxa did not significantly alter the node centrality distribution in the *H. marginatum* and *R. bursa* networks (Table S5).

These findings suggest that despite its low centrality values, the presence of *Rickettsia* in the networks played a crucial role in structuring node centrality traits shared by the microbiota of the two sympatric tick species. In other words, the similarity of node centrality distribution between the networks depended on the presence of *Rickettsia* in the networks.

## 3.3 | Influence of *Rickettsia* on the assembly and robustness of *H. marginatum* and *R. bursa* microbiota

*Rickettsia* plays a crucial role in the node centrality traits shared by the two tick species. To understand the role of *Rickettsia* in the community assembly, we visually inspected the networks of *H. marginatum* (Figure 3a) and *R. bursa* (Figure 3b). Both networks displayed two principal modules with negative co-occurring interactions between their nodes. Other scattered taxa marginally connected or not with the main modules were also found. Although the two networks had a similar number of nodes, the *H. marginatum* network had fewer connected nodes and edges compared to the *R. bursa* network (Table 5). However, both networks had the same proportion of positives and negatives edges (Table 5). Further analysis revealed variations in the correlation between shared taxa in the two networks (Figure 3c).

Removal of *Rickettsia* from the *H. marginatum* microbial network resulted in a loss of two connected nodes but gained of 13 edges, mostly negatives (connected nodes n=95, edges n=344), compared to the network with *Rickettsia* taxon. In contrast, removing *Rickettsia* from the *R. bursa* network increased the number of connected nodes by 6 and the number of edges by 13, mostly positives (connected nodes n=151, edges n=619). Both networks shared 70 nodes and

<footnote>*p < .05.</footnote>

<footnote><sup>a</sup>Jaccard index.</footnote>
<table><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">H. marginatum with vs. without Rickettsia</th><th colspan="3">R. bursa with vs. without Rickettsia</th></tr><tr><th>Jacc<sup>a</sup></th><th>P (≤Jacc)</th><th>P (≥Jacc)</th><th>Jacc</th><th>P (≤Jacc)</th><th>P (≥Jacc)</th></tr></thead><tbody><tr><td>Degree</td><td>0.91</td><td>1</td><td>0***</td><td>0.91</td><td>1</td><td>0***</td></tr><tr><td>Betweenness centrality</td><td>0.82</td><td>1</td><td>0***</td><td>0.84</td><td>1</td><td>0***</td></tr><tr><td>Closeness centrality</td><td>0.91</td><td>1</td><td>0***</td><td>0.91</td><td>1</td><td>0***</td></tr><tr><td>Eigenvector centrality</td><td>0.91</td><td>1</td><td>0***</td><td>0.91</td><td>1</td><td>0***</td></tr><tr><td>Hub taxa</td><td>0.91</td><td>1</td><td>0***</td><td>0.91</td><td>1</td><td>0***</td></tr></tbody></table>

<caption>TABLE 3 Jaccard index for *Hyalomma marginatum* and *Rhipicephalus bursa* networks with and without *Rickettsia* taxon.</caption>

<table><caption>TABLE 4 Jaccard index for *Hyalomma marginatum* and *Rhipicephalus bursa* networks without *Rickettsia* taxon.</caption><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">H. marginatum without Rickettsia vs. R. bursa without Rickettsia</th></tr><tr><th>Jacc<sup>a</sup></th><th>P (≤Jacc)</th><th>P (≥Jacc)</th></tr></thead><tbody><tr><td>Degree</td><td>0.39</td><td>.95</td><td>.07</td></tr><tr><td>Betweenness centrality</td><td>0.32</td><td>.46</td><td>.63</td></tr><tr><td>Closeness centrality</td><td>0.39</td><td>.95</td><td>.07</td></tr><tr><td>Eigenvector centrality</td><td>0.39</td><td>.95</td><td>.07</td></tr><tr><td>Hub taxa</td><td>0.39</td><td>.95</td><td>.07</td></tr></tbody></table>

the *H. marginatum* network had 27 unique nodes, while the *R. bursa* network had 75 unique nodes (Figure S2a). When *Rickettsia* was removed, the two networks shared 68 nodes, and the *H. marginatum* network still had 27 unique nodes, while the *R. bursa* network had 83 unique nodes (Figure S2b).

We extended the comparison to the 'core bacterial microbiota', which consisted of nodes with strong correlation (SparCC ≥ 0.90). The networks of *H. marginatum* (Figure S2c) and *R. bursa* (Figure S2d) shared a core composed of three nodes: *Psychrobacter*, *Atopostipes* and *Jeotgalicoccus*. However, these nodes did not share the same strong co-occurring interaction between the two species. In addition, both *H. marginatum* and *R. bursa* had an extra node connected to three common nodes (*Jeotgalibaca* for *H. marginatum* and *Aerococcaceae* - uncultured for *R. bursa*). Furthermore, *R. bursa* had an extra core of four nodes not shared by *H. marginatum* (*Cutibacterium*, *Lawsonella*, *Desulfobulbaceae* -uncultured and *Sva0081 sediment group*) (Figure S2c,d). Removing *Rickettsia* resulted in the gain of one taxon (*Aerococcaceae* - uncultured) in the 'core bacterial microbiota' of *H. marginatum* (Figure S2e) and the loss of two taxa (*Desulfobulbaceae* - uncultured and *Sva0081 sediment group*) from the 'core bacterial microbiota' of *R. bursa* (Figure S2f). The removal of *Rickettsia* also affected the general structure of the networks by increasing the number of edges.

Next, we compared the robustness of the networks with and without *Rickettsia*. Higher delta values indicated less robustness in networks without *Rickettsia*, meaning fewer nodes were needed to

reach a loss in connectivity of 0.8. Directed removal of nodes from *H. marginatum* and *R. bursa* networks with and without *Rickettsia* resulted in deltas of 2.7% (Figure S3a) and 4.8% (Figure S3b), respectively. Cascading removal of nodes resulted in deltas of 1.8% for *H. marginatum* (Figure S3c) and -0.3% for *R. bursa* (Figure S3d). Random removal of nodes showed deltas of 0.1% for *H. marginatum* (Figure S3e) and 1.2% for *R. bursa* (Figure S3f). Only the direct removal of nodes for *R. bursa* networks with versus without *Rickettsia* significantly different (unpaired t-test, *p* < .01). Therefore, removing *Rickettsia* followed by directed attack had the strongest negative effect on the robustness of both networks by reducing the proportion of nodes required to reach a loss in connectivity of 0.8.

## 3.4 | Rickettsia as a driver of microbial community assembly in *H. marginatum* and *R. bursa*

We conducted an evaluation to assess the potential influence of *Rickettsia* on the arrangement of nodes within the two primary microbial communities (computed modules) identified in *H. marginatum* (Figure 3a) and *R. bursa* (Figure 3b) bacterial microbiota. Comparison of bacterial nodes in the two major modules of the two species networks revealed a significant number of shared taxa. Specifically, one module in *H. marginatum* shared 22 taxa with a module in *R. bursa* (Figure 4a), prompting us to designate these modules as equivalent communities, namely *H. marginatum* (Hm) module 1 (HmM1) and *R. bursa* (Rb) module 1 (RbM1). The other two modules shared 12 nodes (Figure 4a) and were labelled as HmM2 and RbM2.

Each module also had unique nodes, 23 in *H. marginatum* (11 in HmM1 and 12 in HmM2) and 60 in *R. bursa* (22 in RbM1 and 38 in RbM2). Notably, *Rickettsia* was embedded within HmM1 and RbM2. Upon removing *Rickettsia* from the *H. marginatum* network, we observed the absence of 17 taxa (14 in HmM1 and 3 in HmM2) in the two modules. However, the removal of *Rickettsia* did not impact 58 taxa (19 in HmM1 and 22 in HmM2), and 24 taxa (8 in HmM1 and 16 in HmM2) were gained by the modules (Figure 4b). Similarly, removing *Rickettsia* from *R. bursa* network resulted in the disappearance of six taxa (4 in RbM1 and 2 in RbM2), the appearance of five taxa (4 in

<footnote>***p < .0001.</footnote>

<footnote><sup>a</sup>Jaccard index.</footnote>

<footnote><sup>a</sup>Jaccard index.</footnote>
<figure><img src="image_13.png" /><figcaption>FIGURE 3 Co-occurrence networks of *Hyalomma marginatum* and *Rhipicephalus bursa*. Co-occurrence networks of (a) *H. marginatum* and (b) *R. bursa* samples. Nodes correspond to taxa (family or genera level), and connecting edges indicate significant connection between them (*p* < .01**). Only nodes with at least one significant correlation are represented. Node colours are based on modularity class metric and equal colour means modules of co-occurring taxa. Edges represent positive (blue) or negative (red) correlations (SparCC > 0.75 or <-0.75). (c) Differential network between *H. marginatum* and *R. bursa* co-occurrence networks. White nodes represent taxa, edges represent differential of co-occurring connection between the two species. Light green and red means the taxa are positively or negatively connected in both networks. Dark green means the taxa are positively connected in *H. marginatum* and not connected in *R. bursa* networks. Dark blue means the taxa are positively connected in *R. bursa* and not connected in *H. marginatum* networks. Cyan means that taxa are positively connected in *H. marginatum* network and negatively connected in *R. bursa* network. Purple means that taxa are positively connected in *R. bursa* network and negatively connected in *H. marginatum* network.</figcaption></figure>

<table><caption>TABLE 5 Topological features for the perturbed and normal networks of *Hyalomma marginatum* and *Rhipicephalus bursa*.</caption><thead><tr><th>Topological features</th><th>H. marginatum</th><th>R. bursa</th><th>H. marginatum without Rickettsia</th><th>R. bursa without Rickettsia</th></tr></thead><tbody><tr><td>Total nodes</td><td>1256</td><td>1279</td><td>1255</td><td>1278</td></tr><tr><td>Connected nodes</td><td>97</td><td>145</td><td>95</td><td>151</td></tr><tr><td>Edges</td><td>331</td><td>606</td><td>344</td><td>619</td></tr><tr><td>Positives</td><td>237 (72%)</td><td>433 (71%)</td><td>239 (69%)</td><td>447 (72%)</td></tr><tr><td>Negatives</td><td>94 (28%)</td><td>173 (29%)</td><td>105 (31%)</td><td>172 (28%)</td></tr><tr><td>Modularity</td><td>1.15</td><td>1.10</td><td>1.27</td><td>1.06</td></tr><tr><td>Network diameter</td><td>7</td><td>10</td><td>8</td><td>9</td></tr><tr><td>Average degree</td><td>6.83</td><td>8.36</td><td>7.24</td><td>8.20</td></tr><tr><td>Weighted degree</td><td>2.35</td><td>2.86</td><td>2.26</td><td>2.89</td></tr><tr><td>Clustering coefficient</td><td>0.53</td><td>0.59</td><td>0.56</td><td>0.62</td></tr></tbody></table>

RbM1 and 1 in RbM2), while 89 taxa (41 in RbM1 and 48 in RbM2) remained unaffected (Figure 4c). These alterations in the modules, following *Rickettsia* removal, led to similarities between HmM1 and RbM1 (22 shared taxa with *Rickettsia* vs. 23 without *Rickettsia*) (Figure 4a,d), as well as between HmM2 and RbM2 (12 shared taxa with *Rickettsia* vs. 21 without *Rickettsia*; Figure 4a,d). These findings suggested that the composition of HmM1 and RbM1, as well as HmM2 and RbM2, became more alike after removal of *Rickettsia* from the networks.

Next, we examined the differential connectivity of the common taxa between HmM1 and RbM1, and between HmM2 and RbM2, to determine if the taxa present in both modules of the two species exhibited the same co-occurrence pattern regardless of the species. Surprisingly, we observed that most of those connections differed

between two groups (Figure 4e). Although the communities within the major modules shared similar taxa in *H. marginatum* and *R. bursa*, allowing us to draw parallels between HmM1 and RbM1, as well as HmM2 and RbM2, the way the taxa were interconnected differed between the two species. The presence of *Rickettsia* in the communities significantly influenced the assemblage of the communities, leading to divergence between the bacterial microbiota of *H. marginatum* and *R. bursa*.

4 | DISCUSSION

The capacity of TBPs to modulate the diversity, bacterial composition and microbial structure of tick host microbiota has been
<figure><img src="image_14.png" /><figcaption>FIGURE 4 Proportion of unique and shared nodes in the two major communities of *H. marginatum* and *R. bursa* microbiota. (a) Taxa shared between *H. marginatum* (Hm) (with purple outline) module 1 (HmM1) and *R. bursa* (Rb) (with green outline) module 1 (RbM1) in pink and Hm module 2 (HmM2), and Rb module 2 (RbM2) in green. (b) Taxa shared between HmM1 with and without *Rickettsia* (in pink) and HmM2 with and without *Rickettsia* (in green). (c) Taxa shared between RbM1 with and without *Rickettsia* (in pink) and RbM2 with and without *Rickettsia* (in green). (d) Taxa shared between HmM1 and RbM1 without the presence of Rickettsial taxon (in pink), and HmM2 and RbM2 without the presence of Rickettsial taxon (in green). (e) Differential network between the two principal modules of *H. marginatum* and *R. bursa* co-occurrence networks. Node colours represent the two modules. Edges represent strong positive (SparCC > 0.75) or negative (SparCC < -0.75) co-occurrence interactions. Light green and red means the taxa are positively or negatively connected in both networks. Dark green and orange means the taxa are positively or negatively connected in *H. marginatum* and not connected in *R. bursa* networks. Dark blue and dark purple means the taxa are positively or negatively connected in *R. bursa* and not connected in *H. marginatum* networks.</figcaption></figure>

established (Abraham et al., 2017; Estrada-Peña et al., 2020; Hamilton et al., 2021; Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022). However, their impact on community assembly and emergent properties of microbial networks remains unclear. TBPs can potentially restructure the tick host microbial community through microbe-host or microbe-microbe complex (Abraham et al., 2017; Lejal et al., 2021; Narasimhan et al., 2022). In this study, we employed high-throughput pathogen detection, next-generation sequencing (NGS), and an experimental network approach to investigate in silico the impact of the *Rickettsia* taxon on the assembly of *H. marginatum* and *R. bursa* bacterial microbiota.

Among the *Rickettsia* species detectable in this study (*R. conorii*, *R. slovaca*, *R. massiliae*, *R. helvetica*, *R. aeschlimannii*, *R. raoultii* and *R. felis*), *R. aeschlimannii* was found to be overrepresented in both *H. marginatum* (present in 100% of the samples) and *R. bursa* (present in over 50% of the samples, i.e. 54.2%) ticks. None of the other *Rickettsia* species were detected in *H. marginatum*, while *R. slovaca* (4.2%) and *R. raoultii* (20%) were detected at lower occurrence in *R. bursa*. Thus, the network associations observed in this study were assumed to be

primarily associated with the presence of the high-frequently found species *R. aeschlimannii*, yet the possibility of spurious signals resulting from the presence of *R. slovaca* and/or *R. raoultii* in *R. bursa* cannot be excluded. In addition to pathogens, there is also a possibility that *Rickettsia* endosymbionts commonly found in ticks (Nováková & Šmajs, 2018) were present in the ticks and detected by NGS, yet *Francisella*-like and *Midichloria*, rather than *Rickettsia*, are the endosymbionts frequently found in *H. marginatum* (Azagi et al., 2017; Buysse et al., 2021). The high rates of *R. aeschlimannii* infection in *H. marginatum* (Azagi et al., 2017; Wallménius et al., 2014) may be associated with this bacterium being an endosymbiont, although no vertical transmission has been reported for *R. aeschlimannii* in ticks (Nováková & Šmajs, 2018). On the other hand, natural populations of *R. bursa* have been more frequently associated with *Coxiella*-like than *Francisella*-like endosymbionts (Brinkmann et al., 2019; Papa et al., 2017; Raele et al., 2015).

Subsequently, we perturbed the co-occurrence bacterial networks by removing *Rickettsia* in silico and compared the topology, connectivity, local centrality measures, robustness and community assembly between the perturbed and the original
(unperturbed) networks. In silico node removal has been used previously to assess the influence of microorganisms on plant microbiota properties (Agler et al., 2016). Specifically, the removal of hub taxa in the *Arabidopsis thaliana* microbiota affected more edges than removal of non-hub species in the networks (Agler et al., 2016). Furthermore, the significance of two hub taxa was confirmed through host colonisation experiments and interaction assays (Agler et al., 2016), demonstrating the validity of in silico node removal in silico as a tool to predict ecosystem behaviour (Röttjers & Faust, 2018). In our study, we observed that the removal of *Rickettsia* had a notable impact on the assembly of two major bacterial communities. Without *Rickettsia*, the distribution of centrality values and the presence of bacterial nodes in the microbial network of our two tick species were no longer similar. This suggests that *Rickettsia* plays a crucial role in microbiome interactions in both *H. marginatum* and *R. bursa*, which may be partially explained by its interactions with central bacteria in the networks.

However, it is important to acknowledge that while node removal analysis is a useful and powerful methodology, it may not accurately predict real-life system responses since the network represents a simplified and static representation of the ecosystem (Röttjers & Faust, 2018). Nevertheless, we found that the 'core bacterial microbiota' remained conserved in both *H. marginatum* and *R. bursa*, and the removal of *Rickettsia* did not affect the connectivity of the nodes within the core. Similar findings have been observed in different cichlid fish species from various continents, where the 'core bacterial microbiota' remained conserved (Riera & Baldo, 2020). This suggests that core microbiota bacteria are conserved among closely related species due to their critical roles in maintaining the proper functioning of microbial communities (Weese, 2013).

Disturbance of tick microbial networks by the presence of a rickettsial pathogen (*R. helvetica*) was previously reported (Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022). This infection was found to be associated with reduced bacterial diversity in *I. ricinus* ticks and loss in global network connectivity, as well as local connectivity of the *Rickettsia* taxon (Maitre, Wu-Chuang, Aželytė, et al., 2022; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022). However, the structural role of *R. helvetica* in the *R. helvetica*-infected *I. ricinus* network was not tested in that study. In our research, we discovered a positive correlation between *Rickettsia* and the predicted keystone taxon *Cutibacterium* in the undisturbed *H. marginatum* network. Keystone taxa are known to have a major influence on the structure and function of microbiota at specific spatial or temporal contexts (Banerjee et al., 2018; Wu-Chuang, Obregon, Mateos-Hernández, et al., 2022). Certain keystone bacteria have been associated with the recovery of human gut microbiota after antibiotic exposure (Gibbons, 2020), and in another study, four keystone bacteria were found to maintain the functional diversity of the *I. scapularis* microbiota under intense heat stress (Wu-Chuang et al., 2021; Wu-Chuang, Obregon, Estrada-Peña, et al., 2022). Through the privileged connectivity to keystone taxa in *H. marginatum* bacterial microbiota, *Rickettsia* may

exert an influence microbial community assembly and resiliency to attack (i.e. robustness), despite having low centrality. This phenomenon of a low number of associations linked mostly to highly connected hub nodes has been previously demonstrated as a microbial mechanism employed by *Roseofilum reptotaenium* to alter microbiome interactions associated with Black Band Disease in corals (Meyer et al., 2016).

The concept of robustness, which refers to the resistance of a network, can be studied using percolation theory (Cohen et al., 2000). Percolation theory provides insights into how information can flow between nodes in a network (Röttjers & Faust, 2018). In this study, we applied percolation theory to assess network robustness by measuring loss in connectivity following directed, cascading and random attacks (Röttjers & Faust, 2018). Our results showed that the removal of *Rickettsia* slightly decreased network robustness in both tick species across all three types of attacks, except for cascading attack in *R. bursa*. Interestingly, a previous study found that *A. phagocytophilum* infection also reduced the robustness of *I. scapularis* microbial networks, although the impact of pathogen infection on network robustness was marginal compared to that caused by anti-tick immunity (Estrada-Peña et al., 2020). These findings suggest that while pathogen infection (e.g. *A. phagocytophilum*), can be detected in the robustness test, its effect on bacterial community assembly is limited compared with other disturbance factors like anti-tick immunity (Estrada-Peña et al., 2020), or anti-microbiota vaccines (Mateos-Hernández et al., 2021). Furthermore, a significant reduction of network robustness was observed after directed attacks in the microbiota of mice exposed to antibiotics and fed a high-fat diet, as compared to untreated mice on the same diet (Mahana et al., 2016). The authors of this study suggested that specific antibiotics could target keystone taxa leading to ecosystem collapse (Mahana et al., 2016). Thus, network robustness appears to be a property that can reflect microbial community failures in various animal species, ranging from arthropod (Estrada-Peña et al., 2020; Mateos-Hernández et al., 2021) to mammals (Mahana et al., 2016).

Reducing network robustness by anti-microbiota vaccines presents a novel strategy to impact vector microbiota, as demonstrated by Mateos-Hernández et al. (2021). This approach aims to induce infection-refractory states in ticks or other vectors, as highlighted by Maitre, Wu-Chuang, Aželytė, et al. (2022) and Maitre, Wu-Chuang, Mateos-Hernández, et al. (2022), ultimately leading to a reduction or blockade of vector-borne pathogen transmission (Aželytė et al., 2022). Previous research showed that anti-microbiota vaccines that targets keystone taxa can modulate the tick microbiota (Mateos-Hernández et al., 2020, 2021). Our own findings indicate that *Cutibacterium*, *Jeotgalicoccus* and *Staphylococcus* are the key taxa for *H. marginatum* and for *R. bursa*, respectively. *Staphylococcus* is known to be highly abundant in the microbiota of various tick species (Wang et al., 2021), and has been associated with increased tick survival and reproduction (Riera & Baldo, 2020; Weese, 2013). This suggests that targeting *Staphylococcus* with an anti-microbiota vaccine could affect the transmission of rickettsial pathogens in *R. bursa*.
Although previous studies in mice and birds have shown no adverse effect or significant modulation of host microbiota associated with anti-microbiota vaccination (Mateos-Hernández et al., 2020, 2021), it is important to consider that *Cutibacterium* and *Staphylococcus* genera contain species commonly found in the skin microbiota of mammals (Fernández et al., 2020). For example, *Cutibacterium* and *Staphylococcus* have been identified in the microbiota of cattle udder (Gryaznova et al., 2021), and *Jeotgalicoccus*, a genus in the family as *Staphylococcus*, has been found in cattle nasopharynx (Timsit et al., 2018). To ensure safety, it is ideal to formulate the vaccine using protein variants exclusively found in tick-associated bacteria, rather than those present in host-associated microbiota.

Network analysis is a powerful tool that has its limitations in the present study. While microbial networks provide valuable information about the microbiota at the time of tick collection and sample generation; they do not capture the dynamic states of the microbiota that may occur in nature. Consequently, this analysis does not allow assessing the order of bacterial arrival to the tick microbiota. As a result, our analysis cannot conclusively determine whether rickettsial pathogens were acquired during the last blood feeding when ticks were collected or in earlier life stages. Moreover, it is possible that rickettsial pathogens entered the ticks simultaneously with other bacteria during blood feeding. This makes impossible to assess how the community assembly order of arrival affects the microbial community. Furthermore, we are unable determine whether co-occurring bacteria interact directly within the same organs or indirectly from distant tick body locations. This lack of information limits our understanding of the intricate dynamics within the tick microbiota. However, despite these limitations, the field of microbial ecology in ticks can still benefit greatly from the application of network analysis. It provides valuable insights into the interactions and relationships between microbial species within the tick microbiota, allowing researchers to gain a better understanding of tick-borne diseases and potential control strategies.

In summary, we found single and multiple infections associated with low and high bacterial microbiota diversity in two sympatric tick species in Corsica, *H. marginatum* and *R. bursa*. Differences in bacterial diversity and composition may be the summary of stochastic acquisition (i.e. ticks questing in different environments and/or feeding on different animals are exposed to different bacterial species) and colonisation (i.e. exposure to a bacterial species does not guarantee colonisation) events as well as pathogen-associated modulation. Despite differences in diversity and composition, there were conserved patterns of bacterial assembly, since interactions in the 'core bacterial microbiota' and two major communities (i.e. HmM1 and HmM2, and RbM1 and RbM2) were highly conserved in both tick species, revealing that colonisation of different bacterial species in similar context encode ecological information about interactions. Preferential connections of *Rickettsia* in the networks make this taxon a driver of tick microbial community assembly, since *Rickettsia* node removal caused

the re-organisation of connectivity in the two major communities shared by the two tick species. Despite community re-organisation and re-distribution of local centrality measures, network robustness was only slightly affected after *Rickettsia* node removal, suggesting that tick microbiota assembly minimises the centrality of TBPs in the networks to conserve robustness. Co-evolution between vectors, pathogens and microbiota would explain the low evolutionary pressure of pathogen infection on microbiota robustness. Application of network biology will enhance our understanding of microbiota responses to pathogens and endosymbionts. Knowing the intricacies of bacterial community assembly in ticks may prove helpful in determining key players for pathogen colonisation, which may guide interventions such as anti-microbiota vaccines to improve human and animal health.

## AUTHOR CONTRIBUTIONS

ACC and AM conceived the idea and designed experiments. AM, LMH, EPS, AFS, and VC conducted the experiments. ACC, AWC, and DO supervised data analysis. AM, and ACC analysed and/or interpreted the results. AM visualised the results. ACC, SM, JCP, and AF provided resources and input on the experiments and data interpretation. AM, and ACC drafted the first version of the manuscript. AM, AWC, LMH, EPS, AFS, VC, SM, JCP, AF, DO and ACC reviewed, edited, and approved the manuscript.

## ACKNOWLEDGEMENTS

UMR BIPAR is supported by the French Government's Investissement d'Avenir program, Laboratoire d'Excellence 'Integrative Biology of Emerging Infectious Diseases' (grant no. ANR-10-LABX-62-IBEID). Alejandra Wu-Chuang is supported by Programa Nacional de Becas de Postgrado en el Exterior 'Don Carlos Antonio López' (grant no. 205/2018). Apolline Maitre is supported by the 'Collectivité de Corse', grant: 'Formations superieures' (SGCE - RAPPORT No 0300).

## CONFLICT OF INTEREST STATEMENT

The authors declare no competing interests.

## DATA AVAILABILITY STATEMENT

The datasets generated and analysed during the current study are available on the SRA repository (Bioproject No. PRJNA865094).

## BENEFIT-SHARING STATEMENT

Research collaboration was developed with scientists from the Corsica providing genetic samples. All collaborators are included as co-authors. The research addresses a priority concern, in this case the incidence of tick-borne pathogens in Corsica. More broadly, our group is committed to international scientific partnerships, as well as institutional capacity building.

## ORCID

Dasiel Obregón https://orcid.org/0000-0002-5786-1114

Alejandro Cabezas-Cruz https://orcid.org/0000-0002-8660-730X
REFERENCES

Abraham, N. M., Liu, L., Jutras, B. L., Yadav, A. K., Narasimhan, S., Gopalakrishnan, V., Ansari, J. M., Jefferson, K. K., Cava, F., Jacobs-Wagner, C., & Fikrig, E. (2017). Pathogen-mediated manipulation of arthropod microbiota to promote infection. *Proceedings of the National Academy of Sciences of the United States of America*, 114(5), E781-E790.

Aderem, A. (2005). Systems biology: Its practice and challenges. *Cell*, 121(4), 511–513.

Agler, M. T., Ruhe, J., Kroll, S., Morhenn, C., Kim, S.-T., Weigel, D., & Kemen, E. M. (2016). Microbial hub taxa link host and abiotic factors to plant microbiome variation. *PLoS Biology*, 14(1), e1002352.

Azagi, T., Klement, E., Perlman, G., Lustig, Y., Mumcuoglu, K. Y., Apanaskevich, D. A., & Gottlieb, Y. (2017). *Francisella*-like endosymbionts and *Rickettsia* species in local and imported *Hyalomma* ticks. *Applied and Environmental Microbiology*, 83(18), e01302-17.

Aželytė, J., Wu-Chuang, A., Žiegytė, R., Platonova, E., Mateos-Hernandez, L., Maye, J., Obregon, D., Palinauskas, V., & Cabezas-Cruz, A. (2022). Anti-microbiota vaccine reduces Avian Malaria infection within mosquito vectors. *Frontiers in Immunology*, 13, 841835.

Banerjee, S., Schlaeppi, K., & van der Heijden, M. G. A. (2018). Keystone taxa as drivers of microbiome structure and functioning. *Nature Reviews. Microbiology*, 16(9), 567–576.

Bastian, M., Heymann, S., & Jacomy, M. (2009). *Gephi: An open source software for exploring and manipulating networks*. WebAtlas.

Berg, G., Rybakova, D., Fisher, D., Cernava, T., Champonier Vergès, M.-C., Charles, T., Chen, X., Cocolin, L., Eversole, K., Corral, G. H., Kazou, M., Kinkel, L., Lange, L., Lima, N., Loy, A., Macklin, J. A., Maguin, E., Mauchline, T., McClure, R., ... Schloter, M. (2020). Microbiome definition re-visited: Old concepts and new challenges. *Microbiome*, 8, 103.

Binetruy, F., Dupraz, M., Buysse, M., & Duron, O. (2019). Surface sterilization methods impact measures of internal microbial diversity in ticks. *Parasites & Vectors*, 12(1), 268.

Bokulich, N. A., Kaehler, B. D., Rideout, J. R., Dillon, M., Bolyen, E., Knight, R., Huttley, G. A., & Gregory Caporaso, J. (2018). Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2's q2-feature-classifier plugin. *Microbiome*, 6(1), 90.

Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., Al-Ghalith, G. A., Alexander, H., Alm, E. J., Arumugam, M., Asnicar, F., Bai, Y., Bisanz, J. E., Bittinger, K., Brejnrod, A., Brislawn, C. J., Brown, C. T., Callahan, B. J., Caraballo-Rodríguez, A. M., Chase, J., ... Caporaso, J. G. (2019). Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. *Nature Biotechnology*, 37(8), 852–857.

Bray, J. R., & Curtis, J. T. (1957). An ordination of the upland forest communities of Southern Wisconsin. *Ecological Monographs*, 27(4), 325–349.

Brinkmann, A., Hekimoğlu, O., Dincer, E., Hagedorn, P., Nitsche, A., & Ergünay, K. (2019). A metagenomic survey of ticks reveals pathogenic *Rickettsia* and *Francisella/Coxiella*-like endosymbionts in *Anatolia. International Journal of Infectious Diseases*, 19(1), 146.

Buysse, M., Floriano, A. M., Gottlieb, Y., Nardi, T., Comandatore, F., Olivieri, E., Giannetto, A., Palomar, A. M., Makepeace, B. L., Bazzocchi, C., Cafiso, A., Sassera, D., & Duron, O. (2021). A dual endosymbiosis supports nutritional adaptation to hematophagy in the invasive tick *Hyalomma marginatum*. *Life*, 10, e72747.

Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., & Holmes, S. P. (2016). DADA2: High-resolution sample inference from Illumina amplicon data. *Nature Methods*, 13(7), 581–583.

Clayton, K. A., Gall, C. A., Mason, K. L., Scoles, G. A., & Brayton, K. A. (2015). The characterization and manipulation of the bacterial microbiome of the Rocky Mountain wood tick, *Dermacentor andersoni*. *Parasites & Vectors*, 8, 632.

Cohen, R., Erez, K., ben-Avraham, D., & Havlin, S. (2000). Resilience of the internet to random breakdowns. *Physical Review Letters*, 85(21), 4626–4628.

Coyte, K. Z., Rao, C., Rakoff-Nahoum, S., & Foster, K. R. (2021). Ecological rules for the assembly of microbiome communities. *PLoS Biology*, 19(2), e3001116.

Dahmani, M., Davoust, B., Rousseau, F., Raoult, D., Fenollar, F., & Mediannikov, O. (2017). Natural Anaplasmataceae infection in *Rhipicephalus bursa* ticks collected from sheep in the French Basque Country. *Ticks and Tick-Borne Diseases*, 8(1), 18–24.

Debray, R., Herbert, R. A., Jaffe, A. L., Crits-Christoph, A., Power, M. E., & Koskella, B. (2022). Priority effects in microbiome assembly. *Nature Reviews Microbiology*, 20, 109–121.

DeSantis, T. Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E. L., Keller, K., Huber, T., Dalevi, D., Hu, P., & Andersen, G. L. (2006). Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. *Applied and Environmental Microbiology*, 72(7), 5069–5072.

Duron, O., Morel, O., Noël, V., Buysse, M., Binetruy, F., Lancelot, R., Loire, E., Ménard, C., Bouchez, O., Vavre, F., & Vial, L. (2018). Tick-bacteria mutualism depends on B vitamin synthesis pathways. *Current Biology*, 28(12), 1896–1902.e5.

Estrada-Peña, A., Bouattour, A., Camicas, J. L., & Walker, A. R. (2004). *Ticks of veterinary and medical importance: The Mediterranean basin. A guide of identification of species*. University of Zaragoza Press.

Estrada-Peña, A., Cabezas-Cruz, A., & Obregón, D. (2020). Resistance of tick gut microbiome to anti-tick vaccines, pathogen infection and antimicrobial peptides. *Pathogens*, 9(4), 309.

Fernandes, A. D., Macklaim, J. M., Linn, T. G., Reid, G., & Gloor, G. B. (2013). ANOVA-like differential expression (ALDEX) analysis for mixed population RNA-Seq. *PLoS One*, 8(7), e67019.

Fernandes, A. D., Reid, J. N. S., Macklaim, J. M., McMurrrough, T. A., Edgell, D. R., & Gloor, G. B. (2014). Unifying the analysis of high-throughput sequencing datasets: Characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. *Microbiome*, 2, 15.

Fernández, L., Pannaraj, P. S., Rautava, S., & Rodríguez, J. M. (2020). The microbiota of the human mammary ecosystem. *Frontiers in Cellular and Infection Microbiology*, 10, 586667.

Friedman, J., & Alm, E. J. (2012). Inferring correlation networks from genetic survey data. *PLoS Computational Biology*, 8(9), e1002687.

Gibbons, S. M. (2020). Keystone taxa indispensable for microbiome recovery. *Nature Microbiology*, 5(9), 1067–1068.

Grech-Angelini, S., Stachurski, F., Lancelot, R., Boissier, J., Allienne, F., Marco, S., Maestrini, O., & Uilenberg, G. (2016). Ticks (Acari: Ixodidae) infesting cattle and some other domestic and wild hosts on the French Mediterranean Island of Corsica. *Parasites & Vectors*, 9, 582.

Gryaznova, M. V., Syromyatnikov, M. Y., Dvoretskaya, Y. D., Solodskikh, S. A., Klimov, N. T., Mikhalev, V. I., Zimnikov, V. I., Mikhaylov, E. V., & Popov, V. N. (2021). Microbiota of cow's milk with udder pathologies. *Microorganisms*, 9(9), 1974.

Guseva, K., Darcy, S., Simon, E., Alteio, L. V., Montesinos-Navarro, A., & Kaiser, C. (2022). From diversity to complexity: Microbial networks in soils. *Soil Biology and Biochemistry*, 169, 108604.

Hamilton, P. T., Maluenda, E., Sarr, A., Belli, A., Hurry, G., Duron, O., Plantard, O., & Voordouw, M. J. (2021). *Borrelia afzelii* infection in the rodent host has dramatic effects on the bacterial microbiome of *Ixodes ricinus* ticks. *Applied and Environmental Microbiology*, 87(18), e0064121.

Hussain, S., Perveen, N., Hussain, A., Song, B., Aziz, M. U., Zeb, J., Li, J., George, D., Cabezas-Cruz, A., & Sparagano, O. (2022). The symbiotic continuum within ticks: Opportunities for disease control. *Frontiers in Microbiology*, 13, 854803.

Jakimovski, D., Mateska, S., Simin, V., Bogdan, I., Mijatović, D., Estrada-Peña, A., Mateos-Hernández, L., Foucault-Simonin, A., Moutailler,
S., Cabezas-Cruz, A., & Banović, P. (2022). Mediterranean spotted fever-like illness caused by *Rickettsia sibirica* mongolitimonae, North Macedonia, June 2022. *Euro Surveillance*, 27(42), 2200735.

Katoh, K., Misawa, K., Kuma, K., & Miyata, T. (2002). MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. *Nucleic Acids Research*, 30(14), 3059–3066.

Lejal, E., Chiquet, J., Aubert, J., Robin, S., Estrada-Peña, A., Rue, O., Midoux, C., Mariadassou, M., Bailly, X., Cougoul, A., Gasqui, P., Cosson, J. F., Chalvet-Monfray, K., Vayssier-Taussat, M., & Pollet, T. (2021). Temporal patterns in *Ixodes ricinus* microbial communities: An insight into tick-borne microbe interactions. *Microbiome*, 9(1), 153.

Lhomme, S. (2015). Spatial analysis of the structure of technical networks in a risk context. *Cyberge: European Journal of Geography* [Online], *Systems, Modeling, Geostatistics*, 711. https://doi.org/10.4000/cyberge.26763

Mahana, D., Trent, C. M., Kurtz, Z. D., Bokulich, N. A., Battaglia, T., Chung, J., Müller, C. L., Li, H., Bonneau, R. A., & Blaser, M. J. (2016). Antibiotic perturbation of the murine gut microbiome enhances the adiposity, insulin resistance, and liver disease associated with high-fat diet. *Genome Medicine*, 8(1), 48.

Maitre, A., Wu-Chuang, A., Aželytė, J., Palinauskas, V., Mateos-Hernández, L., Obregon, D., Hodžić, A., Moro, C. V., Estrada-Peña, A., Paoli, J.-C., Falchi, A., & Cabezas-Cruz, A. (2022). Vector microbiota manipulation by host antibodies: The forgotten strategy to develop transmission-blocking vaccines. *Parasites & Vectors*, 15(1), 4.

Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Foucault-Simonin, A., Moutailler, S., Paoli, J.-C., & Cabezas-Cruz, A. (2022). *Rickettsia helveticia* infection is associated with microbiome modulation in *Ixodes ricinus* collected from humans in Serbia. *Scientific Reports*, 12(1), 11464.

Martínez, I., Maldonado-Gomez, M. X., Gomes-Neto, J. C., Kittana, H., Ding, H., Schmaltz, R., Joglekar, P., Cardona, R. J., Marsteller, N. L., Kembel, S. W., Benson, A. K., Peterson, D. A., Ramer-Tait, A. E., & Walter, J. (2018). Experimental evaluation of the importance of colonization history in early-life gut microbiota assembly. *eLife*, 7, e36521.

Mateos-Hernández, L., Obregón, D., Maye, J., Borneres, J., Versille, N., de la Fuente, J., Estrada-Peña, A., Hodžić, A., Šimo, L., & Cabezas-Cruz, A. (2020). Anti-tick microbiota vaccine impacts *Ixodes ricinus* performance during feeding. *Vaccine*, 8(4), 702.

Mateos-Hernández, L., Obregón, D., Wu-Chuang, A., Maye, J., Bornères, J., Versillé, N., de la Fuente, J., Díaz-Sánchez, S., Bermúdez-Humarán, L. G., Torres-Maravilla, E., Estrada-Peña, A., Hodžić, A., Šimo, L., & Cabezas-Cruz, A. (2021). Anti-microbiota vaccines modulate the tick microbiome in a taxon-specific manner. *Frontiers in Immunology*, 12, 704621.

Mediannikov, O. Y., Sidelnikov, Y., Ivanov, L., Mokretsova, E., Fournier, P.-E., Tarasevich, I., & Raoult, D. (2004). Acute tick-borne rickettsiosis caused by *Rickettsia heilongjiangensis* in Russian Far East. *Emerging Infectious Diseases*, 10(5), 810–817.

Meyer, J. L., Gunasekera, S. P., Scott, R. M., Paul, V. J., & Teplitski, M. (2016). Microbiome shifts and the inhibition of quorum sensing by black band disease cyanobacteria. *The ISME Journal*, 10(5), 1204–1216.

Michelet, L., Delannoy, S., Devillers, E., Umhang, G., Aspan, A., Juremalm, M., Chirico, J., van der Wal, F. J., Sprong, H., Boye Pihl, T. P., Klitgaard, K., Bødker, R., Fach, P., & Moutailler, S. (2014). High-throughput screening of tick-borne pathogens in Europe. *Frontiers in Cellular and Infection Microbiology*, 4, 103.

Narasimhan, S., Rajeevan, N., Graham, M., Wu, M.-J., DePonte, K., Marion, S., Masson, O., O'Neal, A. J., Pedra, J. H. F., Sonenshine, D. E., & Fikrig, E. (2022). Tick transmission of *Borrelia burgdorferi* to the murine host is not influenced by environmentally acquired midgut microbiota. *Microbiome*, 10(1), 173.

Narasimhan, S., Rajeevan, N., Liu, L., Zhao, Y. O., Heisig, J., Pan, J., Eppler-Epstein, R., Deponte, K., Fish, D., & Fikrig, E. (2014). Gut

microbiota of the tick vector *Ixodes scapularis* modulate colonization of the Lyme disease spirochete. *Cell Host & Microbe*, 15(1), 58–71.

Nováková, M., & Šmajs, D. (2018). *Rickettsial endosymbionts of ticks*. Ticks and Tick-Borne Pathogens. IntechOpen.

Oksanen, J., Simpson, G. L., Blanchet, G., Kindt, R., Legendre, P., Minchin, P. R., O'Hara, R. B., Solymos, P., Stevens, H. M. H., Szöcs, E., Wagner, H. H., Barbour, M., Bedward, M., Bolker, B., Borcard, D., Carvalho, G., Chirico, M., De Cáceres, M., Durand, S., ... Weedom, J. (2021). *Vegan: community ecology package*. R package version 2.6-0.

Papa, A., Tsioka, K., Kontana, A., Papadopoulos, C., & Giadinis, N. (2017). Bacterial pathogens and endosymbionts in ticks. *Ticks and Tick-Borne Diseases*, 8(1), 31–35.

Peschel, S., Müller, C. L., von Mutius, E., Boulesteix, A.-L., & Depner, M. (2021). NetCoMi: Network construction and comparison for microbiome data in R. *Briefings in Bioinformatics*, 22(4), bbaa290.

Pielou, E. C. (1966). The measurement of diversity in different types of biological collections. *Journal of Theoretical Biology*, 13, 131–144.

Price, M. N., Dehal, P. S., & Arkin, A. P. (2010). FastTree 2—Approximately maximum-likelihood trees for large alignments. *PLoS One*, 5(3), e9490.

Raele, D. A., Galante, D., Pugliese, N., De Simone, E., & Cafiero, M. A. (2015). *Coxiella*-like endosymbiont associated to the «Anatolian brown tick» *Rhipicephalus bursa* in Southern Italy. *Microbes and Infection*, 17(11–12), 799–805.

Real, R., & Vargas, J. M. (1996). The probabilistic basis of Jaccard's index of similarity. *Systematic Biology*, 45(3), 380–385.

Ren, T., & Wu, M. (2016). PhyloCore: A phylogenetic approach to identifying core taxa in microbial communities. *Gene*, 593(2), 330–333.

Riera, J. L., & Baldo, L. (2020). Microbial co-occurrence networks of gut microbiota reveal community conservation and diet-associated shifts in cichlid fishes. *Animal Microbiome*, 2(1), 36.

Röttjers, L., & Faust, K. (2018). From hairballs to hypotheses-biological insights from microbial networks. *FEMS Microbiology Reviews*, 42(6), 761–780.

RStudio Team. (2020). RStudio: Integrated development for R. RStudio, PBC. http://www.rstudio.com/

Ruhnau, B. (2000). Eigenvector-centrality—A node-centrality? *Social Networks*, 22(4), 357–365.

Sperling, J. L. H., Fitzgerald, D., Sperling, F. A. H., & Magor, K. E. (2020). Microbiome composition and *Borrelia* detection in *Ixodes scapularis* ticks at the northwestern edge of their range. *Tropical Medicine and Infectious Disease*, 5(4), 173.

Thapa, S., Zhang, Y., & Allen, M. S. (2019). Effects of temperature on bacterial microbiome composition in *Ixodes scapularis* ticks. *MicrobiologyOpen*, 8(5), e00719.

Timsit, E., Workentine, M., van der Meer, F., & Alexander, T. (2018). Distinct bacterial metacommunities inhabit the upper and lower respiratory tracts of healthy feedlot cattle and those diagnosed with bronchopneumonia. *Veterinary Microbiology*, 221, 105–113.

Wallménius, K., Barboutis, C., Fransson, T., Jaenson, T. G. T., Lindgren, P.-E., Nyström, F., Olsen, B., Salaneck, E., & Nilsson, K. (2014). Spotted fever *Rickettsia* species in *Hyalomma* and *Ixodes* ticks infesting migratory birds in the European Mediterranean area. *Parasites & Vectors*, 7, 318.

Wang, S., Hua, X., & Cui, L. (2021). Characterization of microbiota diversity of engorged ticks collected from dogs in China. *Journal of Veterinary Science*, 22(3), e37.

Weese, J. S. (2013). The canine and feline skin microbiome in health and disease. *Veterinary Dermatology*, 24(1), 137–145.e31.

Wei, N., Cao, J., Zhang, H., Zhou, Y., & Zhou, J. (2021). The tick microbiota dysbiosis promote tick-borne pathogen transstadial transmission in a *Babesia microti*-infected mouse model. *Frontiers in Cellular and Infection Microbiology*, 11, 713466.

Werner, J. J., Koren, O., Hugenholtz, P., DeSantis, T. Z., Walters, W. A., Caporaso, J. G., Angenent, L. T., Knight, R., & Ley, R. E. (2012).
Impact of training sets on classification of high-throughput bacter-
rial 16s rRNA gene surveys. *The ISME Journal*, **6**(1), 94–103.

Wu-Chuang, A., Hodžić, A., Mateos-Hernández, L., Estrada-Peña, A., Obregon, D., & Cabezas-Cruz, A. (2021). Current debates and advances in tick microbiome research. *Current Research in Parasitology & Vector-Borne Diseases*, **1**, 100036.

Wu-Chuang, A., Obregon, D., Estrada-Peña, A., & Cabezas-Cruz, A. (2022). Thermostable keystone bacteria maintain the functional diversity of the *Ixodes scapularis* microbiome under heat stress. *Microbial Ecology*, *84*(4), 1224-1235. https://doi.org/10.1007/s00248-021-01929-y

Wu-Chuang, A., Obregon, D., Mateos-Hernández, L., & Cabezas-Cruz, A. (2022). Anti-tick microbiota vaccines: How can this actually work? *Biologia*, **77**, 1555-1562.

Zhong, J. (2012). *Coxiella*-like endosymbionts. *Advances in Experimental Medicine and Biology*, **984**, 365–379.

Zhong, J., Jasinskas, A., & Barbour, A. G. (2007). Antibiotic treatment of the tick vector *Amblyomma americanum* reduced reproductive fitness. *PLoS One*, **2**(5), e405.

SUPPORTING INFORMATION

Additional supporting information can be found online in the
Supporting Information section at the end of this article.

**How to cite this article:** Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Piloto-Sardiñas, E., Foucault-Simonin, A., Cicculli, V., Moutailler, S., Paoli, J.-C., Falchi, A., Obregón, D., & Cabezas-Cruz, A. (2023). Rickettsial pathogens drive microbiota assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks. *Molecular Ecology*, *00*, 1–17. https://doi.org/10.1111/mec.17058