# Differential nested patterns of Anaplasma marginale and Coxiella-like endosymbiont across Rhipicephalus microplus ontogeny

Lianet Abuin-Denis, Elianne Piloto-Sardiñas, Apolline Maitre, Alejandra Wu-Chuang, Lourdes Mateos-Hernández, Patrícia Gonzaga Paulino, Yamil Bello, Frank Ledesma Bravo, Anays Alvarez Gutierrez, Rafmary Rodríguez Fernández, et al.

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

► To cite this version:

Lianet Abuin-Denis, Elianne Piloto-Sardiñas, Apolline Maitre, Alejandra Wu-Chuang, Lourdes Mateos-Hernández, et al.. Differential nested patterns of Anaplasma marginale and Coxiella-like endosymbiont across Rhipicephalus microplus ontogeny. Microbiological Research, 2024, 286, pp.127790. 10.1016/j.micres.2024.127790. hal-04633489

HAL Id: hal-04633489

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

Submitted on 3 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Differential nested patterns of *Anaplasma marginale* and *Coxiella*-like endosymbiont across *Rhipicephalus microplus* ontogeny

Lianet Abuin-Denis<sup>a, b</sup>, Elianne Piloto-Sardiñas<sup>b, c</sup>, Apolline Maitre<sup>b, d, e</sup>, Alejandra Wu-Chuang<sup>b</sup>, Lourdes Mateos-Hernández<sup>b</sup>, Patrícia Gonzaga Paulino<sup>f</sup>, Yamil Bello<sup>a</sup>, Frank Ledesma Bravo<sup>a</sup>, Anays Alvarez Gutierrez<sup>a</sup>, Rafmary Rodríguez Fernández<sup>g</sup>, Alier Fuentes Castillo<sup>g</sup>, Luis Méndez Mellor<sup>g</sup>, Angélique Foucault-Simonin<sup>b</sup>, Dasiel Obregon<sup>h</sup>, Mario Pablo Estrada-García<sup>a</sup>, Alina Rodríguez-Mallon<sup>a, *</sup>, Alejandro Cabezas-Cruz<sup>b, *</sup>

## ARTICLE INFO

Keywords:

Ticks

Tick-borne pathogens

Community assembly

## ABSTRACT

Understanding the intricate ecological interactions within the microbiome of arthropod vectors is crucial for elucidating disease transmission dynamics and developing effective control strategies. In this study, we investigated the ecological roles of *Coxiella*-like endosymbiont (CLE) and *Anaplasma marginale* across larval, nymphal, and adult stages of *Rhipicephalus microplus*. We hypothesized that CLE would show a stable, nested pattern reflecting co-evolution with the tick host, while *A. marginale* would exhibit a more dynamic, non-nested pattern influenced by environmental factors and host immune responses. Our findings revealed a stable, nested pattern characteristic of co-evolutionary mutualism for CLE, occurring in all developmental stages of the tick. Conversely, *A. marginale* exhibited variable occurrence but exerted significant influence on microbial community structure, challenging our initial hypotheses of its non-nested dynamics. Furthermore, *in silico* removal of both microbes from the co-occurrence networks altered network topology, underscoring their central roles in the *R. microplus* microbiome. Notably, competitive interactions between CLE and *A. marginale* were observed in nymphal network, potentially reflecting the impact of CLE on the pathogen transstadial-transmission. These findings shed light on the complex ecological dynamics within tick microbiomes and have implications for disease management strategies.

## 1. Introduction

Ticks are obligate hematophagous ectoparasites that serve as vectors for various pathogens, including bacteria, viruses, and protozoa, causing significant health threats to both humans and animals (de la Fuente et al., 2017). *Rhipicephalus microplus*, commonly known as the cattle tick, is one of the most economically important tick species globally due to its role in transmitting various pathogens including apicomplexan parasites

such as *Babesia bovis*, *Babesia bigemina*, *Theileria equi* (Guimarães et al., 1998), and *Theileria orientalis* (Almazán et al., 2022), and bacterial pathogens including *Anaplasma marginale* (De La Fournière et al., 2023; Piloto-Sardinas et al., 2023) and *Ehrlichia minasensis* (Cabezas-Cruz et al., 2016; Carvalho et al., 2016), an emerging pathogen (Cabezas-Cruz et al., 2019; Moura De Aguiar et al., 2019; Ramírez-Hernández et al., 2022). In addition to pathogens, ticks harbor complex microbial communities (Rojas-Jaimes et al., 2021; Segura et al., 2020), consisting of

<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), Corte 20250, France</footnote>

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

<footnote><sup>f</sup> Department of Epidemiology and Public Health, Federal Rural University of Rio de Janeiro (UFRRJ), Seropedica 23890-000, Brazil</footnote>

<footnote><sup>g</sup> National Laboratory of Parasitology, Ministry of Agriculture, Autopista San Antonio de los Baños, Km 112, San Antonio de los Baños, Artemisa 38100, Cuba</footnote>

<footnote><sup>h</sup> School of Environmental Sciences University of Guelph, Guelph, Ontario N1G 2W1, Canada</footnote>

<footnote>* Corresponding authors.</footnote>

<footnote>E-mail addresses: alina.rodriguez@cigb.edu.cu (A. Rodríguez-Mallon), alejandro.cabezas@vet-alfort.fr (A. Cabezas-Cruz).</footnote>
environmental-acquired mutualistic microbes (Hussain et al., 2022), and transovarially-transmitted endosymbionts (Binetruy et al., 2020), all referred to as **microbiome** (sous rature) (Cabezas-Cruz, 2023). The microbiome plays a critical role in vector survival, fitness (Mesquita et al., 2023), physiological processes and vector competence (Guizzo et al., 2022; Pavanelo et al., 2023). Understanding the assembly and dynamics of these microbial communities is crucial for elucidating vector-pathogen interactions and developing novel strategies for vector-borne disease management, such as anti-microbiota vaccines (Ażelytė et al., 2022; Maitre, Wu-Chuang, Ażelytė, et al., 2022; Mateos-Hernández et al., 2021; Pavanelo et al., 2023; Wu-Chuang et al., 2023).

Nestedness theory (Mariani et al., 2019) provides a valuable framework for understanding the diverse patterns of interactions between symbionts, commensals, and tick-borne pathogens (TBPs) in the tick microbiome. A perfect nesting is defined by a specific structural property: for any pair of nodes, if one node has a smaller degree than the other, then the neighborhood (i.e., set of adjacent nodes) of the less connected node is entirely encompassed within the neighborhood of the more connected node (Cobo-López et al., 2022; Song et al., 2017). In ecological networks, nestedness can also denote a structured pattern where the interactions of one species or taxon constitute a subset of the interactions of another species or taxon (Mariani et al., 2019). Nested patterns have been extensively studied in various ecological systems, ranging from plant-pollinator networks (Song et al., 2017) to host microbiome (Cobo-López et al., 2022), providing valuable insights into the organization and stability of communities. Applying the concept of nestedness to the tick microbiome involves investigating whether specific taxa exhibit a tendency to co-occur in a nested way within the network. This analysis can expose whether certain symbionts or pathogens consistently co-occur with a broader range of taxa, indicating potential ecological relationships. Additionally, it can elucidate whether there are specific subsets of taxa associated with particular environments or conditions, providing insights into the ecological dynamics of the tick microbiome in diverse settings. If a nestedness pattern is present, it suggests a structured organization in the interactions, which has been previously reported in tick microbiome networks (Maitre et al., 2023; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022). This could have implications for understanding the dynamics of symbionts, commensals, and TBPs within the tick microbiome.

In this study, we tracked two key members of the *R. microplus* microbial community, *Coxiella*-like endosymbionts (CLE) (Guizzo et al., 2022; Rialch et al., 2022) and *A. marginale* (De La Fournière et al., 2023; Piloto-Sardiñas et al., 2023), across the tick ontogeny, specifically the succession of life stages from larva to nymph to adult. CLE are intracellular bacteria that have been found in a wide range of arthropod vectors, including ticks (Zhong, 2012). These endosymbionts are believed to confer various fitness benefits to their arthropod hosts, such as enhanced reproduction and protection against other pathogens (Zhong, 2012). On the other hand, *A. marginale* is a Gram-negative intracellular bacterium and the causative agent of bovine anaplasmosis (Aubry and Geale, 2011; Kocan et al., 2015), a severe disease affecting cattle worldwide.

The divergent ecological roles and interactions exhibited by CLE and *A. marginale* within the tick host render the examination of their distinct nested patterns particularly interesting. We expect that the CLE, being mutualistic in nature (Zhong, 2012), will exhibit a nested pattern within the microbial community. These endosymbionts may form a core group of microbial taxa consistently found in ticks across their life stages, reflecting a stable and co-evolved relationship with the host. Conversely, *A. marginale*, being a pathogenic bacterium, may exhibit a less consistent or even non-nested pattern within the tick microbiome. The abundance and occurrence of *A. marginale* could be influenced by external factors such as the presence of infected hosts in the environment, variation in tick exposure to the pathogen, and the vector's immune response (Piloto-Sardiñas et al., 2023). Therefore, we hypothesize that

*A. marginale* will display a more dynamic and transient pattern within the microbial community, reflecting the pathogen's epidemiology and ecological dynamics.

By investigating the differential nested patterns of CLE and *A. marginale* in the *R. microplus* microbial community across the ontogeny of the tick, we aim to gain valuable insights into the complex interactions that shape vector microbiome assembly. This research has the potential to enhance our understanding of vector-pathogen dynamics and may pave the way for novel strategies in managing vector-borne diseases in livestock and public health settings.

## 2. Materials and Methods

### 2.1. Study design

For this study, *Rhipicephalus microplus* ticks from a colony established at the Cuban National Laboratory of Parasitology (Encinosa Guzmán et al., 2021) were fed on a single cattle infected with *Anaplasma marginale*. The ticks, while remaining on the animal, were allowed to progress through their entire life cycle. Ticks from each developmental stage were collected: larvae (n = 13), nymphs (n = 15), and adult males (n = 15) for subsequent analyses. Before DNA extraction, the collected ticks underwent a thorough washing process, which involved two rounds of washing in milliQ sterile water and one round in 70% ethanol. DNA from individual ticks was extracted and purified using the DNA Tissue and Blood Kit (QIAGEN, Germany) following the manufacturer's instructions. *Anaplasma marginale* and *Coxiella*-like endosymbionts (CLE) were detected on individual tick samples through qualitative PCR. High-throughput amplicon sequencing of 16 S rRNA gene was used to compare the microbial community composition and networks were used to study the microbial community assembly.

### 2.2. Detection of *A. marginale* and CLE in tick samples

The presence of *A. marginale* in both the cattle and *R. microplus* ticks was tested using specific primers (Table 1) in PCR reactions conducted in a thermal cycler (MinicyclerTM, MJ Research, Inc., USA). For the design of the primers to amplify the *A. marginale* genes *msp1α* and *msp5*, we made alignments based on sequences reported for these genes in the Genbank (https://www.ncbi.nlm.nih.gov/genbank/) and designed primers using Primer-BLAST (Ye et al., 2012). The PCR reactions were prepared in a final volume of 50 µL, consisting of 5 µL of complementary DNAs, 25 µL of Master Mix 2X (Promega, USA), 50 pmol of each primer and water to complete the volume. The reaction mixtures underwent an initial denaturation of 5 minutes at 94°C, followed by 40 cycles of denaturation-annealing-extension (30 seconds at 94°C, 30 seconds at 56°C, and 2 minutes at 74°C), with a final extension of 5 minutes at 72°C. The presence of CLE in *R. microplus* ticks was assessed by nested PCR using specific primers described by Duron et al. (2014) (Table 1).

### 2.3. 16 S rRNA sequencing for bacterial microbiome analysis

More than 900 ng of DNA at ≥ 20 ng/µL concentration was amplicon sequenced for the bacterial 16 S rRNA gene by Novogene Bioinformatics Technology Co. (London, UK). The sequencing process employed a single lane of the Illumina MiSeq system, generating 251-base paired-end reads from the V4 variable region of the 16 S rRNA gene. Bar-coded universal primers (515 F/806 R) were used for sequencing. The resulting raw 16 S rRNA sequences from tick samples were deposited in the SRA repository under Bioproject PRJNA1034019.

### 2.4. Processing of raw sequences

The raw sequences, represented as demultiplexed fastq files, underwent preprocessing and analysis through the Quantitative Insights Into Microbial Ecology 2 (QIIME2) pipeline (version 2019.1). Briefly, the
<table><caption>**Table 1**

Primer used to identify *A. marginale* and CLE in tick sample.</caption><thead><tr><th>Targeted genes</th><th>Primer sequences (5' - 3')</th><th>Expected amplicon size</th><th>References</th></tr></thead><tbody><tr><td><i>A. marginale</i> msp5</td><td>*F: CCATGGGCAAAGGCATTTTCAGCAAAATCGG<br>**R: CTCGAGAGAATTAAGCATGTGACCGCTG</td><td>555 pb</td><td>This study</td></tr><tr><td><i>A. marginale</i> msp1a</td><td>*F: GCCGAACCTGAGCATGAG<br>**R: CCAAACACGCACAGATAGCA</td><td>101 pb</td><td>This study</td></tr><tr><td><i>Coxiella</i>-like endosymbionts<br>16 S rRNA</td><td>*Cox16SF1: CGTAGGAATCTACCTTRTAGWGG<br>**Cox16SR1: ACTYYCCAACAGCTAGTTCTCA<br>*Cox16SF2: TGAGAACTAGCTGTTGGRRAGT<br>**Cox16SR2: GCCTACCCGCTTCTGGTACAATT</td><td>Cox16SF1/Cox16SR1: 719–813 bp<br>Cox16SF2/Cox16SR2: 624–625 bp</td><td>(Duron et al., 2014)</td></tr></tbody></table>

DADA2 software package (Callahan et al., 2016) was used via q2-dada2 to denoise the fastq files and merge mate reads. Subsequently, all amplicon sequence variants (ASVs) underwent alignment using Multiple Alignment using Fast Fourier Transform (MAFFT) (Katoh, 2002) via q2-alignment, and a phylogeny was constructed using FastTree 2 (Price et al., 2010) via q2-phylogeny. Taxonomy assignments of the ASVs were performed using the q2-feature-classifier plugin with the classify-sklearn naive Bayes taxonomy classifier based on the 16 S rRNA SILVA database (release 138), with a training accuracy of 99%; consequently, only the target sequences fragment was utilized. For visualizing taxonomic profiling by samples, QIIME2 taxa barplot functions were employed.

## 2.5. Diversity indexes and differential taxa abundance

To assess microbial diversity within-sample, we employed alpha diversity metrics, specifically Observed features, Faith's phylogenetic diversity (Faith's PD), and the Pielou's evenness index (DeSantis et al., 2006; Faith, 1992; Pielou, 1966). These metrics were computed using the q2-diversity plugin within the QIIME2 environment. The Observed features metrics counts the numbers of unique features present in a sample, providing insight into the richness of microbial taxa (DeSantis et al., 2006). Faith's phylogenetic diversity index evaluates the evolutionary distinctiveness of species within the community, offering a measure of phylogenetic diversity (Faith, 1992). Pielou's evenness index assesses the even distribution of individuals among species within the community, indicating the balance of species abundances relative to species richness (Pielou, 1966).

High Observed features values suggest different species in a given condition, suggesting a diverse community with a wide range of microbes. Higher evenness values indicate a more equitable distribution of individuals among taxa, while lower values suggest a more skewed or uneven distribution. High phylogenetic diversity values indicate that taxa are more distantly related to each other in evolutionary terms.

In evaluating between-sample diversity, beta diversity serves as a measure of similarity in composition among analyzed communities. For bacterial beta diversity assessment, Bray-Curtis dissimilarity index (Bray and Curtis, 1957) was calculated using QIIME 2 (Bolyen et al., 2019), with statistical significance tested using Past 4 version 4.08 (Hammer et al., 2001). Beta dispersion (variance) was determined utilizing the betadisper function from the Vegan package implemented in R v.4.3.1 (R Core Team, 2023) and performed using RStudio (RStudio Team, 2020).

Unique and shared taxa among the conditions were visualized using Venn diagrams created with an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). To quantify the relative abundance of bacterial taxa across different datasets, centered log-ratio (clr) was computed using the ANOVA-Like Differential Expression (ALDEX2) package (Fernandes et al., 2013), implemented in R v.4.3.1 (R Core Team, 2023) and performed using RStudio (RStudio Team, 2020). Subsequently, the resulting data were utilized to create a heatmap employing the heatmap.2 function implemented in R v.4.3.1 (R Core Team, 2023) and performed using the RStudio (RStudio Team, 2020). Only taxa with significant differences ($p \le 0.05$) were used in the

representation of the differential taxa relative abundance.

## 2.6. Inference of bacterial co-occurrence networks

We constructed co-occurrence networks for various datasets based on taxonomic profiles at the family and genus levels. Networks provide a graphical representation of microbial community assemblies, with nodes representing bacterial taxa and edges denoting significant correlations between taxa. To determine correlation strength we employed the Sparse Correlations for Compositional data (SparCC) method (Friedman and Alm, 2012) in R v.4.3.1 (R Core Team, 2023) and performed using RStudio environment (RStudio Team, 2020). Taxonomic data tables were used to calculate the correlation matrix. In these networks, each node represented a taxon and edges indicated either positive (average weight > 0.50) or negative (average weight < 0.50) correlations between the nodes. The colors of nodes were assigned based on modularity class metric values, and the node size was proportional to the eigenvector centrality of each taxon. The colors of the edges represent positive (blue) or negative correlations (red). The equivalences analysis of the modules was based on the modularity class.

Various network topological features were computed and visualized using Gephi 0.9.5 (Bastian et al., 2009). These features include the number of nodes and edges, network diameter (the shortest path between the two most separated nodes), modularity (indicating the strength of network division into modules), average degree (the average number of edges per node), weighted degree (the sum of edge weights connected to a node), and clustering coefficient (indicating the tendency of nodes to form clusters)(Abuin-Denis et al., 2024),

## 2.7. Local connectivity of *A. marginale* and CLE in the microbial community

To investigate the interactivity of *A. marginale* and CLE within the community, we determined their direct relationships with other bacterial microbiome members. For this purpose, sub-networks were constructed to visualize direct positive and negative associations. The analyses were conducted in Gephi 0.9.5 (Bastian et al., 2009), with the strength of the edges represented using the SparCC weight.

## 2.8. Network comparisons

The network comparisons were conducted using the package Network Construction and Comparison for Microbiome Data (NetCoMi) (Peschel et al., 2021) implemented in R v.4.3.1 (R Core Team, 2023) and performed using RStudio environment (RStudio Team, 2020). This analysis involved constructing and comparing networks with different strengths, categorized as average (SparCC weight > 0.5 / < -0.5) and strong average (SparCC weight > 0.9 / < -0.9) for microbiome data. To assess the similarity in the distribution of local centrality measures across nodes, including degree, betweenness centrality, closeness centrality, and eigenvector centrality, the Jaccard index was calculated for each centrality measure. The Jaccard index measures the similarity between sets of "most central nodes" (nodes with centrality values above

<footnote>* Forward (F), **Reverse (R)</footnote>
the empirical 75 % quartile) in two given networks. It expresses the similarity of both the sets of most central nodes and the sets of hub taxa between the two networks.

Additionally, the Adjusted Rand index (ARI) was calculated to evaluate the dissimilarity of clustering in the networks. ARI values range from -1-1, where negative values indicate lower-than-random clustering, positive values indicate higher-than-random clustering, 1 corresponds to identical clustering, and 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 (https://github.com/ramellose/anuran) implemented in Anaconda Python 3.10 environment. The CAN analysis identifies conserved patterns across networks using a core prevalence threshold 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. Finally, the visualization of the CAN was carried out using Gephi 0.9.5 (Bastian et al., 2009), a tool for exploring and analysis networks.

## 2.9. Keystone taxa identification

Keystone taxa were identified based on three criteria, as outlined in previous studies (Mateos-Hernández et al., 2020, 2021; Wu-Chuang et al., 2023). Firstly, these taxa had to be ubiquitous, signifying that were present across all samples within a particular group. Secondly, their eigenvector centrality had to be higher than 0.75. Eigenvector centrality is a metric for measuring the influence of a node within a network. Lastly, the mean relative abundance of these taxa needed to exceed that of the mean relative abundance of all taxa in the experimental group. To carry out this analysis, eigenvector centrality values were extracted from Gephi 0.9.5 software (Bastian et al., 2009).

## 2.10. Network robustness analysis

To assess the robustness of microbial co-occurrence networks to perturbation, we examined the impact of node removal or addition on network connectivity. The proportion of nodes removal required to achieve a connectivity loss of 0.80 for each network was evaluated. This assessment included random or directed attacks based on betweenness centrality (removing nodes with the highest betweenness centrality values), degree (removing nodes with the highest degree centrality values), and cascading (removing nodes with the highest betweenness centrality values, recalculating betweenness centrality after each removal). The Network Strengths and Weaknesses Analysis (NetSwan) package (Lhomme, 2015) in R v.4.3.1 (R Core Team, 2023) and performed using RStudio environment (RStudio Team, 2020) were employed for the network robustness analysis.

Additionally, a node addition analysis in RStudio was conducted. In each simulation, new nodes were randomly selected and connected to the existing network. The size of the largest connected component (LCC) and the average path length (APL) were calculated. To better estimate the network's robustness, we repeated the simulation with different sets of nodes (100, 300, 500, 700, 1000 nodes added). Finally, the obtained values were graphically represented using GraphPad 9 Prism software (GraphPad Software Inc., San Diego, CA, USA).

## 2.11. Statistical analysis

Differences in alpha diversity between groups were assessed using Kruskal-Wallis test (*p* < 0.05) in QIIME2 (Bolyen et al., 2019). Bray-Curtis dissimilarity index values were compared between groups using a PERMANOVA test (*p* < 0.05). Beta dispersion statistical analyses were performed using ANOVA test (*p* < 0.05). The taxa abundance differences were tested using Kruskal-Wallis test (*p* < 0.01) performed with the ALDEX2 package (Fernandes et al., 2013). For testing similarity in the distribution of local centrality measures between two networks,

two *p*-values P (J ≤ j) and P (J ≥ j) were computed for each local centrality measure (Real and Vargas, 1996). Jaccard index values of 0 and 1 represent minimum and maximum similarity, respectively (Real and Vargas, 1996). Differences were considered significant when *p* < 0.05. The standard error for loss of connectivity was calculated, considering variability, using a threshold of 0.9. The node addition analysis performs Wilcoxon signed-rank tests to determine if the mean size of the LCC and the APL are significantly different from 0. The *p*-values from these tests are adjusted using the Benjamini-Hochberg (BH) procedure to control for multiple comparisons. The removal nodes analysis use bootstrapping to calculate confidence intervals.

# 3. Results

## 3.1. Differential occurrence and abundance of Anaplasma marginale and CLE in the microbiome of different Rhipicephalus microplus developmental stages

We evaluated variations in bacterial diversity across different stages of R. microplus using three metrics of alpha diversity. Our findings showed that larvae have a more even bacterial distribution than adults (Kruskal-Wallis, *p* = 0.022, Fig. 1A) and nymphs (Kruskal-Wallis, *p* = 0.009, Fig. 1A). Additionally, nymphs displayed higher levels of Observed features compared to both larvae (Kruskal-Wallis, *p* = 0.001, Fig. 1A) and adults (Kruskal-Wallis, *p* < 0.018, Fig. 1A). Concerning Faith's index, the nymphs displayed higher values compared with larvae (Kruskal-Wallis, *p* = 0.0029, Fig. 1A) and adults (Kruskal-Wallis, *p* = 0.0018, Fig. 1A). However, when comparing larvae and adults, no significant differences were found in Observed features or Faith's index (Kruskal-Wallis, *p* > 0.05, Fig. 1A). The Bray-Curtis index comparisons indicated significant differences in bacterial compositions between larvae and adults (PERMANOVA, *p* = 0.001), nymphs and adults (PERMANOVA, *p* = 0.001) and larvae and nymphs (PERMANOVA, *p* = 0.002) (Fig. 1B). Significant differences were found also in the beta dispersion (ANOVA test, *p* = 0.001) suggesting high within-group variability. Overall, the bacterial microbiome diversity was markedly different between nymph and the other tick stages.

The compositional analysis revealed a total of 1038 taxa in the three stages, with 725 taxa (69.8 %) shared across all stages (Fig. 1C, Table S1). Unique taxa were identified in larvae (1.4 %, 15/1038), nymphs (5.2 %, 55/1038), and adults (4.0 %, 41/1038) (Fig. 1C, Table S1). Notably, *Coxiella* and *Anaplasma* were among the taxa shared by all stages, being present in at least one sample of each group. Differential relative abundance analysis revealed significant changes in 60 taxa (Fig. 1D). Among those, *Coxiella* exhibited higher abundance than other taxa across all stages, whereas *Anaplasma*, was highly abundant only in nymphs.

The presence of *A. marginale* and *Coxiella*-like endosymbionts (CLE) was assessed in all tick samples using *msp1* and *msp5* genes, and 16 S rRNA, respectively. CLE was present in 100 % (43/43) of the *R. microplus* samples, whereas *A. marginale* was only present in 44 % (19/43) of samples. The frequency of *A. marginale* occurrence differed between stages, being 31 % (4/13) in larvae, 80 % (12/15) in nymphs, and 20 % (3/15) in adults. The contrasting occurrence of *A. marginale* and CLE and relative abundance suggest a differential role of pathogens and symbionts on the tick bacterial microbiome.

## 3.2. Nested patterns of *A. marginale* and CLE in the *R. microplus* microbiome across stages

Bacterial co-occurrence networks (SparCC, weight ≥ 0.5 or ≤ −0.5) were constructed to assess the assembly of the microbiome in *R. microplus* ticks and to trace variations in the nested patterns of *A. marginale* and CLE across tick developmental stages. Visual inspection revealed changes in network topology in larval (Fig. 2A), nymphal (Fig. 2B) and adult (Fig. 2C) stages. There was a reduction in the number
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Comparison of microbial diversity and taxa composition across the stages. (a) Comparison of alpha diversity with evenness, Observed features index and Faith PD for larvae (blue), nymphs (red) and adults (brown). (b) Comparison of beta diversity with Bray Curtis dissimilarity index for larvae (blue), nymphs (red) and adults (brown) samples. Small circles and triangles represent samples; ellipses represent centroid position for each group ($F= 77.2; p = 0.0001^{***}$; stress = 0.562). (c) Venn diagram displaying the comparison of taxa composition for all taxa. Numbers represent the taxa found on each dataset and those shared by the groups. (d) Comparison of relative abundance of complex microbial communities within stages. Each column represents the clr values for bacterial taxa per sample and per group. Each line represents bacterial taxa with significant changes between the datasets. Colors represent the clr value range from -10 (blue) to 10 (red).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Dynamics of microbial communities of *R. microplus* across stages. Co-occurrence network of (a) larvae, (b) nymphal and (c) adult stages. Node colors are based on modularity class metric and equal color means modules of co-occurring taxa. The colors in the edges represent strong positive (blue) or negative (red) correlations (SparCC > 0.50 or < -0.50). Black dots represent *Anaplasma* taxon and white dots *Coxiella* taxon. Larvae (L1–3), nymphs (N1–3) and adults (A1–3) indicated the three major modules in the network.</figcaption></figure>

of communities as stages progressed (Table 2). Concurrently, there was an increase in network complexity, as evidenced by an increase in the number of edges, network diameter, average degree and weighted degree across consecutive stages (Table 2). Notably, positive interactions

between nodes were predominant in the networks (Table 2). Overall complexity increased as the developmental stages progress, with the nymphal stage exhibiting a particular jump in complexity compared to other stages (Table 2).
<table><caption>**Table 2**
Topological features of the taxonomic networks of tick stages.</caption><thead><tr><th>Topological Features</th><th>Larvae</th><th>Nymphs</th><th>Adults</th></tr></thead><tbody><tr><td>Nodes</td><td>566</td><td>604</td><td>545</td></tr><tr><td>Edges</td><td>3150</td><td>7879</td><td>8769</td></tr><tr><td>Positives</td><td>1936 (61.46 %)</td><td>4479 (56.84 %)</td><td>5315 (60 %)</td></tr><tr><td>Negatives</td><td>1214 (38.54 %)</td><td>3400 (43.16 %)</td><td>3544 (40 %)</td></tr><tr><td>Modularity</td><td>1.773</td><td>1.988</td><td>1.435</td></tr><tr><td>Number of Communities</td><td>83</td><td>80</td><td>62</td></tr><tr><td>Network Diameter</td><td>7</td><td>8</td><td>9</td></tr><tr><td>Average Degree</td><td>11.13</td><td>26.08</td><td>32.57</td></tr><tr><td>Weighted Degree</td><td>1.59</td><td>2.594</td><td>4.707</td></tr></tbody></table>

Comparison of nodes within the three major modules—L1, L2, L3 for larvae; N1, N2, N3 for nymphs; and A1, A2, A3 for adults (Fig. 2)—across the three networks revealed substantial sharing of taxa among different developmental stages (Table 3). Taxa were shared significantly within these primary modules across larvae, nymphs, and adults (Table 3). Notably, in all networks CLE and *A. marginale* were placed into distinct modules (Fig. 2). In the larval network, *A. marginale* was present in the thirst largest (L3) community, while CLE was associated with a smaller community (Fig. 2A, Table S2). In the nymphal network, CLE was found in the largest (N1) community whereas *A. marginale* was in the second largest (N2) community (Fig. 2B, Table S2). In adult network, CLE prominent in the A1 community, while *A. marginale* was identified in smaller communities (Fig. 2C, Table S2). This pattern indicated a more stable association of CLE with its community across stages, whereas *A. marginale*'s association with the microbial communities varied.

We then checked whether CLE and *A. marginale* were among taxa with the strongest correlations within the networks (SparCC, weight ≥ 0.9 or ≤ − 0.9). In larval and nymphal networks 29 nodes displayed strong correlations (Fig. 3A). In nymphal and adult networks 59 nodes displayed strong correlation (Fig. 3B) and between larvae and adults, 44 nodes displayed strong correlations (Fig. 3C). In nymphal network, *A. marginale* taxa display a strong correlation with *Pseudoalteromonas*, identified as a keystone taxon (Fig. 3, Table 4). In contrast, CLE did not exhibit stable correlations conserved across stages.

Changes in the nodes interacting with *A. marginale* and CLE were observed across different developmental stages (Fig. 4, Table S3). In the larval stage, *A. marginale* exhibited fewer connections compared to CLE (Fig. 4, Fig. S1, Table S3). In the nymphal network, both *A. marginale* and CLE displayed an increase in associations compared to the larval stage (Fig. 4, Fig. S1, Table S3). In the adult network, *A. marginale* connections decreased while CLE maintained a high level of connectivity (Fig. 4, Fig. S1, Table S3). Notably, no taxon was directly linked to both *A. marginale*, and CLE in the larval network, whereas in the nymphs and adult network, 77 taxa and 2 taxa, respectively, were connected to both *A. marginale* and CLE (Fig. 4, Tables S4-S5). Interestingly, in nymphs, *A. marginale* and CLE showed a direct negative association between them and were directly connected with keystone taxa (Table 4). Specifically, *A. marginale* was found to interact positively with the seven

keystone taxa in nymph, while CLE interacted negatively with six of them (Table 4, Table S3).

From larval to nymphal network, *A. marginale* maintained connections with taxa such as *Atopostipes*, *Enhydrobacter* and *Candidatus Competitibacter RF39*, while CLE was connected to *Lautropia*, *Comamnadaceae*, *Mycoplasma*, *Mogibacterium* and *Porphyromonas* (Fig. 4, Table S4). From nymphs to adult, *A. marginale* remains connected only to *Rikenellaceae RC9_gut_group*, while CLE maintains the same connections with 13 taxa (Fig. 4, Table S5). No taxa were consistently connected to *A. marginale* across all stages, whereas CLE was consistently connected to *Bacillus* (Table S4).

### 3.3. Impact of *A. marginale* and CLE removal on network structure and robustness

To assess the significance of *A. marginale* and CLE in the networks, we removed *in silico* these nodes from the network and then analyzed the changes in network structure. Removal of the node representing *A. marginale* resulted in a decrease in the number of nodes in the larval and adult network, while in the nymphal network the number of nodes increased (Fig. S2, Table 5). The networks without the *A. marginale* node showed reduced number of edges in larvae and nymphs but increased number of edges in adults. A notable effect of *A. marginale* removal was the increase in community numbers in larval and nymphal networks, contrasting with a decrease in adult network. On the other hand, removing CLE had a subtle effect on the overall network structure but caused a greater reduction in the number of connected nodes than the removal of *A. marginale* (Table 5).

Further analysis on the modular composition after *A. marginale* removal revealed that CLE consistently appeared in the largest communities (L1, N1, A1) across all stages. However, without CLE, *A. marginale* was consistently found in the second-largest modules (L2, N2, A2), highlighting its persistent presence across different stages (Fig. S2). This suggests a reciprocal relationship between *A. marginale* and CLE, wherein the removal of one taxon influences the positioning of the other within the modular structure.

Network resilience to perturbations, including node removal and addition, was analyzed across developmental stages. Directed attacks led to faster connectivity loss than random attacks (Fig. 5A), with cascading removal having the most significant impact, followed by betweenness and degree centrality (Fig. 5B, Table 6). The proportion of nodes needed to induce an 80 % loss of connectivity increases as the stages progress, indicating an enhancement in robustness, except for cascading effects in nymphal network (Table 6). Overall, adult network showed greater resistance compared to larval and nymphal network, requiring the removal of a larger fraction of nodes for an 80 % connectivity loss (Fig. 5A, Table 6).

When examining the networks' response to the addition of nodes, distinct trends were observed. Larval network maintained a stable size for the largest connected component (LCC), yet exhibited an increase after the addition of 500 nodes. Nymphal network displayed an increased in LCC size with the addition of nodes, while adult network

<table><caption>**Table 3**
The percentage of taxa shared between different modules.</caption><thead><tr><th></th><th>L1</th><th>L2</th><th>L3</th><th>N1</th><th>N2</th><th>N3</th><th>A1</th><th>A2</th><th>A3</th></tr></thead><tbody><tr><td>L1</td><td>100 %</td><td>0 %</td><td>0 %</td><td>21.72 %</td><td>8.98 %</td><td>2.70 %</td><td>13.24 %</td><td>20.16 %</td><td>0 %</td></tr><tr><td>L2</td><td></td><td>100 %</td><td>0 %</td><td>11.14 %</td><td>18.59 %</td><td>5.08 %</td><td>19.12 %</td><td>10.58 %</td><td>2.32 %</td></tr><tr><td>L3</td><td></td><td></td><td>100 %</td><td>4.33 %</td><td>6.35 %</td><td>21.95 %</td><td>3.23 %</td><td>9.33 %</td><td>10.34 %</td></tr><tr><td>N1</td><td></td><td></td><td></td><td>100 %</td><td>0 %</td><td>0 %</td><td>24.35 %</td><td>21.13 %</td><td>0.78 %</td></tr><tr><td>N2</td><td></td><td></td><td></td><td></td><td>100 %</td><td>0 %</td><td>18.07 %</td><td>13.86 %</td><td>0.51 %</td></tr><tr><td>N3</td><td></td><td></td><td></td><td></td><td></td><td>100 %</td><td>4.62 %</td><td>6.08 %</td><td>10.58 %</td></tr><tr><td>A1</td><td></td><td></td><td></td><td></td><td></td><td></td><td>100 %</td><td>0 %</td><td>0 %</td></tr><tr><td>A2</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>100 %</td><td>0 %</td></tr><tr><td>A3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>100 %</td></tr></tbody></table>

L: larvae, N: nymphs, A: adults
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Strong correlations in tick stage networks. Strong correlation networks (SparCC ≥ 0.90) across tick stages: (a) larvae vs nymphs and (b) nymphs vs adults (c) larvae vs adults networks. Nodes correspond to taxa (family or genera level) and connecting edges indicate significant connection between them. Node colors represent different modules in the networks. Edge colors represent strong positive (green) or negative (red) correlations.</figcaption></figure>
<table><caption>**Table 4**

Keystone taxa in bacterial communities across *R. microplus* stages.</caption><thead><tr><th>Stages</th><th>Keystone taxa</th></tr></thead><tbody><tr><td>Larva</td><td>None found</td></tr><tr><td>Nymph</td><td>Ralstonia**<br>Pseudoalteromonas**<br>Vibrio**<br>Blautia**<br>Cupriavidus**<br>Enterobacteriaceae*<br>Comamonadaceae**</td></tr><tr><td>Adult</td><td>Clostridia UCG-014<br>Clostridium sensu stricto 1</td></tr></tbody></table>

reached a plateau after the addition of 300 nodes. Notably, nymphal network exhibited the highest LCC values, suggesting their ability to maintain the connectivity of a network despite additional nodes (Fig. 5C). Concerning the average path length (APL) larval network exhibited the highest values, followed by nymphal and adult network (Fig. 5D). After the addition of 300 nodes, the three networks exhibited a decrease in the APL values. Larval network reached a plateau after the 300 nodes, whereas the values for nymphal and adult network continued to decrease (Fig. 5D). Adult network showed lower APL, indicating enhanced efficiency in connections and interactions amid disturbances (Fig. 5D).

The networks robustness was also compared after *in silico* removal of *A. marginale* and CLE. In the larval network, a lower fraction of nodes was required to reach an 80 % loss of connectivity when either *A. marginale* or CLE was removed (Fig. S3, Table 7). For the nymphal network, the fraction of nodes needed to reach an 80 % loss of

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Local co-occurrences of *Anaplasma* and *Coxiella*. The direct neighbors of *Anaplasma* and *Coxiella* were identified in the bacterial co-occurrence networks of larvae, nymphs and adults. Black dots represent *Anaplasma* taxon and white dots *Coxiella* taxon. Green dots represent unique taxa connected to *Anaplasma* or *Coxiella*. Red dots represent shared taxa between *Coxiella* and *Anaplasma*. Purples and blues dots represent taxa shared between the stages connected to *Coxiella* or *Anaplasma*, respectively. Blue lines represent positive interactions and red lines represent negative interactions. Numbers inside dots represent the amounts of taxa unique or shared.</figcaption></figure>

<table><caption>**Table 5**

Topological features of the taxonomic networks of tick stages after *in silico* removal of *A. marginale* and CLE.</caption><thead><tr><th rowspan="2">Network features</th><th colspan="3">Network woAna</th><th colspan="3">Network woCox</th></tr><tr><th>Larvae</th><th>Nymphs</th><th>Adults</th><th>Larvae</th><th>Nymphs</th><th>Adults</th></tr></thead><tbody><tr><td>Nodes</td><td>557</td><td>612</td><td>539</td><td>567</td><td>601</td><td>543</td></tr><tr><td>Edges</td><td>3142</td><td>7599</td><td>8884</td><td>3112</td><td>7774</td><td>8785</td></tr><tr><td>Modularity</td><td>1.889</td><td>2.026</td><td>1.441</td><td>1.85</td><td>2.003</td><td>1.426</td></tr><tr><td>Number of Communities</td><td>97</td><td>81</td><td>63</td><td>74</td><td>70</td><td>70</td></tr><tr><td>Network Diameter</td><td>7</td><td>9</td><td>8</td><td>7</td><td>7</td><td>8</td></tr><tr><td>Average Degree</td><td>11.28</td><td>24.83</td><td>32.965</td><td>10.99</td><td>25.87</td><td>32.35</td></tr><tr><td>Weighted Degree</td><td>1.506</td><td>2.526</td><td>4.821</td><td>1.52</td><td>2.72</td><td>4.656</td></tr><tr><td>Triangles number</td><td>7904</td><td>88093</td><td>143276</td><td>7709</td><td>92802</td><td>143061</td></tr></tbody></table>

<footnote>* Interacting only with *A. marginale* ** Interacting with *A. marginale* and CLE</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Comparison of network robustness after removal and addition of nodes. (a)Connectivity loss measured after directed attack, (betweenness, cascading, degree) or random attack for larval, nymph and adult networks. (b). Comparison of connectivity loss measured for the three stages in cascading, betweenness and degree. Comparison of the value reached of (c) LCC size and the (d) APL for larval (blue), nymph (red) and adult (brown) networks, after each node addition was plotted.</figcaption></figure>

<table><caption>**Table 6**
The fraction of nodes removed to reach an 80 % of connectivity loss.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="3">Fraction of nodes removed</th></tr><tr><th>Larvae</th><th>Nymphs</th><th>Adults</th></tr></thead><tbody><tr><td>Cascading</td><td>0.23</td><td>0.28</td><td>0.23</td></tr><tr><td>Betweenness</td><td>0.28</td><td>0.32</td><td>0.36</td></tr><tr><td>Degree</td><td>0.33</td><td>0.38</td><td>0.42</td></tr><tr><td>Random</td><td>0.49</td><td>0.49</td><td>0.50</td></tr></tbody></table>

connectivity was lower based on betweenness and degree centrality after the removal of *A. marginale*. Conversely, after the removal of CLE, cascading, degree, and betweenness centrality exhibited lower values, indicating less impact when CLE was present. In the adult network, the removal of *A. marginale* resulted in a lower fraction of nodes needed to reach an 80 % loss of connectivity for betweenness and degree centrality

<table><caption>**Table 7**
The fraction of nodes removed to reach 80 % of connectivity loss after *in silico* removal of *A. marginale* and CLE.</caption><thead><tr><th rowspan="3">Parameters</th><th colspan="6">Fraction of nodes removed</th></tr><tr><th colspan="3">woAna</th><th colspan="3">woCox</th></tr><tr><th>Larvae</th><th>Nymphs</th><th>Adults</th><th>Larvae</th><th>Nymphs</th><th>Adults</th></tr></thead><tbody><tr><td>Cascading</td><td>0.22</td><td>0.28</td><td>0.24</td><td>0.22</td><td>0.29</td><td>0.22</td></tr><tr><td>Betweenness</td><td>0.27</td><td>0.31</td><td>0.35</td><td>0.27</td><td>0.34</td><td>0.36</td></tr><tr><td>Degree</td><td>0.31</td><td>0.36</td><td>0.41</td><td>0.32</td><td>0.39</td><td>0.42</td></tr><tr><td>Random</td><td>0.49</td><td>0.49</td><td>0.50</td><td>0.48</td><td>0.50</td><td>0.49</td></tr></tbody></table>

compared to when *A. marginale* was present. However, the removal of CLE only impacted the cascading values of the adult network (Fig. S3, Table 7). After simulating the removal of *A. marginale* and CLE, the nymphal network displayed the highest robustness in terms of LCC,
followed by larvae and then adults. Conversely, for APL values, adult network demonstrated higher robustness, followed by nymphs, with larvae being the least robust (Fig. S4).

### 3.4. Impact of *A. marginale* and CLE removal on core features of microbial communities

The Core Association Network (CAN) analysis revealed 42 nodes consistently exhibiting 34 positive correlations across all three tick developmental stages, with only one negative interaction detected (between *Lawsonella* with *Alistipes*) (Fig. 6A). Neither *A. marginale* nor CLE exhibit a consistent partner of correlation throughout the stages. When comparing larval and nymphal stages exclusively, a CAN comprising 101 nodes with 21 edges representing negative correlations was identified. Notably, *A. marginale* exhibited negative interaction with (*Eubacterium*) coprostanoligenes group and positive interaction with *Enhydrobacter*. CLE displayed similar correlation patterns in the CAN resulting from larval and nymphal stages, positively connected with *Bacillus* (Fig. 6B). The CAN of nymphal and adult networks revealed 122 nodes with 166 edges, 41 of which represented negative interactions (Fig. 6C). *Anaplasma marginale* was negatively connected to *Rikenella*-ceae, while CLE was not part of the CAN.

Upon *in silico* removal of *A. marginale*, an increase in the number of nodes (121 nodes) and 163 edges was observed in the CAN including all tick stages. Conversely, the removal of CLE resulted in a CAN with a decreased number of nodes (23 nodes) and edges (14 edges) (Fig. S5A). The number of nodes and edges in the CAN of larval and nymphal networks increased to 228 nodes and 515 edges, and 221 nodes and 486 edges, following the *in silico* removal of *A. marginale* and CLE, respectively (Fig. S5B). The CAN of nymphal and adult networks had lower number of nodes (68 nodes) and edges (67 edges) after *in silico* removal of CLE, while a larger set of 228 nodes and 584 edges were found in the CAN following the removal of *A. marginale* (Fig. S5C).

### 3.5. Impact of CLE and *A. marginale* removal on node centrality distribution

The similarity in distribution of local centrality measures across networks was tested using the Jaccard index comparison. Jaccard index values for degree, hub taxa, betweenness, and eigenvector centrality between larval and nymphal network were significantly higher than expected by random (P(≥Jacc), *p* < 0.01, Table 8), whereas eigenvector centrality was significantly lower than expected by random (P(≥Jacc), *p* < 0.05, Table 8). Similarly, when we compared nymph and adult networks, Jaccard index values for degree, hub taxa, closeness and

eigenvector centrality were significantly higher than expected by random (P (≥Jacc), *p* < 0.05, Table 8), while the Jaccard value for betweenness centrality did not deviate from random distribution (P (≥Jacc), *p* > 0.05, Table 8). This result suggests a remarkable equivalence in the node centrality distribution across these networks. The ARI similarity index revealed the value expected for two random clustering: larval vs nymphal network (ARI = 0.165, *p* = 0.001) and nymphs vs adult networks (ARI = 0.162, *p* = 0.001).

In the context of *A. marginale* and CLE removal, the Jaccard indexes approached 1 for all local centrality measures, emphasizing a striking similarity in the proportion of central nodes within networks without *A. marginale* or CLE across tick stages, with eigenvector centrality exhibiting slight variation (Table 9). A higher similarity was observed between nymph and adult stages after *A. marginale* or CLE removal. The comparison between the network without *A. marginale* or CLE showed a higher ARI value than in networks with *A. marginale* and CLE, which is indicative of a more similar clustering pattern (Table 10).

## 4. Discussion

Our research aimed to explore the distinct ecological roles of *Coxiella*-like endosymbiont (CLE) and *Anaplasma marginale* within the *Rhipicephalus microplus* microbiome, hypothesizing that the mutualistic CLE would demonstrate a stable, nested pattern reflecting co-evolution (Bintruy et al., 2020; Zhong, 2012), whereas the pathogenic *A. marginale* would exhibit a more dynamic (Piloto-Sardinas et al., 2024), non-nested pattern influenced by environmental factors (Cabezas-Cruz et al., 2013; Estrada-Peña et al., 2009) and the immune response of the host (Kenneil et al., 2013; Kocan et al., 2004) and the tick vector (Kalil et al., 2017). The findings from this study lead us to partially accept this hypothesis. While the evidence supports the stable, nested presence of CLE, suggesting a co-evolved mutualism, the data on *A. marginale* indicate a more complex interaction within the microbiome than a purely non-nested, dynamic pattern, as it still plays a significant role in structuring the microbial community.

The higher prevalence of *Anaplasma* in nymphs suggests that the nymphal stage could play a more significant role in pathogen transmission. A similar effect was observed on ticks experimentally infected with *Borrelia afzelii* (Hamilton et al., 2021). Infection with *B. afzelii* appears to increase bacterial microbiome diversity in nymphs by having disproportionately impacts on abundant OTUs, leading to increased community evenness and Shannon diversity index (Hamilton et al., 2021). Furthermore, the identification of keystone taxa, such as *Pseudolateromonas* in nymphs, highlighted their potential roles in modulating the microbiome and influencing pathogen presence. This

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Core Association Networks in tick ontogeny. Taxa shared the same conserved partner throughout (a) the tree stages, (b) larval-nymph and (c) nymph-adult. The nodes represent bacterial taxa. Node colors are based on modularity class metric and edges colors represent strong positive (blue) or negative (red) correlations.</figcaption></figure>
<table><caption>**Table 8**
Jaccard index for larval, nymphal and adult networks.</caption><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">Larvae vs. Nymphs</th><th colspan="3">Nymphs vs. Adults</th></tr><tr><th>Jacc</th><th>P(≤Jacc)</th><th>P(≥Jacc)</th><th>Jacc</th><th>P(≤Jacc)</th><th>P(≥Jacc)</th></tr></thead><tbody><tr><td>Degree centrality</td><td>0.394</td><td>0.986797</td><td>0.017974 *</td><td>0.172</td><td>0.983900</td><td>0.045131 *</td></tr><tr><td>Betweenness centrality</td><td>0.502</td><td>1.000000</td><td>0.000000 **</td><td>0.125</td><td>0.960982</td><td>0.195092</td></tr><tr><td>Closeness centrality</td><td>0.413</td><td>0.997897</td><td>0.003065 **</td><td>0.172</td><td>0.983900</td><td>0.045131 *</td></tr><tr><td>Eigenvector centrality</td><td>0.287</td><td>0.043651 *</td><td>0.966579</td><td>0.172</td><td>0.983900</td><td>0.045131 *</td></tr><tr><td>Hub taxa</td><td>0.551</td><td>1.000000</td><td>0.000000 ***</td><td>0.172</td><td>0.983900</td><td>0.045131 *</td></tr></tbody></table>

Significance codes: **\*: < 0.001, **: < 0.01, *: < 0.05.

<table><caption>**Table 9**
Jaccard index for tick stages networks after *in silico* removal of *A. marginale* or CLE.</caption><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">Larvae (woAna) vs Nymphs (woAna)</th><th colspan="3">Nymphs (woAna) vs Adults (woAna)</th></tr><tr><th>Jacc</th><th>P(≤Jacc)</th><th>P(≥Jacc)</th><th>Jacc</th><th>P(≤Jacc)</th><th>P(≥Jacc)</th></tr></thead><tbody><tr><td>Degree centrality</td><td>0.472</td><td>1.000000</td><td>0.000000 ***</td><td>0.463</td><td>1.000000</td><td>0.000000 ***</td></tr><tr><td>Betweenness centrality</td><td>0.584</td><td>1.000000</td><td>0.000000 ***</td><td>0.520</td><td>1.000000</td><td>0.000000 ***</td></tr><tr><td>Closeness centrality</td><td>0.507</td><td>1.000000</td><td>0.000000 ***</td><td>0.499</td><td>1.000000</td><td>0.000000 ***</td></tr><tr><td>Eigenvector centrality</td><td>0.383</td><td>0.982752</td><td>0.02227 *</td><td>0.387</td><td>0.988003</td><td>0.015699 *</td></tr><tr><td>Hub taxa</td><td>0.639</td><td>1.000000</td><td>0.000000 ***</td><td>0.599</td><td>1.000000</td><td>0.000000 ***</td></tr><tr><td></td><th colspan="3">Larvae (woCox) vs Nymphs (woCox)</th><th colspan="3">Nymphs (woCox) vs Adults (woCox)</th></tr><tr><td>Degree centrality</td><td>0.496</td><td>1.000000</td><td>0.000000 ***</td><td>0.471</td><td>1.000000</td><td>0.000000 ***</td></tr><tr><td>Betweenness centrality</td><td>0.558</td><td>1.000000</td><td>0.000000 ***</td><td>0.534</td><td>1.000000</td><td>0.000000 ***</td></tr><tr><td>Closeness centrality</td><td>0.499</td><td>1.000000</td><td>0.000000 ***</td><td>0.494</td><td>1.000000</td><td>0.000000 ***</td></tr><tr><td>Eigenvector centrality</td><td>0.387</td><td>0.988003</td><td>0.015699 *</td><td>0.394</td><td>0.994518</td><td>0.007373 **</td></tr><tr><td>Hub taxa</td><td>0.640</td><td>1.000000</td><td>0.000000 ***</td><td>0.608</td><td>1.000000</td><td>0.000000 ***</td></tr></tbody></table>

Significance codes: **\*: < 0.001, **: < 0.01, *: < 0.05.

<table><caption>**Table 10**
Network clustering comparisons.</caption><thead><tr><th>Conditions</th><th>Network comparisons</th><th>ARI</th><th>p-value</th></tr></thead><tbody><tr><td rowspan="2">Without Anaplasma</td><td>Larvae vs Nymphs</td><td>0.315</td><td>0</td></tr><tr><td>Nymphs vs. Adults</td><td>0.308</td><td>0</td></tr><tr><td rowspan="2">Without Coxiella</td><td>Larvae vs Nymphs</td><td>0.346</td><td>0</td></tr><tr><td>Nymphs vs Adults</td><td>0.313</td><td>0</td></tr></tbody></table>

undercores the potential importance of targeting nymphal stages in disease control strategies. Nymphs also exhibited significantly higher bacterial diversity than larvae and adults, which contrasts with previous results where tick eggs and larvae exhibited the highest microbial richness and diversity (Ruiling et al., 2019; Zhang et al., 2020). Over the course of tick development, bacterial diversity increased, with adult ticks demonstrating higher diversity than other ticks stages (Gomard et al., 2021). The increase in the microbial diversity of nymphs could be linked to the presence of *Anaplasma*, as pathogen infections have been shown to influence tick microbiome diversity and composition (Abraham et al., 2017; Maitre et al., 2023; Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022).

The ubiquitous presence of CLE across all developmental stages of *R. microplus*, detected in 100 % of samples, supports a stable, potentially co-evolved mutualistic relationship, which is consistent with a nested microbiome pattern that aligns with our hypothesis. Additionally, the substantial decrease in the number of connected nodes upon CLE's removal demonstrates that CLE is a central connector within the microbiome. Nodes that were connected through CLE likely rely on it either directly or indirectly. Nested structure of symbionts has been described in various systems (McCutcheon and von Dohlen, 2011), and involves mainly nutritional symbioses in obligate blood feeders (Duron and Gottlieb, 2020; McCutcheon and von Dohlen, 2011). CLE is known to provide essential B vitamins and cofactors critical for tick development, fecundity, and metabolism (Duron et al., 2017; Guizzo et al., 2017; Smith et al., 2015; Zhang et al., 2017; Zhong et al., 2007). Additionally, studies have demonstrated that reducing CLE levels

through antibiotic treatments in ticks like *R. microplus*, *Rhipicephalus sanguineus*, and *Haemaphysalis longicornis* significantly compromises their blood intake (Ben-Yosef et al., 2020; Guizzo et al., 2017; Zhang et al., 2017). This mutualistic symbiosis is crucial for the survival of both the microbe and its tick host, and therefore a stable nesting in the community is expected. This was further validated by bacterial co-occurrence networks that show CLE's consistent associations within significant microbial communities across tick stages, suggesting a co-evolved, nested pattern. These findings are consistent with previous research suggesting a close and stable relationship between ticks and endosymbionts that is maintained across various stages of tick development (Andreotti et al., 2011; Binetruy et al., 2020; Clayton et al., 2015; Lalzar et al., 2012).

Conversely, the variable occurrence of *A. marginale* across stages but its notable influence on network complexity (see below), especially in nymph stages, suggests a dynamic yet structurally impactful presence. Several factors can influence the variable presence of *A. marginale* in tick vectors. For example, tick immunity (Kalil et al., 2017; Mercado-Curiel et al., 2014), microbiome dynamics (Gall et al., 2016; Piloto-Sardiñas et al., 2024), and protective symbionts (Gall et al., 2016) can affect susceptibility of tick vectors to *A. marginale* infection. Moreover, a robust immune response in persistently infected cattle can influence tick fitness and limit the ticks' ability to acquire *A. marginale* (Kenneil et al., 2013; Kocan et al., 2004). Despite variations in *A. marginale* occurrence in ticks (Obregón et al., 2019; Piloto-Sardiñas et al., 2023; Cicculli et al., 2020; Scoles et al., 2007; Shimada et al., 2004; Xu et al., 2023), here we showed that this pathogen engages in structured microbial interactions, challenging the initial hypothesis of a non-nested, purely dynamic pattern. In contrast, a previous study noted strong temporal dynamics in *A. marginale*'s nesting patterns within the microbial communities of engorged adult female *R. microplus* ticks, with each time point across three years showing varying levels of modularity and connectivity of this pathogen with other member of the tick microbiome (Piloto-Sardiñas et al., 2024). This indicates that tick ontogeny and temporal variations of the microbiome can affect *A. marginale*'s nesting patterns within the tick microbiome.
We observed a direct negative interaction between CLE and *A. marginale* in nymphs, suggesting microbial competition. Such competitive interactions could significantly impact pathogen transmission dynamics. For example, Gall et al. (2016) examined the bacterial microbiome of two populations of *Dermacentor andersoni* with historically different susceptibilities to *A. marginale*. They found that ticks more susceptible to *A. marginale* harbored a higher proportion of *Arsenophonus* spp. and *Rickettsia peacockii*, while those less susceptible had a greater abundance of *Francisella*-like endosymbionts (FLEs), and *Rickettsia bellii* (Gall et al., 2016). Given that *R. bellii* was not detected in ticks with high *A. marginale* prevalence, the authors proposed that the *R. bellii* endosymbiont may negatively impact tick acquisition of *A. marginale* (Gall et al., 2016). In *Amblyomma maculatum*, *Coxiella* operational taxonomic units (OTUs), believe to be CLE, were also negatively correlated with *Rickettsia* (Adegoke et al., 2022). Similar negative associations were found between *Rickettsia helvetica* and *Candidatus* Midichloria mitochondrii and *Spiroplasma ixodetis* in *I. ricinus* ticks collected in the field (Krawczyk et al., 2022), as well as between *R. helvetica* and Enterobacteriaceae, Comamonadaceae, and Bacillus in *R. helvetica*-infected *I. ricinus* ticks collected from humans (Maitre, Wu-Chuang, Mateos-Hernández, et al., 2022).

Previous studies on the microbiome of *I. ricinus* a negative association was also found between *Rickettsia* and *Spiroplasma* (Aivelo et al., 2019, Lejal et al., 2021). The ubiquitous presence of CLE across all samples and the low infection rate of *A. marginale* in our study suggest that, influenced by biotic or abiotic factors, CLE could inhibit *A. marginale* infection in *R. microplus*. Beyond nutrient provisioning, mutualistic symbionts may confer protection against pathogens (Brownlie and Johnson, 2009; Smith and Ashby, 2023). Mechanisms of protection may involve immunomodulation, wherein protective symbionts prime the arthropod immune system to continually produce low quantities of antimicrobial peptides, thereby either inhibiting the multiplication of other microbes (Fogaça et al., 2021; Pavanelo et al., 2023). Additional mechanisms such as resource competition, direct interference (Brownlie and Johnson, 2009; Smith and Ashby, 2023), or priority effects (Maitre et al., 2023) could also account for the observed negative association between CLE and *A. marginale*.

We perturbed bacterial co-occurrence networks by in silico removal of *Anaplasma* and CLE, comparing topology, connectivity, centrality measures, robustness, and community assembly with original networks. *In silico* node removal has been used previously to assess the influence of microorganisms on tick (Maitre et al., 2023; Piloto-Sardiñas et al., 2024), and plant (Agler et al., 2016) microbiota properties. For example, *Arabidopsis thaliana* microbiota hub taxa removal affected more network edges than non-hub species, and the biological relevance of network-predicted hub taxa was confirmed by host colonization experiments (Agler et al., 2016), which validates *in silico* removal as a tool to predict ecosystem behaviour. Additionally, *Rickettsia* removal from tick networks notably impacted bacterial community assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* microbiota, altering centrality values and bacterial node presence (Maitre et al., 2023). In our study, the removal of *A. marginale* and CLE from the networks led to distinct changes in network structure across different developmental stages of the tick, highlighting their influence on microbiome complexity and connectivity. Specifically, the removal of *A. marginale* resulted in a decrease in node numbers in larval and adult networks but an increase in nymphal network. This suggests that *A. marginale* plays varying roles across stages, potentially suppressing some bacterial species particularly in nymphs where its removal increases network complexity. Conversely, the removal of CLE generally caused a greater reduction in the number of connected nodes across all stages, indicating its critical role in maintaining network integrity and connectivity. This underscores CLE's central role as a connector within the microbiome, facilitating interactions among various microbial taxa. Interestingly, the analysis showed that the modular composition of the networks was significantly altered upon the removal of these nodes. Without CLE, *A. marginale*

appeared in larger community clusters than when CLE was present, possibly indicating that CLE's presence suppresses the expansion or connectivity of *A. marginale* within the microbiome. This agrees with the negative interaction found between CLE and *A. marginale* and as mentioned above, it could be due to direct competition or modulation of the microbial environment by CLE that is unfavorable to *A. marginale*.

## 5. Conclusion

This study unravels distinct nested patterns of *A. marginale* and CLE in the bacterial communities of *R. microplus*. The observed differences in nested patterns could be indicative of the secondary nature of *A. marginale* acquisition during feeding, contrasting with the evolutionarily stable relationship between CLE and ticks. This suggests that the acquisition dynamics of *A. marginale*, likely occurring during feeding on infected hosts, could contribute to its transient and dynamic presence within the tick microbiome. Conversely, the consistently high abundance and stable connections of CLE across developmental stages emphasize a long-standing and evolutionarily stable association with the tick microbiome. Understanding these dynamics is crucial for developing targeted interventions. For instance, the stable relationship between ticks and CLE might offer targets for disrupting the viability of tick populations, thereby reducing their capacity as disease vectors.

## Funding statement

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

## CRediT authorship contribution statement

**Alina Rodríguez-Mallon:** Writing – review & editing, Supervision, Conceptualization. **Lourdes Mateos-Hernández:** Writing – review & editing, Data curation. **Mario Pablo Estrada-García:** Writing – review & editing. **Yamil Bello:** Writing – review & editing, Investigation. **Alejandro Cabezas-Cruz:** Writing – review & editing, Writing – original draft, Supervision, Conceptualization. **Patrícia Gonzaga Paulino:** Writing – review & editing, Formal analysis. **Anays Alvarez Gutierrez:** Writing – review & editing, Investigation. **Frank Ledesma Bravo:** Writing – review & editing, Investigation. **Alier Fuentes Castillo:** Writing – review & editing, Investigation. **Rafmary Rodríguez Fernández:** Writing – review & editing, Investigation. **Angélique Foucault-Simonin:** Writing – review & editing, Investigation. **Elianne Piloto-Sardiñas:** Writing – review & editing, Formal analysis. **Luis Méndez Mellor:** Writing – review & editing, Investigation. **Lianet Abuin-Denis:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Formal analysis, Conceptualization. **Alejandra Wu-Chuang:** Writing – review & editing, Supervision, Methodology. **Dasiel Obregon:** Writing – review & editing, Software. **Apolline Maitre:** Writing – review & editing, Supervision, Methodology.

## Declaration of Competing Interest

The authors declare no conflict of interest.

## Data Availability

Raw data used in the study was uploaded to public repository, SRA: Bioproject PRJNA1034019
# Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi.org/10.1016/j.micres.2024.127790.

## References

Abraham, N.M., Liu, L., Jutras, B.L., Yadav, A.K., Narasimhan, S., Gopalakrishnan, V., Ansari, J.M., Jefferson, K.K., Cava, F., Jacobs-Wagner, C., Fikrig, E., 2017. Pathogen-mediated manipulation of arthropod microbiota to promote infection. Proc. Natl. Acad. Sci. **114** (5) https://doi.org/10.1073/pnas.1613422114.

Abuin-Denis, L., Piloto-Sardinas, E., Maître, A., Wu-Chuang, A., Mateos-Hernández, L., Obregon, D., Corona-González, B., Fogaça, A.C., Palinauskas, V., Azélytè, J., Rodríguez-Mallon, A., Cabezas-Cruz, A., 2024. Exploring the impact of Anaplasma phagocytophilum on colonization resistance of Ixodes scapularis microbiota using network node manipulation. Curr. Res. Parasitol. Vector-Borne Dis. **5**, 100177 https://doi.org/10.1016/j.crpvbd.2024.100177.

Adegoke, A., Kumar, D., Budachetik, K., Karim, S., 2022. Hematophagy and tick-borne Rickettsial pathogen shape the microbial community structure and predicted functions within the tick vector, Amblyomma maculatum. Front. Cell. Infect. Microbiol. **12**, 1037387 https://doi.org/10.3389/fcimb.2022.1037387.

Agler, M.T., Ruhe, J., Kroll, S., Morhenn, C., Kim, S.-T., Weigel, D., Kemen, E.M., 2016. Microbial Hub Taxa Link Host and Abiotic Factors to Plant Microbiome Variation. PLOS Biol. **14** (1), e1002352 https://doi.org/10.1371/journal.pbio.1002352.

Aiveloo, T., Norberg, A., Tsir�in, 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.

Almazán, C., Scimeca, R.C., Reichard, M.V., Mosqueda, J., 2022. Babesiosis and Theileriosis in North America. Pathogens **11** (2), 168. https://doi.org/10.3390/pathogens11020168.

Aubry, P., Geale, D.W., 2011. A Review of Bovine Anaplasmosis: Review of Bovine Anaplasmosis. Transound. Emerg. Dis. **58** (1), 1-30. https://doi.org/10.1111/j.1865-1682.2010.01173.x.

Azélytè, J., Wu-Chuang, A., Ziegytè, R., Platonova, E., Mateos-Hernandez, L., Maye, J., Obregon, D., Palinauskas, V., Cabezas-Cruz, A., 2022. Anti-Microbiota Vaccine Reduces Avian Malaria Infection Within Mosquito Vectors. Front. Immunol. **13**, 841835 https://doi.org/10.3389/fimmu.2022.841835.

Bastian, M., Heymann, S., & Jacomy, M. (2009). Gephi: An Open Source Software for Exploring and Manipulating Networks. WebAtlas.

Binetruy, F., Buyse, M., Lejarre, Q., Barosi, R., Villa, M., Rahola, N., Paupy, C., Ayala, D., Duron, O., 2020. Microbial community structure reveals instability of nutritional symbiosis during the evolutionary radiation of *Amblyomma* ticks. Mol. Ecol. **29** (5), 1016-1029. https://doi.org/10.1111/mec.15373.

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

Bray, J.R., Curtis, J.T., 1957. An Ordination of the Upland Forest Communities of Southern Wisconsin. Ecol. Monogr. **27** (4), 325-349. https://doi.org/10.2307/1942268.

Brownlie, J.C., Johnson, K.N., 2009. Symbiont-mediated protection in insect hosts. Trends Microbiol. **17** (8), 348-354. https://doi.org/10.1016/j.tim.2009.05.005.

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

Cabezas-Cruz, A., Passos, L.M.F., Lis, K., Kennel, R., Verdés, J.J., Ferrolohn, J., Tonk, M., Pohl, A.E., Grubhoffer, L., Zweygarth, E., Shkap, V., Ribeiro, M.F.B., Estrada-Pena, A., Kocan, K.M., De La Fuente, J., 2013. Functional and Immunological Relevance of Anaplasma marginale Major Surface Protein 1a Sequence and Structural Analysis. PLoS ONE **8** (6), e65243. https://doi.org/10.1371/journal.pone.0065243.

Cabezas-Cruz, A., Zweygarth, E., Vancová, M., Broniszewska, M., Grubhoffer, L., Passos, L.M.F., Ribeiro, M.F.B., Alberdi, P., De La Fuente, J., 2016. *Ehrlichia minasensis* sp. Nov., isolated from the tick *Rhipicephalus microplus*. Int. J. Syst. Evolut. Microbiol. **66** (3), 1426-1430. https://doi.org/10.1099/ijsem.0.000895.

Cabezas-Cruz, A., Zweygarth, E., Aguiar, D.M., 2019. *Ehrlichia minasensis*, an old demon with a new name. Ticks Tick. Borne Dis. **10** (4), 828-829. https://doi.org/10.1016/j.ttbdis.2019.03.018.

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** (7), 581-583. https://doi.org/10.1038/nmeth.3869.

Carvalho, I.T.S., Melo, A.L.T., Freitas, L.C., Verçoza, R.V., Alves, A.S., Costa, J.S., Chitarra, C.S., Nakazato, L., Dutra, V., Pacheco, R.C., Aguiar, D.M., 2016. Minimum infection rate of *Ehrlichia minasensis* in *Rhipicephalus microplus* and *Amblyomma sculptum* ticks in Brazil. Ticks Tick. Borne Dis. **7** (5), 849-852. https://doi.org/10.1016/j.ttbdis.2016.04.004.

Cobo-López, S., Gupta, V.K., Sung, J., Guimerà, R., Sales-Pardo, M., 2022. Stochastic block models reveal a robust nested pattern in healthy human gut microbiomes. PNAS Nexus **1** (3), pgac055. https://doi.org/10.1093/pnasexus/pgac055.

De La Fournière, S., Guillemi, E.C., Paoletta, M.S., Pérez, A., Obregón, D., Cabezas-Cruz, A., Sarmiento, N.F., Farber, M.D., 2023. Transovarial Transmission of *Anaplasma marginale* in *Rhipicephalus (Boophilus) microplus* Ticks Results in a

Bottleneck for Strain Diversity. Pathogens **12** (8), 1010. https://doi.org/10.3390/pathogens12081010.

DeSantis, T.Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E.L., Keller, K., Huber, T., Dalevi, D., Hu, P., Andersen, G.L., 2006. Greengenes, A. Chimera-Checked 16S rRNA Gene Database and Workbench Compatible with ARB. Appl. Environ. Microbiol. **72** (7), 5069-5072. https://doi.org/10.1128/AEM.03006-05.

Duron, O., Gottlieb, Y., 2020. Convergence of Nutritional Symbionts in Obligate Blood Feeders. Trends Parasitol. **36** (10), 816-825. https://doi.org/10.1016/j.pt.2020.07.007.

Duron, O., Jourdain, E., McCoy, K.D., 2014. Diversity and global distribution of the *Coxiella* intracellular bacterium in seabird ticks. Ticks Tick. Borne Dis. **5** (5), 557-563. https://doi.org/10.1016/j.ttbdis.2014.04.003.

Encinosa Guzmán, P.E., Fernández Cuétara, C., Cano Argüelles, A.L., Fuentes Castillo, A., García Martínez, Y., Rodríguez Fernández, R., Fernández Afonso, Y., Bello Soto, Y., González Alfaro, Y., Méndez, L., Díaz García, A., Estrada, M.P., Rodríguez-Mallon, A., 2021. Characterization of two Cuban colonies of *Rhipicephalus microplus* ticks. Vet. Parasitol.: Reg. Stud. Rep. **25**, 100591 https://doi.org/10.1016/j.vpsr.2021.100591.

Estrada-Pena, A., Naranjo, V., Acevedo-Whitehouse, K., Mangold, A.J., Kocan, K.M., De La Fuente, J., 2009. Phylogeographic analysis reveals association of tick-borne pathogen, Anaplasma marginale, MSP1a sequences with ecological traits affecting tick vector performance. BMC Biol. **7** (1), 57. https://doi.org/10.1186/1741-7007-7-57.

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

Fernandes, A.D., Macklain, J.M., Linn, T.G., Reid, G., Gloor, G.B., 2013. ANOVA-Like Differential Expression (ALDEX) Analysis for Mixed Population RNA-Seq PLoS ONE **8** (7), e67019. https://doi.org/10.1371/journal.pone.0067019.

Fogaça, A.C., Sousa, G., Pavanelo, D.B., Esteves, E., Martins, L.A., Urbanová, V., Kopáček, P., Daffre, S., 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.

Friedman, J., Alm, E.J., 2012. Inferring Correlation Networks from Genomic Survey Data. PLoS Comput. Biol. **8** (9), e1002687 https://doi.org/10.1371/journal.pcbi.1002687.

de la Fuente, J., Antunes, S., Bonnet, S., Cabezas-Cruz, A., Domingos, A.G., Estrada-Pena, A., Johnson, N., Kocan, K.M., Mansfield, K.L., Nijhof, A.M., Papa, A., Rudenko, N., Villar, M., Alberdi, P., Torina, A., Ayllón, N., Vancova, M., Golovchenko, M., Grubhoffer, L., Rego, R.O.M., 2017. Tick-Pathogen Interactions and Vector Competence: Identification of Molecular Drivers for Tick-Borne Diseases. Front. Cell. Infect. Microbiol. **7** https://doi.org/10.3389/fcimb.2017.00114.

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** (8), 1846-1855. https://doi.org/10.1038/ismej.2015.266.

Gomard, Y., Flores, O., Vittecoq, M., Blanchon, T., Toty, C., Duron, O., Mavingui, P., Tortosa, P., McCoy, K.D., 2021. Changes in Bacterial Diversity, Composition and Interactions During the Development of the Seabird Tick *Ornithodorus maritimus* (Argasidae). Microb. Ecol. **81** (3), 770-783. https://doi.org/10.1007/s00248-020-01611-9.

Guimarães, A.M., Lima, J.D., Ribeiro, M.F.B., 1998. Sporogony and experimental transmission of Babesia equi by *Boophilus microplus*. Parasitol. Res. **84** (4), 323-327. https://doi.org/10.1007/s004360050404.

Guizzo, M.G., Tirloni, L., Gonzalez, S.A., Farber, M.D., Braz, G., Parizi, L.F., Davidav E Silva, L.A., Da Silva Vaz, I., Oliveira, P.L., 2022. *Coxiella* Endosymbiont of *Rhipicephalus microplus* Modulates Tick Physiology With a Major Impact in Blood Feeding Capacity. Front. Microbiol. **13**, 868575 https://doi.org/10.3389/fmicb.2022.868575.

Hamilton, P.T., Maluenda, S., Sarr, A., Belli, A., Hurry, G., Duron, O., Plantard, O., Voordouw, M.J., 2021. *Borrelia afzelii* Infection in the Rodent Host Has Dramatic Effects on the Bacterial Microbiome of *Ixodes ricinus* Ticks. Appl. Environ. Microbiol. **87** (18), e00641-21 https://doi.org/10.1128/AEM.00641-21.

Hammer, Ø., Harper, D.A.T., Ryan, P.D., 2001. PAST: paleontologicalstatistics software package for education and data analysis. Palaeontol. Electron.

Kalil, S.P., Rosa, R.D.D., Capelli-Peixoto, J., Pohl, P.C., Oliveira, P.L.D., Fogaça, A.C., Daffre, S., 2017. Immune-related redox metabolism of embryonic cells of the tick *Rhipicephalus microplus* (BME26) in response to infection with Anaplasma marginale. Parasites Vectors **10** (613), 613. https://doi.org/10.1186/s13071-017-2575-9.

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

Kenneil, R., Shkap, V., Leibovich, B., Zweygarth, E., Pfister, K., Ribeiro, M.F.B., Passos, L. M.F., 2013. Cross-Protection Between Geographically Distinct *Anaplasma marginale* Isolates Appears to be Constrained by Limited Antibody Responses. Transbound. Emerg. Dis. **60**, 97-104. https://doi.org/10.1111/tbed.12125.

Kocan, K.M., De La Fuente, J., Blouin, E.F., Garcia-Garcia, J.C., 2004. *Anaplasma marginale* (Rickettsiales: Anaplasmataceae): recent advances in defining host-pathogen adaptations of a tick-borne rickettsiae. Parasitology **129** (S1), S285-S300. https://doi.org/10.1016/S0031182003004700.

Kocan, K.M., De La Fuente, F.J., Cabezas-Cruz, A., 2015. The genus *Anaplasma*: New challenges after reclassification: -EN- -FR- Les bactéries du genre *Anaplasma*: nouveaux défis suite à la reclassification -ES- Nuevas dificultades tras la reclassificación del género *Anaplasma*. Rev. Sci. Et. Tech. De. l'OIE **34** (2), 577-586. https://doi.org/10.20506/rst.34.2.2381.

Krawczyk, A.I., Röttjers, S., Coimbra-Dores, M.J., Heylen, D., Fonville, M., Takken, W., Faust, K., Sprong, H., 2022. Tick microbial associations at the crossroad of horizontal
and vertical transmission pathways. Parasites Vectors **15** (1), 380. https://doi.org/10.1186/s13071-022-05519-w.

Lejal, E., Chiquet, J., Aubert, J., Robin, S., Estrada-Pena, A., Rue, O., Midoux, C., Mariadassou, M., Bailly, X., Cougoul, A., Gasqui, P., Cosson, J.F., Chalvet-Monfray, K., Vayssier-Taussat, M., Pollet, T., 2021. Temporal patterns in *Ixodes ricinus* microbial communities: An insight into tick-borne microbe interactions. *Microbiome* 9 (1), 153. 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. https://doi.org/10.4000/cybergeo.26763.

Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Foucault-Simonin, A., Moutailler, S., Paoli, J.-C., Falchi, A., Díaz-Sánchez, A.A., Banović, P., Obregón, D., Cabezas-Cruz, A., 2022. *Rickettsia helvetica* infection is associated with microbiome modulation in *Ixodes ricinus* collected from humans in Serbia. Sci. Rep. **12** (1), 11464 https://doi.org/10.1038/s41598-022-15681-x.

Maitre, A., Wu-Chuang, A., Azelyté, J., Palinauskas, V., Mateos-Hernández, L., Obregon, D., Hodžić, A., Valiente Moro, C., Estrada-Pena, A., Paoli, J.-C., Falchi, A., Cabezas-Cruz, A., 2022. Vector microbiota manipulation by host antibodies: The forgotten strategy to develop transmission-blocking vaccines. Parasites Vectors **15** (1), 4. https://doi.org/10.1186/s13071-021-05122-5.

Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Piloto-Sardiñas, E., Foucault-Simonin, A., Cicculli, V., Moutailler, S., Paoli, J., Falchi, A., Obregón, D., Cabezas-Cruz, A., 2023. Rickettsial pathogens drive microbiota assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks. Mol. Ecol. **32** (16), 4660-4676. https://doi.org/10.1111/mec.17058.

Mariani, M.S., Ren, Z.-M., Bascompte, J., Tessone, C.J., 2019. Nestedness in complex networks: Observation, emergence, and implications. Phys. Rep. **813**, 1-90. https://doi.org/10.1016/j.physrep.2019.04.001.

Mateos-Hernández, L., Obregón, D., Maye, J., Borneres, J., Versille, N., de la Fuente, J., Estrada-Pena, A., Hodžić, A., Simo, L., Cabezas-Cruz, A., 2020. Anti-Tick Microbiota Vaccine Impacts *Ixodes ricinus* Performance during Feeding. Vaccines **8** (4), 702. https://doi.org/10.3390/vaccines8040702.

Mateos-Hernández, L., Obregón, D., Wu-Chuang, A., Maye, J., Bornères, J., Versillé, N., de la Fuente, J., Díaz-Sánchez, S., Bermúdez-Humarán, L.G., Torres-Maravilla, E., Estrada-Pena, A., Hodžić, A., Simo, L., Cabezas-Cruz, A., 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.

McCutcheon, J.P., von Dohlen, C.D., 2011. An Interdependent Metabolic Patchwork in the Nested Symbiosis of Mealybugs. Curr. Biol. **21** (16), 1366-1372. https://doi.org/10.1016/j.cub.2011.06.051.

Mercado-Curiel, R.F., Ávila-Ramírez, M.L., Palmer, G.H., Brayton, K.A., 2014. Identification of *Rhipicephalus microplus* Genes That Modulate the Infection Rate of the Rickettsia Anaplasma marginale. PLoS ONE **9** (3), e91062. https://doi.org/10.1371/journal.pone.0091062.

Mesquita, E., Da Costa, D.P., Meirelles, L.N., Camargo, M.G., Corrêa, T.A., Bittencourt, V. R.E.P., Da Silva Coelho, I., Santos, H.A., Humber, R.A., Golo, P.S., 2023. Entomopathogenic fungus treatment changes the gut bacterial diversity of *Rhipicephalus microplus* ticks. Parasites Vectors **16** (1), 185. https://doi.org/10.1186/s13071-023-05790-5.

Moura De Aguiar, D., Pessoa Araújo Junior, J., Nakazato, L., Bard, E., Aguilar-Bultet, L., Vorimore, F., Leonidovich Popov, V., Moleta Colodel, E., Cabezas-Cruz, A., 2019. Isolation and Characterization of a Novel Pathogenic Strain of *Ehrlichia minasensis*. Microorganisms **7** (11), 528. https://doi.org/10.3390/microorganisms7110528.

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. Archn. Sci. **2**, 1297733 https://doi.org/10.3389/frchs.2023.1297733.

Pielou, E.C., 1966. Species-diversity and pattern-diversity in the study of ecological succession. J. Theor. Biol. **10** (2), 370-383. https://doi.org/10.1016/0022-5193(66)90133-0.

Piloto-Sardiñas, E., Foucault-Simonin, A., Wu-Chuang, A., Mateos-Hernández, L., Marrero-Perera, R., Abuin-Denis, L., Roblejo-Arias, L., Díaz-Corona, C., Zająć, Z., Kulisz, J., Woźniak, A., Moutailler, S., Corona-González, B., Cabezas-Cruz, A., 2023.

Dynamics of Infections in Cattle and *Rhipicephalus microplus*: A Preliminary Study. Pathogens **12** (8), 998. https://doi.org/10.3390/pathogens12080998.

Piloto-Sardiñas, E., Abuin-Denis, L., Maitre, A., Foucault-Simonin, A., Corona-González, B., Díaz-Corona, C., Roblejo-Arias, L., Mateos-Hernández, L., Marrero-Perera, R., Obregon, D., Svobodová, K., Wu-Chuang, A., Cabezas-Cruz, A., 2024. Dynamic nesting of *Anaplasma marginale* in the microbial communities of *Rhipicephalus microplus*. Ecol. Evol. **14** (4), e11228 https://doi.org/10.1002/eece3.11228.

Price, M.N., Dehal, P.S., Arkin, A.P., 2010. FastTree 2 – Approximately Maximum-Likelihood Trees for Large Alignments. PLoS ONE **5** (3), e9490. https://doi.org/10.1371/journal.pone.0009490.

R Core Team (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing [Computer software]. https://www.R-project.org/.

Ramírez-Hernández, A., Arroyave, E., Faccini-Martínez, Á.A., Martínez-Diaz, H.C., Betancourt-Ruiz, P., Olaya-M, L.-A., Forero-Becerra, E.G., Hidalgo, M., Blanton, L.S., Walker, D.H., 2022. Emerging Tickborne Bacteria in Cattle from Colombia. Emerg. Infect. Dis. **28** (10), 2109-2111. https://doi.org/10.3201/eid2810.220657.

Real, R., Vargas, J.M., 1996. The Probabilistic Basis of Jaccard's Index of Similarity. Syst. Biol. **45** (3), 380-385. https://doi.org/10.1093/sysbio/45.3.380.

Rialch, A., Sankar, M., Silamparasan, M., Madhusoodan, A.P., Kharayat, N.S., Gautam, S., Gurav, A.R., Thankappan, S., 2022. Molecular detection of *Coxiella*-like endosymbionts in *Rhipicephalus microplus* from north India. Vet. Parasitol.: Reg. Stud. Rep. **36**, 100803 https://doi.org/10.1016/j.vprsr.2022.100803.

Rojas-Jaimes, J., Lindo-Seminario, D., Correa-Núnez, G., Diringer, B., 2021. Characterization of the bacterial microbiome of *Rhipicephalus (Boophilus) microplus* collected from Pecari tajacu "Sajino" Madre de Dios, Peru. Sci. Rep. **11** (1), 6661. https://doi.org/10.1038/s41598-021-86177-3.

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

RStudio Team (2020). RStudio: Integrated development for R. (http://www.rstudio.org/).

Ruiling, Z., Zhendong, H., Guangfu, Y., Zhong, Z., 2019. Characterization of the bacterial community in *Haemaphysalis longicornis* (Acari: Ixodidae) throughout developmental stages. Exp. Appl. Acarol. **77** (2), 173-186. https://doi.org/10.1007/s10493-019-00393-7.

Segura, J.A., Isaza, J.P., Botero, L.E., Alzate, J.F., Gutiérrez, L.A., 2020. Assessment of bacterial diversity of *Rhipicephalus microplus* ticks from two livestock agroecosystems in Antioquia, Colombia. PLOS ONE **15** (7), e0234005. https://doi.org/10.1371/journal.pone.0234005.

Smith, C.A., Ashby, B., 2023. Tolerance-conferring defensive symbionts and the evolution of parasite virulence. Evol. Lett. **7** (4), 262-272. https://doi.org/10.1093/evlett/grad015.

Song, C., Rohr, R.P., Saavedra, S., 2017. Why are some plant-pollinator networks more nested than others? J. Anim. Ecol. **86** (6), 1417-1424. https://doi.org/10.1111/1365-2656.12749.

Wu-Chuang, A., Mateos-Hernandez, L., Maitre, A., Rego, R.O.M., Síma, R., Porcelli, S., Rakotobe, S., Foucault-Simonin, A., Moutailler, S., Palinauskas, V., Azelyté, J., Simo, L., Obregon, D., Cabezas-Cruz, A., 2023. Microbiota perturbation by anti-microbiota vaccine reduces the colonization of *Borrelia afzelii* in *Ixodes ricinus*. Microbiome **11** (1), 151. https://doi.org/10.1186/s40168-023-01599-7.

Ye, J., Coulouris, G., Zaretskaya, I., Cutcutache, I., Rozen, S., Madden, T.L., 2012. Primer-BLAST: A tool to design target-specific primers for polymerase chain reaction. BMC Bioinforma. **13** (1), 134. https://doi.org/10.1186/1471-2105-13-134.

Zhang, R., Yu, G., Huang, Z., Zhang, Z., 2020. Microbiota assessment across different developmental stages of *Dermacentor silvarum* (Acari: Ixodidae) revealed stage-specific signatures. Ticks Tick. Borne Dis. **11** (2), 101321 https://doi.org/10.1016/j.ttbdis.2019.101321.

Zhong, J., 2012. Coxiella-like Endosymbionts. In: Toman, R., Heinzen, R.A., Samuel, J.E., Mege, J.-L. (Eds.), Coxiella burnetii: Recent Advances and New Perspectives in Research of the Q Fever Bacterium, Vol. 984. Springer Netherlands, pp. 365-379. https://doi.org/10.1007/978-94-007-4315-1_18.