# Exploring the impact of Anaplasma phagocytophilum on colonization resistance of Ixodes scapularis microbiota using network node manipulation

Lianet Abuin-Denis, Elianne Piloto-Sardiñas, Apolline Maitre, Alejandra

Wu-Chuang, Lourdes Mateos-Hernández, Dasiel Obregon, Belkis

Corona-González, Andréa Cristina Fogaça, Vaidas Palinauskas, Justė Aželytė, et al.

► To cite this version:

Lianet Abuin-Denis, Elianne Piloto-Sardiñas, Apolline Maitre, Alejandra Wu-Chuang, Lourdes Mateos-Hernández, et al. Exploring the impact of Anaplasma phagocytophilum on colonization resistance of Ixodes scapularis microbiota using network node manipulation. Current Research in Parasitology & Vector-Borne Diseases, 2024, 5, pp.100177. 10.1016/j.crpvbd.2024.100177. hal-04635884

HAL Id: hal-04635884

https://hal.science/hal-04635884v1

Submitted on 4 Jul 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Check for
updates

# Exploring the impact of *Anaplasma phagocytophilum* on colonization resistance of *Ixodes scapularis* microbiota using network node manipulation

Lianet Abuin-Denis<sup>a, b</sup>, Elianne Piloto-Sardiñas<sup>b, c</sup>, Apolline Maître<sup>b, d, e</sup>, Alejandra Wu-Chuang<sup>b</sup>, Lourdes Mateos-Hernández<sup>b</sup>, Dasiel Obregon<sup>f</sup>, Belkis Corona-González<sup>c</sup>, Andréa Cristina Fogaça<sup>g</sup>, Vaidas Palinauskas<sup>h</sup>, Justè Aželytè<sup>h</sup>, Alina Rodríguez-Mallon<sup>a</sup>, Alejandro Cabezas-Cruz<sup>b, *</sup>

## ARTICLE INFO

### Keywords:

Ticks

Tick-borne pathogens

Community assembly

Colonization resistance

## ABSTRACT

Upon ingestion from an infected host, tick-borne pathogens (TBPs) have to overcome colonization resistance, a defense mechanism by which tick microbiota prevent microbial invasions. Previous studies have shown that the pathogen *Anaplasma phagocytophilum* alters the microbiota composition of the nymphs of *Ixodes scapularis*, but its impact on tick colonization resistance remains unclear. We analyzed tick microbiome genetic data using published Illumina 16S rRNA sequences, assessing microbial diversity within ticks (alpha diversity) through species richness, evenness, and phylogenetic diversity. We compared microbial communities in ticks with and without infection with *A. phagocytophilum* (beta diversity) using the Bray-Curtis index. We also built co-occurrence networks and used node manipulation to study the impact of *A. phagocytophilum* on microbial assembly and network robustness, crucial for colonization resistance. We examined network robustness by altering its connectivity, observing changes in the largest connected component (LCC) and the average path length (APL). Our findings revealed that infection with *A. phagocytophilum* does not significantly alter the overall microbial diversity in ticks. Despite a decrease in the number of nodes and connections within the microbial networks of infected ticks, certain core microbes remained consistently interconnected, suggesting a functional role. The network of infected ticks showed a heightened vulnerability to node removal, with smaller LCC and longer APL, indicating reduced resilience compared to the network of uninfected ticks. Interestingly, adding nodes to the network of infected ticks led to an increase in LCC and a decrease in APL, suggesting a recovery in network robustness, a trend not observed in networks of uninfected ticks. This improvement in network robustness upon node addition hints that infection with *A. phagocytophilum* might lower ticks' resistance to colonization, potentially facilitating further microbial invasions. We conclude that the compromised colonization resistance observed in tick microbiota following infection with *A. phagocytophilum* may facilitate co-infection in natural tick populations.

<footnote><sup>a</sup> Animal Biotechnology Department, Center for Genetic Engineering and Biotechnology, Avenue 31 between 158 and 190, P.O. Box 6162, Havana, 10600, Cuba</footnote>

<footnote><sup>b</sup> ANSES, INRAE, Ecole Nationale Vétérinaire d'Alfort, UMR BIPAR, Laboratoire de Santé Animale, Maisons-Alfort, F-94700, France</footnote>

<footnote><sup>c</sup> Direction of Animal Health, National Center for Animal and Plant Health, Carretera de Tapaste y Autopista Nacional, Apartado Postal 10, San José de las Lajas, Mayabeque, 32700, Cuba</footnote>

<footnote><sup>d</sup> INRAE, UR 0045 Laboratoire de Recherches sur le Développement de l'Elevage (SELMET-LRDE), 20250, Corte, France</footnote>

<footnote><sup>e</sup> EA 7310, Laboratoire de Virologie, Université de Corse, Corte, France</footnote>

<footnote><sup>f</sup> School of Environmental Sciences, University of Guelph, Guelph, ON, Canada</footnote>

<footnote><sup>g</sup> Department of Parasitology, Institute of Biomedical Sciences, University of São Paulo, São Paulo, 05508-000, SP, Brazil</footnote>

<footnote><sup>h</sup> Nature Research Centre, Akademijos 2, Vilnius, Lithuania</footnote>

<footnote>* Corresponding author.
E-mail address: alejandro.cabezas@vet-alfort.fr (A. Cabezas-Cruz).</footnote>
# 1. Introduction

Ticks are obligate hematophagous ectoparasites that transmit infectious agents, including bacteria (such as *Coxiella burnetti*, *Rickettsia helvetica*, *Borrelia burgdorferi*, and *Anaplasma marginale*), viruses (like severe fever with thrombocytopenia syndrome virus, tick-borne encephalitis virus, and Crimean-Congo haemorrhagic fever virus), and parasites (such as *Babesia microti* and *Theileria orientalis*) to terrestrial vertebrates (de la Fuente et al., 2017). Once acquired from an infected host, these tick-borne pathogens (TBPs) encounter several barriers, such as the peritrophic membrane, the dityrosine network, and tick immunity (Kurokawa et al., 2020; Fogaça et al., 2021; Kitsou et al., 2021). In addition to pathogens, ticks also carry non-pathogenic microorganisms (Bonnet et al., 2017), which include commensal microbes acquired from the environment, as well as transovarially-transmitted endosymbionts (Binetruy et al., 2020; Hussain et al., 2022), collectively referred to as microbiota (sous rature) (Cabezas-Cruz, 2023).

These non-pathogenic microorganisms gradually assemble over time, progressing from a state of low diversity to form richer multic species communities that can have a substantial effect on the structure, organization, and function of the tick microbiota (Díaz-Sánchez et al., 2019; Obregón et al., 2019; Estrada-Pena et al., 2020a). In addition, non-pathogenic microorganisms may also play a role in driving the transmission of TBPs, which has significant implications for both human and animal health (Narasimhan et al., 2014; Gall et al., 2016; Wei et al., 2021; Wu-Chuang et al., 2021). There is empirical evidence for the presence of colonization resistance within the tick microbiota. Colonization resistance is the phenomenon where established microbial communities prevent the invasion and establishment of new, often pathogenic, species (Mullineaux-Sanders et al., 2018; Ducarmon et al., 2019; Stacy et al., 2021; Karita et al., 2022). For example, in *Derma-centor andersoni* ticks, the resident microbiota influences the acquisition levels of certain pathogens (Gall et al., 2016). Particularly, *Rickettsia bellii*, a non-pathogenic microorganism associated with *D. andersoni*, negatively correlates with *Anaplasma marginale* acquisition, indicating an antagonistic interaction (Gall et al., 2016). Additionally, experiments with *Haemaphysalis longicornis* ticks show that adult ticks emerging from nymphs treated with antibiotics exhibit a disrupted microbiota (Wei et al., 2021). These adult ticks had a higher infection rate of *Babesia microti* (44.7%) compared to control ticks (24.2%) (Wei et al., 2021), indicating that a healthy microbiota may play a crucial role in managing pathogen loads. This highlights that colonization resistance serves as a natural barrier against the establishment of TBPs in ticks.

Tick-transmitted bacterial pathogens likely evolved mechanisms to overcome colonization resistance by resident microbiota. For example, *Borrelia burgdorferi* colonization in *I. scapularis* increases the expression of several tick gut genes including *pixr*. Abrogation of *PIXR* function results in alterations in the gut microbiome, metabolome, and immune responses affecting the spirochete entering the tick gut (Narasimhan et al., 2014). Other mechanisms might involve triggering the expression of tick proteins with anti-microbial activity, such as *Ixodes scapularis* antifreeze glycoprotein (IAFGP), perturbing the tick gut microbiota (Abraham et al., 2017). This alteration affects the capacity of bacteria to form biofilms, influences the integrity of the peritrophic matrix, and reduces barriers to *Anaplasma phagocytophilum* colonization in the tick (Abraham et al., 2017). Overall, these mechanisms show an adaptation by TBPs to manipulate the tick microbial environment, effectively reducing colonization resistance and allowing for their proliferation and transmission. However, it remains unclear whether TBP-mediated modulation of tick microbiota results in community traits associated with reduced colonization resistance.

Altered microbial interactions and a weaker network structure have been proposed to lower colonization resistance (Shade et al., 2012). Disruptions in how these microbial communities are structured may reduce their overall stability and resilience, possibly affecting their capacity to fend off pathogens (Shade et al., 2012). A recent study by

Maitre et al. (2022) found that an *R. helvetica* infection in the tick *Ixodes ricinus* significantly diminishes the diversity and connectivity within the tick's microbiota network, indicating a decrease in colonization resistance due to the infection Cohen and Havlin (2009).

In this study, we investigated colonization resistance of *I. scapularis* nymph microbiota infected with *A. phagocytophilum*. *Anaplasma phagocytophilum* is an obligate intracellular bacterium transmitted by ticks, causing human granulocytic anaplasmosis (Kocan et al., 2015). We hypothesized that alteration of the tick microbiota by *A. phagocytophilum* results in reduced colonization resistance, and diminished ability of the resident microbiota in the tick gut to prevent colonization by other microorganisms, including pathogens. To assess this hypothesis, we used published 16S rRNA amplicon sequencing data (Abraham et al., 2017) to compare the community assembly and robustness of *A. phagocytophilum*-infected and uninfected ticks, using a network approach. Co-occurrence networks (Estrada-Pena et al., 2020b; Mateos-Hernández et al., 2020, 2021, 2023; Maitre et al., 2023), with nodes representing individual microbial taxa and edges representing their interactions (Faust and Raes, 2012; Röttjers and Faust, 2018), were used to assess the impact of *A. phagocytophilum* on the tick microbial communities. Our findings indicate that infection with *A. phagocytophilum* reduced network robustness, potentially compromising colonization resistance. Interestingly, we observed a robustness recovery with the addition of new nodes in the *A. phagocytophilum*-infected network, suggesting that after infection, further microbial invasions may recover colonization resistance. These results are relevant for understanding the microbiota dynamics and responses to pathogen infections in ticks.

# 2. Materials and methods

## 2.1. Original datasets

We evaluated the effect of bacterial infection on tick microbiota using publicly available datasets (Fig. 1). The selected studies utilized barcoded universal primers to amplify the V4 hypervariable regions of the 16S rRNA gene, followed by sequencing using the Illumina MiSeq system. The raw data were obtained from Abraham et al. (2017). In their study, the authors examined the changes in gut microbiota composition of *I. scapularis* nymphs fed on mice (C3H/SCID), which were experimentally infected with *A. phagocytophilum* (strain NCH-1). DNA was extracted from the guts of individual fed nymphs. The V4 variable region of the bacterial 16S rRNA was amplified from the genomic DNA of each sample using 12-base barcoded primer sets. In our study, we referred to *A. phagocytophilum* (Ap) tick groups as Ap-infected (n = 22) and Ap-uninfected (n = 10).

## 2.2. Processing of original raw sequences

We conducted taxonomic profiling of the 16S rRNA gene sequence datasets from the previously referenced studies. The sequences were sourced from the SRA repository: PRJNA353730. For sequence data processing, we used the pipelines in Quantitative Insights Into Microbial Ecology 2 Software (QIIME2) (version 2021.4.0) (Bolyen et al., 2019), as described in our previous publications (Estrada-Pena et al., 2020b; Wu-Chuang et al., 2023). In summary, the demultiplexed fastq files were denoised utilizing the DADA2 algorithm (version 1.18.0) as described by Callahan et al. (2016). Following this, all amplicon sequence variants (ASVs) were aligned employing the MAFFT tool (Katoh, 2002) through the q2-alignment interface. A phylogenetic tree was then generated using FastTree (Price et al., 2010) via the q2-phylogeny plugin. Diversity metrics, both phylogenetic and non-phylogenetic, were calculated based on a rarefied table of ASV abundances using the core-metrics-phylogenetic function within the q2-diversity plugin. Taxonomic annotation of the ASVs was carried out using the q2-feature-classifier's classify-sklearn naive Bayes taxonomy classifier,
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Experimental design. The data used in this study were downloaded from the SRA database and the raw sequences were processed. Bacterial composition and abundance were analyzed. The alpha- and beta diversity of the two datasets were compared, and co-occurrence networks were inferred to assess the structure of the microbial communities.</figcaption></figure>

referencing the 16S rRNA SILVA database (release 132) (Quast et al., 2013) for alignment.

## 2.3. Diversity indices

To test for differences in bacterial diversity between Ap-uninfected and Ap-infected samples, we conducted analyses of alpha and beta diversity. For alpha diversity, we used the metrics “observed features” (DeSantis et al., 2006) and Faith’s phylogenetic diversity index (Faith, 1992). Evenness was assessed using Pielou’s evenness index (Pielou, 1966). The observed features metrics is the number of unique features present in a sample (DeSantis et al., 2006) and Faith’s phylogenetic diversity index measures the cumulative evolutionary history represented in the community (Faith, 1992). Pielou’s evenness index measures how evenly individuals are distributed among species within a community, indicating the balance of species abundances relative to species richness (Pielou, 1966). On the other hand, beta diversity measures the diversity between samples and examines the similarity in composition between analyzed communities. Beta diversity was evaluated using the principal coordinates analysis (PCoA), based on the Bray-Curtis dissimilarity index (Bray and Curtis, 1957). Unique and shared taxa among the conditions were visualized using Venn diagrams created with an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) and the *Upset* package implemented in R v.4.3.1 (R Core Team, 2023). Analyses were performed using the RStudio Integrated Development Environment (IDE) v.2023.03.0-daily+82. pro2 (RStudio Team, 2020).

## 2.4. Inference of bacterial co-occurrence networks

Co-occurrence networks were employed to analyze bacterial interaction within microbial communities, both in Ap-uninfected and Ap-

infected groups. These networks visually depict the relationships among interacting microbes, where an edge connects two bacteria if their abundances demonstrate significant correlation across samples under the same condition or treatment. The networks were constructed with taxonomic data tables at both family and genus levels (Supplementary file S2: Table S1), using the Sparse Correlations for Compositional data (SparCC) algorithm (Friedman and Alm, 2012), as implemented (Supplementary file S1) in the SpiecEasi R package (Kurtz et al., 2015). Only significant correlations between taxa, both negative and positive (SparCC, weight ≥ 0.5 or ≤ -0.5) were represented as edges. Furthermore, we also analyzed the strongest correlations (SparCC ≥ 0.9 or ≤ -0.9). Network metrics used in this study (Table 1) were computed using Gephi 0.9.5 (Bastian et al., 2009), an open-source software that transforms co-occurrence data in a graph.

## 2.5. Network comparisons

To compare networks, a statistical estimation analysis was conducted using Network Construction and Comparison with Microbes (*NetCoMi*) package (Peschel et al., 2021) in R v.4.3.1 (R Core Team, 2023) (Supplementary file S1), and performed using RStudio (RStudio Team, 2020). *NetCoMi* offers tools for networks alignment, which involves matching nodes (microbial taxa) and edges (co-occurrence relationships) between networks based on their topological properties helping to identify corresponding features between them, even if they are not identical. To assess similarities in the distribution of local centrality measures across nodes, i.e. degree-, betweenness-, closeness-, and eigenvector centrality (Table 1), between two networks, we computed the Jaccard index for each centrality measure. The Jaccard index measures the similarity between sets of “most central nodes”, i.e. nodes with a centrality value above the empirical 75th quartile, in the two
<table><caption>**Table 1**
Key network metrics and their ecological interpretations in microbial community analysis.</caption><thead><tr><td>Measure</td><td>Definition</td><td>Ecological interpretation</td></tr></thead><tbody><tr><td>Network diameter</td><td>The network diameter measures the shortest path between the two most distant nodes in a network. It provides insight into the maximum number of edges required to connect any two nodes.</td><td>In microbial networks, a large diameter might suggest that some microbial interactions are spread over larger "distances", implying less direct interaction between certain microbe taxa, potentially due to ecological or functional divergence.</td></tr><tr><td>Modularity</td><td>Modularity assesses the strength of network division into modules or communities. It quantifies the extent to which nodes in the same module are more densely connected to each other compared to nodes in different modules.</td><td>High modularity in a microbial network suggests that the network is divided into distinct communities or clusters, each possibly representing different ecological niches or functional groups within the microbiome.</td></tr><tr><td>Average degree</td><td>The average degree of a network determines the average number of edges per node. It reflects the overall connectivity of the network by measuring how many connections each node has on average.</td><td>In a microbial context, a higher average degree suggests a community where numerous microbe taxa interact with several others, which might imply a robust ecological network where many species are involved in maintaining community structure and function.</td></tr><tr><td>Weighted degree</td><td>The weighted degree takes into account the sum of edge weights connected to a node. It considers the strength or correlation intensity associated with each edge, providing a more nuanced view of node connectivity.</td><td>Considering the intensity or strength of connections, this metric provides insights into which microbe taxa play central roles based on the strength of their interactions.</td></tr><tr><td>Clustering coefficient</td><td>The clustering coefficient captures the tendency of nodes to form clusters or tightly interconnected groups. It measures the extent to which a node's neighbors are connected, indicating the presence of local clustering or community structures within the network.</td><td>A high clustering coefficient in microbial networks indicates a tendency of microbe taxa to form tightly knit groups, suggesting the presence of cooperative clusters or consortia.</td></tr><tr><td>Number of communities</td><td>Partitioning of nodes (representing entities or elements) within a network into distinct groups or communities based on their structural connectivity.</td><td>In microbial networks, each community might represent a different ecological niche or a group of microbe taxa performing similar functions.</td></tr><tr><td>Number of triangles</td><td>The number of triangles that exist within a given network. A triangle in a network is a set of three nodes that are all connected.</td><td>In microbial communities, triangles (three microbe taxa all interacting with each other) might suggest robust sub-communities that could stabilize the network against perturbations by providing redundant paths for interaction.</td></tr><tr><td>Largest Connected Component (LCC)</td><td>The LCC represents the main connected structure of the network. It identifies the largest subset of nodes that are mutually reachable through edges.</td><td>The size of the LCC in a microbial network indicates the core structure of connectivity and can reflect the main functional and structural backbone of the community.</td></tr><tr><td>Average Path Length (APL)</td><td>The average path length measures the efficiency of information flow within the network. It calculates the average number of steps required to travel between any two nodes in the network, indicating how quickly information can spread through the network.</td><td>In microbiological networks, the average path length relates to the efficiency of material or signal transfer across the network.</td></tr><tr><td>Degree centrality</td><td>Degree centrality measures the number of edges connected to a node, indicating the importance or influence of a node based on the number of connections it has.</td><td>In microbial networks, a high degree centrality indicates a species with numerous interactions, which could suggest a generalist species that engages with many different partners or a keystone species that plays a critical role in maintaining community structure and stability.</td></tr><tr><td>Betweenness centrality</td><td>Betweenness centrality quantifies the extent to which a node lies on the shortest paths between other nodes. It identifies nodes that act as intermediaries or bridges, facilitating communication in the network.</td><td>In ecological terms, species with high betweenness centrality may be those that link otherwise disparate groups of organisms, facilitating important ecological processes such as energy or material transfer across the community.</td></tr><tr><td>Closeness centrality</td><td>Closeness centrality measures how close a node is to all other nodes in the network. It reflects the efficiency of information flow from a node to other nodes, considering the shortest path lengths.</td><td>Species with high closeness centrality can quickly spread effects (either beneficial such as nutrients or detrimental such as pathogens) throughout the network, indicating their efficiency in influencing the community dynamics.</td></tr><tr><td>Eigenvector centrality</td><td>Eigenvector centrality considers both the local and global importance of a node. It assigns a centrality score to a node based on the centrality of its neighboring nodes, indicating its overall influence.</td><td>In a microbial community, such a metric would highlight species that are not only well-connected but also connected to other significant species, reinforcing their role in maintaining or disrupting complex community structures.</td></tr><tr><td>Hub taxa</td><td>Hub taxa are nodes in a network that exhibit high connectivity or act as hubs, having numerous connections with other nodes. They play a crucial role in maintaining network structure and information flow.</td><td>These taxa may represent species that provide essential ecosystem services, such as keystone species in ecological networks or core microbiota in host-associated microbial communities, whose loss might lead to drastic changes in network structure and function.</td></tr><tr><td>Core Association Network (CAN)</td><td>This network model identifies and visualizes the core set of interactions or associations between nodes (taxa) that are consistently present across multiple samples or conditions.</td><td>The CAN is useful for highlighting interactions that are critical to community structure and function, regardless of external conditions or perturbations.</td></tr></tbody></table>

networks. It expresses the similarity of the sets of most central nodes and the sets of hub taxa (highly connected nodes) between the two networks.

The Adjusted Rand Index (ARI) was also calculated to test the dissimilarity of clustering in the networks. Negative and positive ARI values range from -1 to 1, where values below 0 indicate lower than random clustering, values above 0 indicate higher than random clustering, a value of 1 corresponds to identical clustering, and a value of 0 indicates dissimilar clustering (Peschel et al., 2021).

Additionally, the Core Association Network (CAN) analysis (Röttjers et al., 2021), was performed using the Anuran software, implemented in Python environment (https://github.com/ramellose/anuran). The CAN analysis identifies conserved patterns across networks (Table 1), using a core prevalence threshold of 0.8. This approach utilizes null models to generate random networks and assesses the properties of these networks, allowing the identification of patterns in groups of networks. Our hypothesis was that CAN does not differentiate the microbial networks

of Ap-uninfected and Ap-infected *I. scapularis* nymphs, indicating that infection does not lead to altered core associations in the microbial network. CAN visualization was carried out using Gephi 0.9.5 (Bastian et al., 2009).

## 2.6. Network robustness analysis

We evaluated the robustness of microbial co-occurrence networks to perturbation by measuring the impact of removing or adding nodes on network connectivity. To assess this, we simulated the proportion of node removal required to reach a loss in connectivity of 0.80 for each network using random or directed attacks. For the directed attack, we employed three strategies: betweenness centrality, degree centrality, and cascading. In the betweenness centrality approach, we removed nodes with the highest betweenness centrality values. In the degree centrality approach, we removed nodes with the highest degree
centrality values. In the cascading approach, we first removed nodes with the highest betweenness centrality values, recalculating betweenness centrality after each node removal. To perform the network robustness analysis, we utilized the Network Strengths and Weaknesses Analysis (NetSwan) package (Lhomme, 2015) in R v.4.3.1 (R Core Team, 2023) (Supplementary file S1), performed using the RStudio (RStudio Team, 2020).

Additionally, we conducted a node addition analysis based on the method described by Freitas et al. (2020) in R v.4.3.1 (R Core Team, 2023) (Supplementary file S1), and performed using the RStudio (RStudio Team, 2020). In this analysis, new nodes were randomly selected and connected to the existing network. We then calculated the size of the largest connected component (LCC) and the average path length (APL) (Table 1). To obtain a more accurate estimate of the network's robustness, we repeated the simulation with different sets of nodes, adding 100, 300, 500, 700, and 1000 nodes. The obtained values were plotted using GraphPad Prism 8.0.1 to visualize the results.

## 2.7. Statistical analysis

Differences in alpha diversity between groups were assessed using the Kruskal-Wallis test (*P* < 0.05) in QIIME2 (Bolyen et al., 2019). The Bray-Curtis dissimilarity index values were compared between groups using a PERMANOVA test (*P* < 0.05). A PERMANOVA test (*P* < 0.05) was performed to analyze beta dispersion. Additionally, we calculated beta dispersion (variance) using the betadisper function from the *Vegan* package in R v.4.3.1 (R Core Team, 2023) (Supplementary file S1), and performed using RStudio (RStudio Team, 2020).

The differences in bacterial taxa abundance between the Ap-uninfected and Ap-infected groups were performed with the ANOVA-Like Differential Expression (ALDEX2) method (Fernandes et al., 2013) implemented in R v.4.3.1 (R Core Team, 2023), and performed using RStudio (RStudio Team, 2020). Relative abundance was measured as centered 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., 2013). The resulting data were used to construct the heatmap with the heatmap.2 function (Supplementary file S1), implemented in R v.4.3.1 (R Core Team, 2023), and performed using RStudio (RStudio Team, 2020). The comparisons were performed with Welch's *t*-test (*P* ≤ 0.05).

To test the similarity in the distribution of local centrality measures between two networks, two *P*-values, namely *P* (J ≤ j) and *P* (J ≥ j), were computed for each local centrality measure (Real and Vargas, 1996). These *P*-values represent the probability that the observed Jaccard index (J) value is either “less than or equal to” or “greater than or equal to” the expected Jaccard value at random (j). Differences were considered significant when *P* < 0.05.

The standard error for loss of connectivity was calculated, considering variability, using a threshold of 0.975. Additionally, the node addition analysis employed Wilcoxon signed-rank tests to determine if the mean size of the LCC and the APL differed significantly from 0. The *P*-values from these tests were adjusted using the Benjamini-Hochberg (BH) procedure to control for multiple comparisons. Additionally, bootstrapping was performed to obtain confidence intervals for the variables. Significance was determined at a threshold of *P* < 0.05.

# 3. Results

## 3.1. Impact of *A. phagocytophilum* on bacterial diversity and composition of tick microbiota

The results indicated that there were no significant differences in alpha diversity (Fig. 2A) and beta diversity (Fig. 2B) between the Ap-uninfected and Ap-infected conditions (Kruskal-Wallis, *P* > 0.05; PERMANOVA, *P* > 0.05). In terms of microbiota composition, a total of 409

bacterial taxa were identified (Supplementary file S2: Table S2). Of these, 4 bacterial taxa (0.98%) were exclusive to the Ap-uninfected group, 29 bacterial taxa (7.10%) were exclusive to the Ap-infected group, and 376 bacterial taxa (91.9%) were shared between the two groups (Fig. 2C). Next, we conducted a differential abundance analysis to identify changes in specific taxa between Ap-uninfected and Ap-infected groups. We observed a significantly different abundance of three taxa, *Lactococcus*, *Pantoea* and *Mycobacterium*, across the two conditions (Fig. 2D).

## 3.2. Impact of *A. phagocytophilum* on bacterial community assembly

We evaluated the impact of *A. phagocytophilum* infection on tick bacterial communities using co-occurrence networks. Visual examination of the networks (SparCC, weight ≥ 0.5 or ≤ -0.5) of Ap-uninfected (Fig. 3A) and Ap-infected (Fig. 3B) revealed topological differences. Out of 409 bacterial taxa that were found in the microbiota, there were 377 (92.2%) nodes in the Ap-uninfected network and 265 (64.8%) nodes were present in the Ap-infected network (Table 2). Analysis of the topological features of the networks revealed a higher number of nodes and edges in the Ap-uninfected network (Fig. 3A, Table 2) than in the Ap-infected network (Fig. 3B, Table 2). Ap-infected network presented a reduction of 30% (265 out of 377) of the number of nodes and 82% (826 out of 4605) of the number of edges compared with Ap-uninfected network (Table 2). Furthermore, in both network conditions, positive interactions between the nodes were predominant (Table 2). Upon infection, the microbial assembly underwent reorganization, resulting in a decrease in the modularity and an increase in the number of distinct communities (Table 2). There were 243 taxa (61%, total 399) shared between the two networks, while 134 (33.6%) and 22 (5.4%) were unique to the Ap-uninfected and Ap-infected networks, respectively (Fig. 3C, Table S3).

The Jaccard index was used to evaluate the similarities in selected local network centrality measures between Ap-uninfected and Ap-infected networks. The results showed that the Jaccard index of all local centrality measures was significantly lower (*P* ≤ Jacc) than expected by random for the comparisons of Ap-uninfected vs Ap-infected networks (Table 3). The ARI similarity index was employed to assess clustering differences between networks. The ARI supported the low similarity in clustering between Ap-infected and Ap-uninfected networks (ARI = -0.001, *P* = 0.891). The Core Association Network (CAN) revealed 40 core-interacting nodes (Fig. 3D, Supplementary file S2: Table S4), which represents a reduction of 89.4% and 85.0% in the number of nodes in Ap-uninfected and Ap-infected networks, respectively. This supported that infection causes a major impact in the tick microbial network, but also that a core set of interactions critical to community structure and function remains regardless of perturbation. Among these core-associated nodes, only positive interactions were found.

Finally, we examined the presence of strong correlations between network nodes. In the Ap-uninfected network, we identified 34 nodes with the strongest connections (Fig. 3E), and in the Ap-infected network 18 nodes were identified as the strongest connected (Fig. 3F). In the Ap-infected network, the nodes formed a major cluster, but this pattern was not observed in the Ap-uninfected network. Overall, these findings suggest significant alterations in the assembly of tick bacterial communities due to the presence of pathogens, which is supported by the observed topological differences, dissimilarity in local network centrality measures, and network clustering.

## 3.3. Impact of *A. phagocytophilum* on the network robustness

One crucial aspect of networks is their ability to withstand perturbations, including the removal or addition of nodes. In terms of connectivity loss, the removal of nodes through direct attack had a more significant impact on all networks compared to random attacks
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Impact of *A. phagocytophilum* infection on microbiome diversity. A Comparison of the alpha diversity by observed features (Kruskal-Wallis test, *P* = 0.31), Pielou's evenness index (Kruskal-Wallis test, *P* = 0.90) and Faith's index (Kruskal-Wallis test, *P* = 0.46). B Comparison of beta diversity for Ap-uninfected vs Ap-infected using the Bray-Curtis dissimilarity index (PERMANOVA, *F* = 32.84, *P* = 0.087, stress = 0.1507). Small circles and triangles in the principal coordinate analysis (PCoA) plot represent samples; ellipses indicate 95% confidence intervals. ANOVA test was performed and showed that the beta dispersion of the samples is not significantly different (*P* = 0.849). C Venn diagram showing the number of unique or shared taxa between Ap-uninfected and Ap-infected ticks. D Dendrogram heatmap resulting from the heatmap.2 functions implemented in R (R Core Team, 2023), and performed in RStudio IDE (RStudio Team, 2020). The taxa were clustered based on relative abundance (calculated as clr-transformed values). Each column represents the clr-values for bacterial taxa per sample and group. Lines represent bacterial taxa with significant changes between the two datasets. Color represents the clr-value.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Networks representing community assemblies. Co-occurrences networks of Ap-uninfected (A) and Ap-infected (B). Nodes correspond to taxa (family or genus level), and connecting edges indicate a correlation between them. Only nodes with at least one significant correlation are represented. Node colors are based on modularity class metric and equal color means modules of co-occurring taxa. Node sizes are based on eigenvector centrality. Edges represent positive (blue) or negative (red) correlations (SparCC ≥ 0.5 or ≤ -0.5). C Number of unique or shared taxa in the Ap-uninfected and Ap-infected network. D Core Association Network (CAN) between Ap-uninfected and Ap-infected groups. E, F Strong correlation networks (weight > 0.9 or < -0.9) across the Ap-uninfected (E) and Ap-infected networks (F). Nodes correspond to taxa (family or genus level) and connecting edges indicate significant correlation between them. Edges represent positive (green) correlations; edge widths are proportional to correlation coefficients.</figcaption></figure>
<caption>**Table 2**</caption>

Topological features of the taxonomic networks of Ap-uninfected and Ap-infected groups.

<table><caption>**Table 3**</caption><thead><tr><th rowspan="2">Network features</th><th colspan="2">Anaplasma phagocytophilum</th></tr><tr><th>Ap-uninfected</th><th>Ap-infected</th></tr></thead><tbody><tr><td>No. of nodes</td><td>377</td><td>265</td></tr><tr><td>No. of edges</td><td>4605</td><td>826</td></tr><tr><td>Positive interactions</td><td>2844 (61.75%)</td><td>764 (92.49%)</td></tr><tr><td>Negative interactions</td><td>1762 (38.26%)</td><td>62 (7.50%)</td></tr><tr><td>Modularity</td><td>1.753</td><td>0.4</td></tr><tr><td>Number of communities</td><td>6</td><td>62</td></tr><tr><td>Network diameter</td><td>6</td><td>18</td></tr><tr><td>Average degree</td><td>24.43</td><td>4.089</td></tr><tr><td>Weighted degree</td><td>4.038</td><td>2.305</td></tr><tr><td>Clustering coefficient (Triangles method)</td><td>0.436</td><td>0.699</td></tr><tr><td>No. of triangles</td><td>20,375</td><td>3910</td></tr><tr><td>Connectivity</td><td>1</td><td>24</td></tr></tbody></table>

Jaccard index for infected and uninfected tick networks.

<table><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">Jaccard index and statistical significance of differences</th></tr><tr><th>Jacc</th><th>P (≤ Jacc)</th><th>P (≥ Jacc)</th></tr></thead><tbody><tr><td>Degree centrality</td><td>0.12</td><td>0.001***</td><td>1.00</td></tr><tr><td>Betweenness centrality</td><td>0.14</td><td>0.001***</td><td>1.00</td></tr><tr><td>Closeness centrality</td><td>0.14</td><td>0.001***</td><td>1.00</td></tr><tr><td>Eigenvector centrality</td><td>0.14</td><td>0.001***</td><td>1.00</td></tr><tr><td>Hub taxa</td><td>0.32</td><td>0.46</td><td>0.58</td></tr></tbody></table>

Notes: *P* (≤ Jacc) and *P* (≥ Jacc) refer to the probabilities of observed Jaccard similarity coefficient (Jacc) that is either less than or equal to, or greater than or equal to, a Jacc calculated for randomly generated networks, respectively. Significance code: ***P = 0.001.

(Fig. 4A). Cascading removal had the most profound effect on network connectivity (Fig. 4B). In the Ap-infected network, 0.05 fraction of nodes removed was enough to achieve a connectivity loss of 0.80 (Fig. 4B). In

contrast, the Ap-uninfected network required the removal of a larger fraction, specifically 0.38 of nodes, to reach the same level of connectivity loss (Fig. 4B). The results also showed that the networks of uninfected ticks maintain a larger largest connected component (LCC) (Fig. 4C) and have shorter average path length (APL) on average (Fig. 4D) compared to the infected ones. This suggests that the infection reduces network robustness.

However, after adding nodes, the Ap-infected network showed a notable improvement in the LCC (Fig. 4C) and a reduction in the APL (Fig. 4D), indicating a recovery of robustness. In contrast, the Ap-uninfected network showed an increase in the average path length after adding nodes (Fig. 4D), which implies a slight decrease in robustness. Overall, these findings point out the vulnerabilities in microbial communities of ticks infected with *A. phagocytophilum* as well as their ability to regain robustness after eventual microbial invasions.

## 4. Discussion

The colonization resistance of tick microbiota may play a pivotal role in shaping the dynamics of tick-borne diseases. In this study, we propose that infection with *A. phagocytophilum* significantly restructures the tick microbial community, potentially leading to reduced colonization resistance. To explore this hypothesis, we employed a network-based approach, with a specific focus on network robustness, to assess the impact of *A. phagocytophilum* on colonization resistance.

We observed no significant differences in alpha diversity metrics, including richness, evenness, and phylogenetic diversity, or in beta diversity using the Bray-Curtis dissimilarity index between the Ap-uninfected and Ap-infected groups. However, in contrast to our findings, Abraham et al. (2017) observed beta diversity differences between the Ap-uninfected and Ap-infected groups using the weighted UniFrac metric. This difference underscores the impact of employing different analytical pipelines in the analysis of microbial community data. For instance, while we employed amplicon sequence variants (ASVs), Abraham et al. (2017) used operational taxonomic units (OTUs). ASVs

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Comparison of network robustness after node addition and removal. A Connectivity loss was measured in different attack scenarios for Ap-uninfected and Ap-infected networks: betweenness (red), cascading (green), degree (orange), and random (blue). B Cascading removal representation for Ap-uninfected vs Ap-infected. Additionally, the impact of node addition on network robustness was evaluated by employing two measures largest connected component (LCC) and average path length (APL) (C, D). For each network, a total of 1000 nodes were added. The value reached for LCC size (C) and the APL (D) for *A. phagocytophilum* after each node addition was plotted. Tick groups (uninfected and infected) are represented by blue and red, respectively.</figcaption></figure>
capture finer-scale diversity due to their resolution at the sequence level, whereas OTUs may lump together closely related sequences into the same unit, potentially obscuring underlying diversity patterns (Chiarello et al., 2022). This choice of pipeline can also influence alpha- and beta diversities, altering the ecological signal detected, particularly regarding presence/absence indices like richness index and unweighted UniFrac (Chiarello et al., 2022). This discrepancy in diversity estimates may be associated with the lower evenness of relative abundances and a higher number of rare units observed with OTUs compared to ASVs (Chiarello et al., 2022). Furthermore, the correlation between microbial communities varied across different beta diversity indices (Nearing et al., 2018; Chiarello et al., 2022). For instance, weighted UniFrac emphasizes phylogenetic relationships between microbial taxa, revealing patterns of community assembly driven by evolutionary processes. This metric is particularly useful for understanding how closely related taxa are distributed across different samples (Parks and Beiko, 2013). On the other hand, Bray-Curtis focuses on the abundance and presence of taxa, providing insights into compositional differences between samples (Ricotta and Podani, 2017). It is valuable for identifying which taxa contribute most to dissimilarities between communities and can highlight ecological factors shaping community composition. The contrasting results obtained using Bray-Curtis dissimilarity and weighted UniFrac metrics highlight the substantial impact of employing different analytical pipelines in microbial community analysis. The choice between ASVs and OTUs as well as the selection of diversity indices can significantly influence the ecological signals detected, particularly in relation to presence/absence metrics and phylogenetic relationships.

The lack of significant changes in alpha diversity suggests that overall microbial richness and evenness may remain relatively stable despite the presence of *A. phagocythophilum*. However, the observed change in composition (Fig. 2C) after the infection suggests that *A. phagocythophilum* could affect the occurrence of specific taxa in the tick microbiota, which may have important implications for colonization resistance. Colonization resistance is not dependent on a single species but rather on the associations of multiple bacteria cohesively living in a community (Spragge et al., 2023). These findings suggest that while the overall diversity of the microbiota may not be impacted, specific microbial communities crucial for colonization resistance could be compromised by infection with *A. phagocythophilum*, as evidenced by the present network analysis.

The observed composition changes could be associated with immunodulation mediated by *A. phagocytophilum* infection (Abraham et al., 2017). *Ixodes scapularis* antifreeze glycoprotein (IAFGP), a protein induced by *A. phagocytophilum* in ticks, selectively binds the terminal d-alanine residue of Gram-positive bacteria and inhibits biofilm formation among Gram-positive pathogens like *Staphylococcus aureus*, correlating with a reduction in Gram-positive biofilm-forming species, including enterococci, during *A. phagocytophilum* infection (Abraham et al., 2017). The differential role of IAFGP on Gram-positive versus Gram-negative bacteria provides additional insights into the observed effects of IAFGP on the tick microbiota and specific genera, particularly in the context of *A. phagocytophilum* infection in ticks. Notably, silencing *iafgp* impaired *A. phagocytophilum* colonization of the tick gut (Abraham et al., 2017), suggesting that reduced colonization resistance mediated by IAFGP is essential for pathogen infection. Despite the stability in the overall number of observed features, the altered microbial composition may compromise functional redundancy (Estrada-Peña et al., 2020a), potentially weakening the tick microbiota's ability to resist colonization by pathogens.

As previously reported for other intracellular pathogens (Maitre et al., 2022, 2023), infection reduces network complexity, which agrees with our findings of *A. phagocytophilum* infection being linked to a reduction in the number of nodes and edges in the microbial networks. The varying proportions of unique and shared nodes between the Ap-uninfected and Ap-infected networks suggest that infection changes

the composition of the bacterial community resulting in changes in network structure, and interactions among bacterial species within the community assembly. Nevertheless, the presence of shared nodes between the Ap-uninfected and Ap-infected networks, indicates that certain microbial taxa persist and maintain their interactions despite the presence of pathogens. These shared taxa may play important roles in maintaining the stability of the community (Seal et al., 2021; Paulino et al., 2023).

The observation that only positive interactions remained in the core network suggests a cooperative relationship among the microbial taxa (Chow et al., 2014; Fountain-Jones et al., 2023). Positive interactions can indicate functional associations where organisms perform similar or complementary functions or interactions shaped by interspecies cross-feeding (Fuhrman and Steele, 2008; Eiler et al., 2012; Chow et al., 2014; Lejal et al., 2021). Such cooperative interactions are thought to contribute to the overall stability and functionality of the microbial community, potentially enhancing colonization resistance (He et al., 2014; Aivelo et al., 2019; Oña and Kost, 2022). Negative interactions can reflect competition and niche partitioning among microorganisms (Engel and Moran, 2013). While negative interactions may disrupt the stability of the network to some extent, their presence suggests ongoing ecological dynamics within the community that could benefit the tick. For example, antagonistic interactions between gut microorganisms in insects can potentially have protective functions against pathogens (Engel and Moran, 2013).

Microbial communities, like those within tick populations, play a crucial role in the ecosystem balance, affecting both disease transmission and resistance (Wu-Chuang et al., 2021, 2023; Pavanelo et al., 2023). The stability of microbial networks in the face of disturbances - such as infections - can significantly influence their function and resilience (Shade et al., 2012; Estrada-Peña et al., 2020b). In our study, using network analysis and principles of percolation theory, we assess the microbial networks' resistance to various types of disruptions, including node removals (Cohen et al., 2000, 2001) and additions (Cohen and Havlin, 2009), and how these disruptions affect network connectivity and efficiency. Network robustness reflects a system's ability to maintain its connectivity and function despite disturbances (Cohen et al., 2000, 2001; Cohen and Havlin, 2009). Our findings indicate that infection with *A. phagocytophilum* compromises the network robustness, demonstrated by a more significant loss of connectivity following node removals. This suggests that the infection makes the microbial community more vulnerable to disturbances, with a potential decrease in colonization resistance.

On the other hand, after the introduction of new nodes (simulating the addition of new microbial species), we observed that infected tick networks exhibited lower LCC and higher APL values. The LCC is a measure of network cohesiveness, representing the largest group of interconnected nodes (Barabási and Pósfai, 2016; Kitsak et al., 2018), whereas APL indicates the network's efficiency, with shorter paths meaning quicker and more efficient communication between nodes (Barabási and Pósfai, 2016). The changes observed in these metrics suggest that, while *A. phagocytophilum* infection initially disrupts network structure, making it less cohesive and efficient, the addition of new nodes helps to somewhat mitigate these effects by enhancing connectivity and communication efficiency. This is in agreement with previous reports showing that less robust networks are more susceptible to new associations (Kwon and Cho, 2008; Scheffer et al., 2012; de Morais and Antunes, 2019), making it easier for potential pathogens or commensals to establish and persist within the microbiota. Nevertheless, it is important to acknowledge the limitations inherent to the use of simulated data, such as the inclusion of random nodes. While this approach enables us to explore potential shifts in network dynamics and theoretical resilience, it does not directly replicate biological realities. Consequently, the outcomes should be interpreted as indicative of possible structural responses to hypothetical changes rather than precise predictions of microbial interactions.
The compromised colonization resistance of tick microbiota caused by *A. phagocytophilum* rendering them more susceptible to pathogenic invasion may have implications in the natural environment where tick populations often harbor multiple tick-borne pathogens (TBPs). Co-infection with multiple pathogens is a common phenomenon in ticks (Nieto et al., 2018; Civitello et al., 2010; Moutailler et al., 2016; Hoffmann et al., 2023). It has been suggested that the presence of one pathogen may facilitate the establishment or proliferation of another within the tick microbiota (Gall et al., 2016; Sun et al., 2020). Our results suggest a novel mechanism by which one pathogen infection may decrease the colonization resistance of tick microbiota, favouring subsequent invasion events, that may include TBPs and/or commensals.

Using insights from human microbiota research (Spragge et al., 2023), we can understand potential mechanisms by which compromised colonization resistance in ticks contributes to co-infection dynamics. In humans, diverse microbial communities limit pathogen growth by consuming nutrients that pathogens need, a principle known as nutrient blocking (Spragge et al., 2023). Key species within these communities, especially those closely related to pathogens, are crucial because they increase the overlap in nutrient use, making it harder for pathogens to find the resources they need to grow. Applying this concept to ticks, a compromised colonization resistance - due to reduced robustness and weak community assembly - lessens nutrient blocking, making ticks more susceptible to co-infections. This is because a less diverse microbiome cannot effectively consume all available nutrients, leaving more resources for pathogens. Moreover, the lack of key species that directly compete with pathogens for nutrients may further ease the establishment and spread of these pathogens within the tick.

This may have an impact on disease ecology, as for example, *I. ricinus* and *I. scapularis* ticks are capable of harboring multiple pathogens, thereby increasing the likelihood of co-transmission to humans or animals (Lou et al., 2017). This phenomenon may apply to TBPs other than *A. phagocytophilum*. For example, infection with *B. burgdorferi* has been shown to enhance the transmission of *Babesia microti*, by increasing its basic reproduction number above the threshold for persistence (Diuk--Wasser et al., 2016). Furthermore, co-infection of *B. burgdorferi* and *Rickettsia* spp. in *Ixodes* nymphs results in bacterial replication rates higher than in single infections (Raulf et al., 2018; Sun et al., 2020). Further research could test empirically whether enhanced transmission of *B. microti* is due to a decreased colonization resistance induced by a primary infection with *B. burgdorferi*. While co-infection has been reported to enhance pathogen transmission, it is important to note that pathogens interacting within a population may exhibit complex dynamics, including mutual promotion, competition, or independence (Lou et al., 2017; Sun et al., 2020). Therefore, beyond pairwise interactions, facilitation, and competition summatory effects resulting in reduced colonization resistance may determine specific microbe invasion outcomes.

## 5. Conclusions

Our study revealed that infection with *A. phagocytophilum* alters the microbial community within ticks, potentially affecting their resistance to colonization by other pathogens. While overall microbial diversity remains stable, the composition and network structure of the microbiome showed significant changes post-infection. Our network analysis indicates that *A. phagocytophilum* infection reduces the complexity and robustness of the microbial network, possibly making ticks more susceptible to further pathogenic invasions. These changes suggest a decrease in the microbiota's ability to prevent additional infections, highlighting a potential mechanism for increased co-infection rates in tick populations. This study contributes to understanding how tick-borne diseases might spread more easily due to changes in tick microbiota caused by infection. It suggests that maintaining a healthy and diverse microbial community within ticks could be crucial for controlling the transmission of tick-borne diseases. Further research is

necessary to explore the broader implications of these findings across different tick species and pathogens.

## Funding

BIPAR was funded 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 was supported by the Programa Nacional de Becas de Postgrado en el Exterior "Don Carlos Antonio López" (Grant No. 205/2018). Apolline Maître is supported by the "Collectivité de Corse", grant: "Formations superieures" (SGCE-RAPPORT No. 0300).

## Ethical approval

Not applicable.

## CRediT authorship contribution statement

**Lianet Abuin-Denis:** Investigation, Formal analysis, Visualization, Writing – original draft, Writing – review & editing. **Elianne Piloto-Sardinas:** Formal analysis, Visualization, Writing – original draft, Writing – review & editing. **Apolline Maître:** Visualization, Supervision, Writing – review & editing. **Alejandra Wu-Chuang:** Investigation, Supervision, Formal analysis, Writing – review & editing. **Lourdes Mateos-Hernández:** Investigation, Data curation, Formal analysis, Writing – review & editing. **Dasiel Obregón:** Investigation, Methodology, Writing – review & editing. **Belkis Corona-González:** Supervision, Writing – review & editing. **Andréa Cristina Fogaça:** Supervision, Writing – review & editing. **Vaidas Palinauskas:** Supervision, Methodology, Writing – review & editing. **Justè Azelytè:** Formal analysis, Visualization, Writing – review & editing. **Alina Rodríguez-Mallon:** Supervision, Writing – review & editing. **Alejandro Cabezas-Cruz:** Supervision, Visualization, Resources, Writing – original draft, Writing – review & editing. All authors read and approved the final manuscript.

## Declaration of competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Given their role as Guest Editor, Alejandro Cabezas-Cruz had no involvement in the peer review of this article and has no access to information regarding its peer review. Full responsibility for the editorial process for this was delegated to Professor Aneta Kostadinova (Editor-in-Chief).

## Data availability

The data supporting the conclusions of this article are included within the article and its supplementary files. All the datasets shown in the present study can be found at the SRA repository https://www.ncbi.nlm.nih.gov/sra (accession numbers: PRJNA35373 and PRJNA803003).

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.crpvbd.2024.100177.

## References

* Abraham, N.M., Liu, L., Jutras, B.L., Yadav, A.K., Narasimhan, S., Gopalakrishnan, V., et al., 2017. Pathogen-mediated manipulation of arthropod microbiota to promote infection. Proc. Natl. Acad. Sci. USA 114, E781-E790. https://doi.org/10.1073/pnas.1613422114.

* Aivelo, T., Norberg, A., Tschirren, B., 2019. Bacterial microbiota composition of *Ixodes ricinus* ticks: The role of environmental variation, tick characteristics and microbial interactions. PeerJ 7, e8217. https://doi.org/10.7717/peerj.8217.
Barabási, A.-L., Pósfai, M., 2016. Network Science. Cambridge University Press, Cambridge.

Bastian, M., Heymann, S., Jacomy, M., 2009. Gephi: an open source software for exploring and manipulating networks. Third International AAAI Conference on Weblogs and Social Media, Vol. 3., https://doi.org/10.1609/icwsm.v3i1.13937. No. 1.

Binetruy, F., Buysse, M., Lejarre, Q., Barosi, R., Villa, M., Rahola, N., et al., 2020. Microbial community structure reveals instability of nutritional symbiosis during the evolutionary radiation of *Amblyomma* ticks. Mol. Ecol. 29, 1016–1029. https://doi.org/10.1111/med.15373.

Bolyen, E., Rideout, J.R., Dillon, M.R., Bokulich, N.A., Abnet, C.C., Al-Ghalith, G.A., et al., 2019. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857. https://doi.org/10.1038/s41587-019-0209-9.

Bonnet, S.I., Binetruy, F., Hernández-Jargún, A.M., Duron, O., 2017. The tick microbiome: Why non-pathogenic microorganisms matter in tick biology and pathogen transmission. Front. Cell. Infect. Microbiol. 7, 236. https://doi.org/10.3389/fcimb.2017.00236.

Bray, J.R., Curtis, J.T., 1957. An ordination of the upland forest communities of southern Wisconsin. Ecol. Monogr. 27, 325–349. https://doi.org/10.2307/1942268.

Cabezas-Cruz, A., 2023. Grand challenges in arachnid microbiome and diseases. Front. Arachnid Sci. 2, 1215831. https://doi.org/10.3389/frchs.2023.1215831.

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. Nat. Methods 13, 581–583. https://doi.org/10.1038/nmeth.3869.

Chiarello, M., McCauley, M., Villéger, S., Jackson, C.R., 2022. Ranking the biases: The choice of OTUs vs. ASVs in 16S rRNA amplicon data analysis has stronger effects on diversity measures than rarefaction and OTU identity threshold. PLoS One 17, e0264443. https://doi.org/10.1371/journal.pone.0264443.

Chow, C.-E.T., Kim, D.Y., Sachdeva, R., Caron, D.A., Fuhrman, J.A., 2014. Top-down controls on bacterial community structure: Microbial network analysis of bacteria, T4-like viruses and protists. ISME J. 8, 816–829. https://doi.org/10.1038/ismej.2013.199.

Civitello, D.J., Rynkiewicz, E., Clay, K., 2010. Meta-analysis of co-infections in ticks. Israel J. Ecol. Evolut. 56, 417–431. https://doi.org/10.1560/IJEE.56.3.4.417.

Cohen, R., Erezn, K., Ben-Avraham, D., Havlin, S., 2000. Resilience of the Internet to random breakdowns. Phys. Rev. Lett. 85, 4626–4628. https://doi.org/10.1103/PhysRevLett.85.4626.

Cohen, R., Erezn, K., Ben-Avraham, D., Havlin, S., 2001. Breakdown of the Internet under intentional attack. Phys. Rev. Lett. 86, 3682–3685. https://doi.org/10.1103/PhysRevLett.86.3682.

Cohen, R., Havlin, S., 2009. Percolation in complex networks. In: Meyers, R. (Ed.), Encyclopedia of Complexity and Systems Science. Springer, New York, NY. https://doi.org/10.1007/978-0-387-30440-3_383.

de la Fuente, J., Antunes, S., Bonnet, S., Cabezas-Cruz, A., Domingos, A.G., Estrada-Peña, A., et al., 2017. Tick-pathogen interactions and vector competence: Identification of molecular drivers for tick-borne diseases. Front. Cell. Infect. Microbiol. 7, 114. https://doi.org/10.3389/fcimb.2017.00114.

De Morais, U.L., Antunes, L., 2019. A look at the way we look at complex networks robustness and resilience. arXiv. https://doi.org/10.48550/ArXiv.1909.06448.

DeSantis, T.Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E.L., Keller, K., et al., 2006. Greeninges, a Chimeria-checked 16S rRNA gene database and workbench compatible with ARB. Appl. Environ. Microbiol. 72, 5069–5072. https://doi.org/10.1128/AEM.03006-05.

Díaz-Sánchez, S., Estrada-Peña, A., Cabezas-Cruz, A., de la Fuente, J., 2019. Evolutionary insights into the tick hologenome. Trends Parasitol. 35, 725–737. https://doi.org/10.1016/j.pt.2019.06.014.

Diuk-Wasser, M.A., Vannier, E., Krause, P.J., 2016. Coinfection by *Ixodes* tick-borne pathogens: Ecological, epidemiological, and clinical consequences. Trends Parasitol. 32, 30–42. https://doi.org/10.1016/j.pt.2015.09.008.

Ducarmon, Q.R., Zwickitt, R.D., Hornung, B.V.H., van Schaik, W., Young, V.B., Kuijper, E.J., 2019. Gut microbiota and colonization resistance against bacterial enteric infection. Microbiol. Mol. Biol. Rev. 83, e00007 https://doi.org/10.1128/MMBR.00007-19, 19.

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

Engel, P., Moran, N.A., 2013. The gut microbiota of insects - diversity in structure and function. FEMS Microbiol. Rev. 37, 699–735. https://doi.org/10.1111/1574-6976.12025.

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

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

Faith, D.P., 1992. Conservation evaluation and phylogenetic diversity. Biol. Conserv. 61, 1–10. https://doi.org/10.1016/0006-3207(92)91201-3.

Faust, K., Raes, J., 2012. Microbial interactions: From networks to models. Nat. Rev. Microbiol. 10, 538–550. https://doi.org/10.1038/nrmicro2832.

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, e67019. https://doi.org/10.1371/journal.pone.0067019.

Fogaça, A.C., Sousa, G., Pavanelo, D.B., Esteves, E., Martins, L.A., Urbanová, V., et al., 2021. Tick immune system: What is known, the interconnections, the gaps, and the

challenges. Front. Immunol. 12, 628054 https://doi.org/10.3389/fimmu.2021.628054.

Fountain-Jones, N.M., Khoo, B.S., Rau, A., Berman, J.D., Burton, E.N., Oliver, J.D., 2023. Positive associations matter: Microbial relationships drive tick microbiome composition. Mol. Ecol. 32, 4078–4092. https://doi.org/10.1111/mec.16985.

Freitas, S., Yang, D., Kumar, S., Tong, H., Chau, D.H., 2020. Evaluating graph vulnerability and robustness using TIGER. https://doi.org/10.48550/ARXIV.2006.05648.

Friedman, J., Alm, E.J., 2012. Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8, e1002687 https://doi.org/10.1371/journal.pcbi.1002687.

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

Gall, C.A., Reif, K.E., Scoles, G.A., Mason, K.L., Mousel, M., Noh, S.M., Brayton, K.A., 2016. The bacterial microbiome of *Dermacentor andersoni* ticks influences pathogen susceptibility. ISME J. 10, 1846–1855. https://doi.org/10.1038/ismej.2015.266.

He, X., McLean, J.S., Guo, L., Lux, R., Shi, W., 2014. The social structure of microbial community involved in colonization resistance. ISME J. 8, 564–574. https://doi.org/10.1038/ismej.2013.172.

Hoffmann, A., Müller, T., Fingerle, V., Silaghi, C., Noll, M., 2023. Co-infection of potential tick-borne pathogens of the order Rickettsiales and *Borrelia burgdorferi* s.l. and their link to season and area in Germany. Microorganisms 11, 157. https://doi.org/10.3390/microorganisms11010157.

Hussain, S., Perveen, N., Hussain, A., Song, B., Aziz, M.U., Zeb, J., Li, J., et al., 2022. The symbiotic continuum within ticks: Opportunities for disease control. Front. Microbiol. 13, 854803 https://doi.org/10.3389/fmicb.2022.854803.

Karita, Y., Limmer, D.T., Hallatschek, O., 2022. Scale-dependent tipping points of bacterial colonization resistance. Proc. Natl. Acad. Sci. USA 119, e2115496119. https://doi.org/10.1073/pnas.2115496119.

Katoh, K., 2002. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30, 3059–3066. https://doi.org/10.1093/nar/gk436.

Kitsak, M., Ganin, A.A., Eisenberg, D.A., Krapivsky, P.L., Krioukov, D., Alderson, D.L., Linkov, I., 2018. Stability of a giant connected component in a complex network. Phys. Rev. Lett. 97, 012309 https://doi.org/10.1103/PhysRevE.97.012309.

Kitsou, C., Foor, S.D., Dutta, S., Bista, S., Pal, U., 2021. Tick gut barriers impacting tick-microbe interactions and pathogen persistence. Mol. Microbiol. 116, 1241–1248. https://doi.org/10.1111/mmi.14822.

Kocan, M., De La Fuente, F.J., Cabezas-Cruz, A., 2015. The genus *Anaplasma*: New challenges after reclassification. Rev. Sci. Tech. 34, 577–586. https://doi.org/10.20506/rtst.34.2.2381.

Kurokawa, C., Lynn, G.E., Pedra, J.H.F., Pal, U., Narasimhan, S., Fikrig, E., 2020. Interactions between *Borrelia burgdorferi* and ticks. Nat. Rev. Microbiol. 18, 587–600. https://doi.org/10.1038/s41579-020-04000-5.

Kurtz, Z.D., Müller, C.L., Miraldi, E.R., Littman, D.R., Blaser, M.J., Bonneau, R.A., 2015. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput. Biol. 11, e1004226 https://doi.org/10.1371/journal.pcbi.1004226.

Kwon, Y.-K., Cho, K.-H., 2008. Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics. Bioinformatics 24, 987–994. https://doi.org/10.1093/bioinformatics/btn060.

Lejal, E., Chiquet, J., Aubert, J., Robin, S., Estrada-Peña, A., Rue, O., et al., 2021. Temporal patterns in *Ixodes ricinus* microbial communities: An insight into tick-borne microbiome interactions. Microbiome 9, 153. https://doi.org/10.1186/s40168-021-01051-8.

Lhomme, S., 2015. Analyse spatiale de la structure des réseaux techniques dans un contexte de risques. Cybergeo 2015, 711. https://doi.org/10.4000/cybergeo.26763.

Lou, Y., Liu, L., Gao, D., 2017. Modelling co-infection of *Ixodes* tick-borne pathogens. Math. Biosci. Eng. 14, 1301–1316. https://doi.org/10.3934/mbe.2017067.

Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Foucault-Simonin, A., Moutailler, S., Paoli, J.-C., et al., 2022. *Rickettsia helvetica* infection is associated with microbiome modulation in *Ixodes ricinus* collected from humans in Serbia. Sci. Rep. 12, 11464 https://doi.org/10.1038/s41598-022-15681-x.

Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Piloto-Sardiñas, E., Foucault-Simonin, A., Cicculi, V., et al., 2023. Rickettsial pathogens drive microbiota assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks. Mol. Ecol. 32, 4660–4676. https://doi.org/10.1111/mec.17058.

Mateos-Hernández, L., Maitre, A., Abuin-Denis, L., Obregon, D., Martin, E., Luis, P., et al., 2023. Hierarchical shift of the *Aedes albopictus* microbiota caused by antimicrobiota vaccine increases fecundity and egg-hatching rate in female mosquitoes. FEMS Microbiol. Ecol. 99, fiad140 https://doi.org/10.1093/femsec/fiad140.

Mateos-Hernández, L., Obregón, D., Maye, J., Borneres, J., Versille, N., de la Fuente, J., et al., 2020. Anti-tick microbiota vaccine impacts *Ixodes ricinus* performance during feeding. Vaccines 8, 702. https://doi.org/10.3390/vaccines8040702.

Mateos-Hernández, L., Obregón, D., Wu-Chuang, A., Maye, J., Bornéres, J., Versillé, N., et al., 2021. Anti-microbiota vaccines modulate the tick microbiome in a taxon-specific manner. Front. Immunol. 12, 704621 https://doi.org/10.3389/fimmu.2021.704621.

Moutailler, S., Valiente Moro, C., Vaumourin, E., Michelet, L., Tran, F.H., Devillers, E., et al., 2016. Co-infection of ticks: The rule rather than the exception. PLoS Neglected Trop. Dis. 10, e0004539 https://doi.org/10.1371/journal.pntd.0004539.

Mullineaux-Sanders, C., Suez, J., Elinav, E., Frankel, G., 2018. Sieving through gut models of colonization resistance. Nat. Microbiol. 3, 132–140. https://doi.org/10.1038/s41564-017-0095-1.

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

Nearing, J.T., Douglas, G.M., Comeau, A.M., Langille, M.G.I., 2018. Denoising the
denoisers: An independent evaluation of microbiome sequence error-correction
approaches. PeerJ 6, e5364. https://doi.org/10.7717/peerj.5364.

Nieto, N.C., Tanner Porter, W., Wachara, J.C., Lowrey, T.J., Martin, L., Motyka, P.J., Salkeld, D.J., 2018. Using citizen science to describe the prevalence and distribution of tick bite and exposure to tick-borne diseases in the United States. PLoS One 13, e0199644. https://doi.org/10.1371/journal.pone.0199644.

Obregón, D., Bard, E., Abrial, D., Estrada-Peña, A., Cabezas-Cruz, A., 2019. Sex-specific linkages between taxonomic and functional profiles of tick gut microbiomes. Front. Cell. Infect. Microbiol. 9, 298. https://doi.org/10.3389/fcimb.2019.00298.

Oña, L., Kost, C., 2022. Cooperation increases robustness to ecological disturbance in microbial cross-feeding networks. Ecol. Lett. 25, 1410-1420. https://doi.org/10.1111/ele.14006.

Parks, D.H., Beiko, R.G., 2013. Measures of phylogenetic differentiation provide robust and complementary insights into microbial communities. ISME J. 7, 173-183. https://doi.org/10.1038/ismej.2012.88.

Paulino, P.G., Abuin-Denis, L., Maitre, A., Piloto-Sardiñas, E., Obregon, D., Santos, H.A., Cabezas-Cruz, A., 2023. Dissecting the impact of *Anaplasma phagocytophilum* infection on functional networks and community stability of the tick microbiome. Int. Microbiol. https://doi.org/10.1007/s10123-023-00473-8 (Online ahead of print).

Pavanelo, D.B., Piloto-Sardiñas, E., Maitre, A., Abuin-Denis, L., Kopáček, P., Cabezas-Cruz, A., Fogaça, A.C., 2023. Arthropod microbiota: Shaping pathogen establishment and enabling control. Front. Arachnid Sci. 2, 1297733 https://doi.org/10.3389/frchs.2023.1297733.

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 Bioinf. 22, bbaa290. https://doi.org/10.1093/bib/bbaa290.

Pielou, E.C., 1966. The measurement of diversity in different types of biological collections. J. Theor. Biol. 13, 131-144. https://doi.org/10.1016/0022-5193(66)90013-0.

Price, M.N., Dehal, P.S., Arkin, A.P., 2010. FastTree 2 – approximately maximum-likelihood trees for large alignments. PLoS One 5, e9490. https://doi.org/10.1371/journal.pone.0009490.

Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., et al., 2013. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41 (Database issue), D590-D596. https://doi.org/10.1093/nar/gks1219.

R Core Team, 2023. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

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

Real, R., Vargas, J.M., 1996. The probabilistic basis of Jaccard's index of similarity. Syst. Biol. 45, 380-385. https://doi.org/10.1093/sysbio/45.3.380.

Ricotta, C., Podani, J., 2017. On some properties of the Bray-Curtis dissimilarity and their ecological meaning. Ecol. Complex. 31, 201-205. https://doi.org/10.1016/j.ecocom.2017.07.003.

Röttjers, L., Faust, K., 2018. From hairballs to hypotheses - biological insights from microbial networks. FEMS Microbiol. Rev. 42, 761-780. https://doi.org/10.1093/femsre/fuy030.

Röttjers, L., Vandeputte, D., Raes, J., Faust, K., 2021. Null-model-based network comparison reveals core associations. ISME Communications 1, 36. https://doi.org/10.1038/s43705-021-00036-w.

RStudio Team, 2020. RStudio: Integrated development for R. RStudio, PBC, Boston, MA, USA. http://www.rstudio.org/.

Scheffer, M., Carpenter, S.R., Lenton, T.M., Bascompte, J., Brock, W., Dakos, V., et al., 2012. Anticipating critical transitions. Science 338, 344-348. https://doi.org/10.1126/science.1225244.

Seal, S., Dharmarajan, G., Khan, I., 2021. Evolution of pathogen tolerance and emerging infections: A missing experimental paradigm. Elife 10, e68874. https://doi.org/10.7554/eLife.68874.

Shade, A., Peter, H., Allison, S.D., Baho, D.L., Berga, M., Bürgmann, H., et al., 2012. Fundamentals of microbial community resistance and resilience. Front. Microbiol. 3, 417. https://doi.org/10.3389/fmicb.2012.00417.

Spragge, F., Bakkeren, E., Jahn, M.T., Araujo, E., Pearson, C.F., Wang, X., et al., 2023. Microbiome diversity protects against pathogens by nutrient blocking. Science 382, eadj3502. https://doi.org/10.1126/science.adj3502.

Stacy, A., Andrade-Oliveira, V., McCulloch, J.A., Hild, B., Oh, J.H., Perez-Chaparro, P.J., et al., 2021. Infection trains the host for microbiota-enhanced resistance to pathogens. Cell 184, 615-627.e17. https://doi.org/10.1016/j.cell.2020.12.011.

Sun, B., Zhang, X., Tosato, M., 2020. Effects of coinfection on the dynamics of two pathogens in a tick-host infection model. Complexity 2020, 5615173. https://doi.org/10.1155/2020/5615173.

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. Front. Cell. Infect. Microbiol. 11, 713466 https://doi.org/10.3389/fcimb.2021.713466.

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. Curr. Res. Parasitol. Vector Borne Dis. 1, 100036 https://doi.org/10.1016/j.crpvbd.2021.100036.

Wu-Chuang, A., Mateos-Hernandez, L., Maitre, A., Rego, R.O.M., Síma, R., Porcelli, S., et al., 2023. Microbiota perturbation by anti-microbiota vaccine reduces the colonization of *Borrelia afzelii* in *Ixodes ricinus*. Microbiome 11, 151. https://doi.org/10.1186/s40168-023-01599-7.