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

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 amplicon 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.

While the distinction between human pathogens and other nonpathogenic 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, 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; 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. 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; 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;. 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 . In contrast to pathogen infection, anti-tick immunity significantly reduced the network robustness . 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;.
This suggests that TBP infection induces minimal changes on the tick bacterial microbiota robustness , and that once established; pathogens occupy marginal positions in the bacterial community (Maitre, Wu-Chuang, Aželytė, 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'.
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.

| 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

| 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.

| 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 putted to 0 in the ASVs table, and the following tests were realised with the same procedure.

| 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).

| Inference of bacterial co-occurrence networks
Here we used co-occurrence network to compared 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.

| 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).

| 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).

| 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: (i) ubiquitousness, (ii) high eigenvector centrality (≥0.75), (iii) high relative abundance (clr value above the average).

| 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.

| 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). Ehrlichia minasensis (H. marginatum, 4.3%, total 23) were found at lower frequencies.

| Rickettsial pathogens infection and microbiota diversity in H. marginatum and R. bursa
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

| 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.  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 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) (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.

| 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). F I G U R E 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  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 ad-  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.

| Rickettsia as a driver of microbial community assembly in H. marginatum and R. bursa
We conducted an evaluation to assess the potential influence of in HmM2) were gained by the modules (Figure 4b). Similarly, remov-  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.

| DISCUSS ION
The capacity of TBPs to modulate the diversity, bacterial composition and microbial structure of tick host microbiota has been  established (Abraham et al., 2017;Estrada-Peña et al., 2020;Hamilton et al., 2021;Maitre, Wu-Chuang, Aželytė, 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 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 re- 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;. 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;. 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, 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 microbemicrobe 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 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 (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, 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 loca- tions. 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 summatory 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 pathogenassociated 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.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasets generated and analysed during the current study are available on the SRA repository (Bioproject No. PRJNA865094).

B EN EFIT-S H A R I N G S TATEM ENT
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.