# Dynamic nesting of Anaplasma marginale in the microbial communities of Rhipicephalus microplus

Elianne Piloto-Sardiñas, Lianet Abuin-Denis, Apolline Maitre, Angélique Foucault-Simonin, Belkis Corona-González, Cristian Díaz-Corona, Lisset Roblejo-Arias, Lourdes Mateos-Hernández, Roxana Marrero-Perera, Dasiel Obregon, et al.

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

## To cite this version:

Elianne Piloto-Sardiñas, Lianet Abuin-Denis, Apolline Maitre, Angélique Foucault-Simonin, Belkis Corona-González, et al. Dynamic nesting of Anaplasma marginale in the microbial communities of Rhipicephalus microplus. Ecology and Evolution, 2024, 14 (4), 10.1002/ece3.11228. hal-04557489

HAL Id: hal-04557489

https://hal.science/hal-04557489

Submitted on 24 Apr 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>

## RESEARCH ARTICLE

Ecology and Evolution
Open Access

WILEY

# Dynamic nesting of *Anaplasma marginale* in the microbial communities of *Rhipicephalus microplus*

Elianne Piloto-Sardiñas<fnref n="1" />,<fnref n="2" />

Lianet Abuin-Denis<fnref n="1" />,<fnref n="3" />

Apolline Maitre<fnref n="1" />,<fnref n="4" />,<fnref n="5" />

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

Belkis Corona-González<fnref n="2" />

Cristian Díaz-Corona<fnref n="2" />

Lisset Roblejo-Arias<fnref n="2" />

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

Roxana Marrero-Perera<fnref n="2" />

Dasiel Obregon<fnref n="6" /> ID

Karolína Svobodová<fnref n="7" />

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

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

### Correspondence

Alejandro Cabezas-Cruz and Elianne Piloto-Sardiñas, ANSES, INRAE, Ecole Nationale Vétérinaire d'Alfort, UMR BIPAR, Laboratoire de Santé Animale, Maisons-Alfort, France.

Email: alejandro.cabezas@vet-alfort.fr and elianne9409@gmail.com

### Funding information

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

### Abstract

Interactions within the tick microbiome involving symbionts, commensals, and tick-borne pathogens (TBPs) play a pivotal role in disease ecology. This study explored temporal changes in the microbiome of *Rhipicephalus microplus*, an important cattle tick vector, focusing on its interaction with *Anaplasma marginale*. To overcome limitations inherent in sampling methods relying on questing ticks, which may not consistently reflect pathogen presence due to variations in exposure to infected hosts in nature, our study focused on ticks fed on chronically infected cattle. This approach ensures continuous pathogen exposure, providing a more comprehensive understanding of the nesting patterns of *A. marginale* in the *R. microplus* microbiome. Using next-generation sequencing, microbiome dynamics were characterized over 2 years, revealing significant shifts in diversity, composition, and abundance. *Anaplasma marginale* exhibited varying associations, with its increased abundance correlating with reduced microbial diversity. Co-occurrence networks demonstrated *Anaplasma's* evolving role, transitioning from diverse connections to keystone taxa status. An integrative approach involving *in silico* node removal unveils the impact of *Anaplasma* on network stability, highlighting its role in conferring robustness to the microbial community. This study provides insights into the intricate interplay between the tick microbiome and *A. marginale*, shedding light on potential avenues for controlling bovine anaplasmosis through microbiome manipulation.

### KEYWORDS

*Anaplasma marginale*, microbiome dynamics, nesting, networks, *Rhipicephalus microplus*, ticks

### TAXONOMY CLASSIFICATION

Community ecology

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

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

<footnote n="3">³Animal Biotechnology Department, Center for Genetic Engineering and Biotechnology, Havana, Cuba</footnote>

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

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

<footnote n="6">⁶School of Environmental Sciences, University of Guelph, Guelph, Ontario, Canada</footnote>

<footnote n="7">⁷Faculty of Science, University of South Bohemia, Ceske Budejovice, Czech Republic</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
© 2024 The Authors. *Ecology and Evolution* published by John Wiley & Sons Ltd.</footnote>
<figure><img src="image_2.png" /></figure>

## RESEARCH ARTICLE

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

WILEY

# Dynamic nesting of *Anaplasma marginale* in the microbial communities of *Rhipicephalus microplus*

Elianne Piloto-Sardiñas<fnref n="1" />,<fnref n="2" /> | Lianet Abuin-Denis<fnref n="1" />,<fnref n="3" /> | Apolline Maitre<fnref n="1" />,<fnref n="4" />,<fnref n="5" /> |

Angélique Foucault-Simonin<fnref n="1" /> | Belkis Corona-González<fnref n="2" /> | Cristian Díaz-Corona<fnref n="2" /> |

Lisset Roblejo-Arias<fnref n="2" /> | Lourdes Mateos-Hernández<fnref n="1" /> | Roxana Marrero-Perera<fnref n="2" /> |

Dasiel Obregon<fnref n="6" /> ID | Karolína Svobodová<fnref n="7" /> | Alejandra Wu-Chuang<fnref n="1" /> |

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

### Correspondence

Alejandro Cabezas-Cruz and Elianne Piloto-Sardiñas, ANSES, INRAE, Ecole Nationale Vétérinaire d'Alfort, UMR BIPAR, Laboratoire de Santé Animale, Maisons-Alfort, France.

Email: alejandro.cabezas@vet-alfort.fr and elianne9409@gmail.com

### Funding information

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

### Abstract

Interactions within the tick microbiome involving symbionts, commensals, and tick-borne pathogens (TBPs) play a pivotal role in disease ecology. This study explored temporal changes in the microbiome of *Rhipicephalus microplus*, an important cattle tick vector, focusing on its interaction with *Anaplasma marginale*. To overcome limitations inherent in sampling methods relying on questing ticks, which may not consistently reflect pathogen presence due to variations in exposure to infected hosts in nature, our study focused on ticks fed on chronically infected cattle. This approach ensures continuous pathogen exposure, providing a more comprehensive understanding of the nesting patterns of *A. marginale* in the *R. microplus* microbiome. Using next-generation sequencing, microbiome dynamics were characterized over 2 years, revealing significant shifts in diversity, composition, and abundance. *Anaplasma marginale* exhibited varying associations, with its increased abundance correlating with reduced microbial diversity. Co-occurrence networks demonstrated *Anaplasma's* evolving role, transitioning from diverse connections to keystone taxa status. An integrative approach involving *in silico* node removal unveils the impact of *Anaplasma* on network stability, highlighting its role in conferring robustness to the microbial community. This study provides insights into the intricate interplay between the tick microbiome and *A. marginale*, shedding light on potential avenues for controlling bovine anaplasmosis through microbiome manipulation.

### KEYWORDS

*Anaplasma marginale*, microbiome dynamics, nesting, networks, *Rhipicephalus microplus*, ticks

### TAXONOMY CLASSIFICATION

Community ecology

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

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

<footnote n="3">³Animal Biotechnology Department, Center for Genetic Engineering and Biotechnology, Havana, Cuba</footnote>

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

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

<footnote n="6">⁶School of Environmental Sciences, University of Guelph, Guelph, Ontario, Canada</footnote>

<footnote n="7">⁷Faculty of Science, University of South Bohemia, Ceske Budejovice, Czech Republic</footnote>

<footnote>This is an open access article under the terms of the [Creative Commons Attribution](https://creativecommons.org/licenses/by) License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
© 2024 The Authors. *Ecology and Evolution* published by John Wiley & Sons Ltd.</footnote>
# 1 | INTRODUCTION

Interactions between symbionts, commensals, and tick-borne pathogens (TBPs) within the tick microbiome can potentially shape disease ecology. Symbiont-microbiome interactions in questing *Ixodes ricinus* can be dynamic, with variations observed in the prevalence and distribution of tick symbionts across different forest sites (Krawczyk et al., 2022). The strongest determinants of microbiome clustering were found to be the abundance and prevalence of specific symbionts, such as *Rickettsia* and *Rickettsiella*. The proportions of these symbionts varied between geographically close forest sites, suggesting a potential spatial scale influencing their distribution. The observed variations in the prevalence of tick symbionts were not consistent with horizontally transmitted pathogens such as *Borrelia afzelii*, *Borrelia garinii*, *Anaplasma phagocytophilum*, and *Neorickettsia mikurensis*, which showed more random patterns across geographically close forest sites (Krawczyk et al., 2022). A finding supported by previous studies in different tick species including *I. ricinus* (Lejal et al., 2019; Zając et al., 2023), *Dermacentor reticulatus* (Zając et al., 2023), and *Rhipicephalus microplus* (Piloto-Sardiñas, Foucault-Simonin, et al., 2023). This distinction suggests that the factors influencing symbiont prevalence may differ from those affecting horizontally transmitted pathogens, which are mainly determined by local vertebrate communities (Takumi et al., 2019).

Microbe-microbe associations within the tick microbiome, particularly those involving TBPs and other nonpathogenic bacteria, exhibit dynamic patterns over time (Lejal et al., 2021). The temporal dynamics of the *I. ricinus* microbiome and its impact on microbiome-pathogen interactions were evaluated in questing ticks collected during three consecutive years in a peri-urban forest in France (Lejal et al., 2021). Results revealed temporal variations in the microbiome, with distinct clusters of tick samples collected during different months. *Wolbachia*, *Arsenophonus*, *Spiroplasma*, and *Pseudomonas* were identified as drivers of certain clusters, indicating their role in shaping temporal variations (Lejal et al., 2021). Notably, comparisons of tick samples positive for specific TBPs (*Rickettsia*, *Borrelia*, and *Anaplasma*) with TBP-negative samples demonstrated significantly higher abundance of relevant operational taxonomic units (OTUs) in TBP-positive samples. Network analyses revealed varying correlation patterns between TBP-positive and TBP-negative samples, with *Borrelia*-positive samples showing both positive and negative correlations with specific OTUs, including those associated with environmental and pathogenic/symbiotic genera (Lejal et al., 2021).

Temporal changes in the vector microbiome may affect the ability of pathogens to persist in vectors, transmit to hosts, and cause disease, as the transmission of vector-borne pathogens often involves complex interactions between pathogens and the microbiome (Abraham et al., 2017; Maitre et al., 2022, 2023; Narasimhan et al., 2014, 2017). While the studies by Krawczyk et al. (2022), and Lejal et al. (2021) provided valuable insights into the dynamic nature of symbiont-microbiome interactions and the temporal dynamics of TBPs within questing ticks, it is essential to acknowledge a crucial limitation inherent to this sampling approach. The absence of a

pathogen in a questing tick may be linked to the fact that the host on which the tick fed did not harbor the pathogen in the first place (Takumi et al., 2019). This inherent limitation poses a challenge in disentangling whether the absence of a TBP in a tick is due to the absence of the pathogen in the previous host or influenced by the tick microbiome.

This sampling strategy also prevents assessing nesting dynamics of TBPs in the tick microbiome, as dosage and frequency of pathogen exposure are difficult to control in natural settings. To address the challenge associated with the limitations of questing tick sampling and better elucidate TBP-microbiome interactions in nature, we propose exploring systems involving ticks feeding on chronically infected hosts. In such scenarios, ticks would encounter a continuous presence of the pathogen across time. This approach, exemplified by systems such as *R. microplus* feeding on cattle chronically infected with *Anaplasma marginale*, provides an unique opportunity to assess the impact of the microbiome on TBP dynamics under more controlled conditions.

*Rhipicephalus microplus*, commonly known as the cattle tick or tropical cattle tick, is a significant ectoparasite that infests cattle and other livestock. This tick species holds great importance due to its capacity to transmit various pathogens, with *A. marginale* being one of the most notable (De La Fourniere et al., 2023; Pereira et al., 2022). *Anaplasma marginale* is a bacterium that causes bovine anaplasmosis, a disease characterized by anemia, fever, and other clinical symptoms in cattle (Salinas-Estrella et al., 2022). This disease can lead to significant economic losses due to decreased productivity, increased veterinary costs, and even livestock mortality (Rodríguez et al., 2009; Salinas-Estrella et al., 2022). *Rhipicephalus microplus* plays a crucial role in the transmission of *A. marginale*, as it acts as a vector by feeding on infected cattle and subsequently transmitting the bacteria to susceptible animals during subsequent feedings (Zivkovic et al., 2010). In a recent study, it was observed that while *A. marginale* consistently infected all cattle across different sampling periods, its presence was not uniformly detected in all *R. microplus* infesting the cattle (Piloto-Sardiñas, Foucault-Simonin, et al., 2023). This indicates nonlinearity between tick infestation rate and pathogen prevalence in ticks (Ostfeld & Keesing, 2023), likely under strong influence of the tick microbiome (Tonk-Rügen et al., 2023).

The present study aimed to investigate whether temporal fluctuations in the microbial communities of the cattle tick *R. microplus* could potentially disturb or alter the interactions between *A. marginale* and the tick microbial communities and as a consequence shape the impact of the pathogen on the microbial community assembly. To achieve this, *R. microplus* samples confirmed for infection with *A. marginale* were used (Piloto-Sardiñas, Foucault-Simonin, et al., 2023), and next-generation sequencing (NGS) was performed for microbiome characterization. Employing an experimental network approach, an *in silico* node removal technique was applied, a strategy previously employed to investigate the impact of *Rickettsia* pathogens on the microbiome assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks (Maitre et al., 2023). By simulating the absence of specific TBPs such as *Anaplasma* sp. in *silico*, the objective
# 1 | INTRODUCTION

Interactions between symbionts, commensals, and tick-borne pathogens (TBPs) within the tick microbiome can potentially shape disease ecology. Symbiont-microbiome interactions in questing *Ixodes ricinus* can be dynamic, with variations observed in the prevalence and distribution of tick symbionts across different forest sites (Krawczyk et al., 2022). The strongest determinants of microbiome clustering were found to be the abundance and prevalence of specific symbionts, such as *Rickettsia* and *Rickettsiella*. The proportions of these symbionts varied between geographically close forest sites, suggesting a potential spatial scale influencing their distribution. The observed variations in the prevalence of tick symbionts were not consistent with horizontally transmitted pathogens such as *Borrelia afzelii*, *Borrelia garinii*, *Anaplasma phagocytophilum*, and *Neorickettsia mikurensis*, which showed more random patterns across geographically close forest sites (Krawczyk et al., 2022). A finding supported by previous studies in different tick species including *I. ricinus* (Lejal et al., 2019; Zając et al., 2023), *Dermacentor reticulatus* (Zając et al., 2023), and *Rhipicephalus microplus* (Piloto-Sardiñas, Foucault-Simonin, et al., 2023). This distinction suggests that the factors influencing symbiont prevalence may differ from those affecting horizontally transmitted pathogens, which are mainly determined by local vertebrate communities (Takumi et al., 2019).

Microbe-microbe associations within the tick microbiome, particularly those involving TBPs and other nonpathogenic bacteria, exhibit dynamic patterns over time (Lejal et al., 2021). The temporal dynamics of the *I. ricinus* microbiome and its impact on microbiome-pathogen interactions were evaluated in questing ticks collected during three consecutive years in a peri-urban forest in France (Lejal et al., 2021). Results revealed temporal variations in the microbiome, with distinct clusters of tick samples collected during different months. *Wolbachia*, *Arsenophonus*, *Spiroplasma*, and *Pseudomonas* were identified as drivers of certain clusters, indicating their role in shaping temporal variations (Lejal et al., 2021). Notably, comparisons of tick samples positive for specific TBPs (*Rickettsia*, *Borrelia*, and *Anaplasma*) with TBP-negative samples demonstrated significantly higher abundance of relevant operational taxonomic units (OTUs) in TBP-positive samples. Network analyses revealed varying correlation patterns between TBP-positive and TBP-negative samples, with *Borrelia*-positive samples showing both positive and negative correlations with specific OTUs, including those associated with environmental and pathogenic/symbiotic genera (Lejal et al., 2021).

Temporal changes in the vector microbiome may affect the ability of pathogens to persist in vectors, transmit to hosts, and cause disease, as the transmission of vector-borne pathogens often involves complex interactions between pathogens and the microbiome (Abraham et al., 2017; Maitre et al., 2022, 2023; Narasimhan et al., 2014, 2017). While the studies by Krawczyk et al. (2022), and Lejal et al. (2021) provided valuable insights into the dynamic nature of symbiont-microbiome interactions and the temporal dynamics of TBPs within questing ticks, it is essential to acknowledge a crucial limitation inherent to this sampling approach. The absence of a

pathogen in a questing tick may be linked to the fact that the host on which the tick fed did not harbor the pathogen in the first place (Takumi et al., 2019). This inherent limitation poses a challenge in disentangling whether the absence of a TBP in a tick is due to the absence of the pathogen in the previous host or influenced by the tick microbiome.

This sampling strategy also prevents assessing nesting dynamics of TBPs in the tick microbiome, as dosage and frequency of pathogen exposure are difficult to control in natural settings. To address the challenge associated with the limitations of questing tick sampling and better elucidate TBP-microbiome interactions in nature, we propose exploring systems involving ticks feeding on chronically infected hosts. In such scenarios, ticks would encounter a continuous presence of the pathogen across time. This approach, exemplified by systems such as *R. microplus* feeding on cattle chronically infected with *Anaplasma marginale*, provides an unique opportunity to assess the impact of the microbiome on TBP dynamics under more controlled conditions.

*Rhipicephalus microplus*, commonly known as the cattle tick or tropical cattle tick, is a significant ectoparasite that infests cattle and other livestock. This tick species holds great importance due to its capacity to transmit various pathogens, with *A. marginale* being one of the most notable (De La Fourniere et al., 2023; Pereira et al., 2022). *Anaplasma marginale* is a bacterium that causes bovine anaplasmosis, a disease characterized by anemia, fever, and other clinical symptoms in cattle (Salinas-Estrella et al., 2022). This disease can lead to significant economic losses due to decreased productivity, increased veterinary costs, and even livestock mortality (Rodríguez et al., 2009; Salinas-Estrella et al., 2022). *Rhipicephalus microplus* plays a crucial role in the transmission of *A. marginale*, as it acts as a vector by feeding on infected cattle and subsequently transmitting the bacteria to susceptible animals during subsequent feedings (Zivkovic et al., 2010). In a recent study, it was observed that while *A. marginale* consistently infected all cattle across different sampling periods, its presence was not uniformly detected in all *R. microplus* infesting the cattle (Piloto-Sardiñas, Foucault-Simonin, et al., 2023). This indicates nonlinearity between tick infestation rate and pathogen prevalence in ticks (Ostfeld & Keesing, 2023), likely under strong influence of the tick microbiome (Tonk-Rügen et al., 2023).

The present study aimed to investigate whether temporal fluctuations in the microbial communities of the cattle tick *R. microplus* could potentially disturb or alter the interactions between *A. marginale* and the tick microbial communities and as a consequence shape the impact of the pathogen on the microbial community assembly. To achieve this, *R. microplus* samples confirmed for infection with *A. marginale* were used (Piloto-Sardiñas, Foucault-Simonin, et al., 2023), and next-generation sequencing (NGS) was performed for microbiome characterization. Employing an experimental network approach, an *in silico* node removal technique was applied, a strategy previously employed to investigate the impact of *Rickettsia* pathogens on the microbiome assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks (Maitre et al., 2023). By simulating the absence of specific TBPs such as *Anaplasma* sp. in *silico*, the objective
was to evaluate the impact on clustering patterns, microbial composition and abundance of diverse taxa, community assembly, and network robustness over time. With this integrative approach, we aimed to uncover connections between *A. marginale* presence, temporal dynamics of the microbiome, and network structure, potentially identifying key microbial taxa to be used in anti-microbiota vaccines for the control of bovine anaplasmosis.

# 2 | MATERIALS AND METHODS

## 2.1 | Study design and tick samples

Tick samples collected from eight bovines on a farm in Mayabeque province, Cuba (Piloto-Sardiñas, Foucault-Simonin, et al., 2023) at three time points: July 2020 (J-20), September 2020 (S-20), and March 2021 (M-21) were included in this study. Engorged adult female ticks were manually collected from the same animals at different time points and morphologically identified as *R. microplus* using standardized taxonomic keys (Estrada-Peña et al., 2004; Piloto-Sardiñas, Foucault-Simonin, et al., 2023). Tick-borne pathogens (TBPs) were detected in individual tick samples through high-throughput real-time microfluidic PCR method (Piloto-Sardiñas, Foucault-Simonin, et al., 2023). This PCR method allows the detection of 27 bacterial species (belonging to the bacterial genera *Borrelia*, *Anaplasma*, *Ehrlichia*, *Rickettsia*, and *Mycoplasma*), 7 parasite species (such as *Babesia* and *Hepatozoon*), 5 bacterial genera, and 3 parasites taxa (Apicomplexa, *Theileria* and *Hepatozoon*) (Grech-Angelini et al., 2020; Michelet et al., 2014). The tested pathogens, target genes, and primer sequences used for amplification are shown in Table S1 (Gondard et al., 2020). Tick samples showing a low level of engorgement and single *A. marginale* infection (Piloto-Sardiñas, Foucault-Simonin, et al., 2023) were selected for microbiome sequencing.

Before DNA extraction, the collected ticks underwent a washing process, which involved two rounds of washing in miliQ sterile water and one round in 70% ethanol. It is worth noting that ethanol, rather than bleach, was used for washing to intentionally include both internal and external tick microbiome in our analysis, as we consider tick surface microbes to be part of the tick's microbiome. Following the washing process, the ticks were preserved in 70% ethanol and stored at -80°C until further processing. For the extraction of total DNA, the homogenization of whole ticks was performed on a MagNA Lyser instrument (Roche Molecular Diagnostics, Rotkreuz, Switzerland) at a speed of 5000rpm for 5 cycles of 60s each. Total DNA extraction was performed using the Wizard Genomic DNA Purification kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. The DNA samples were eluted in 60µL of DNA Rehydration Solution. The used of Colibri Microvolume Spectrophotometer (Titertek-Berthold, Pforzheim, Germany) allowed determining the quantitative and qualitative assessment of DNA extraction. Reagent extraction controls were set in DNA extraction process, using the same conditions as for the samples but using water as template. DNA amplification

was then performed on the extraction control in the same conditions as for any other sample.

## 2.2 | 16S rRNA amplicon sequencing and processing of raw sequences

A single lane of the Illumina MiSeq system was used to generate 251-base paired-end reads from variable region V4 of the 16S rRNA gene using barcoded universal primers (515F/806R) in ticks. The paired 16S rRNA raw sequences obtained from the J-20 (n=7), S-20 (n=7), and M-21 (n=8) samples were deposited in the SRA repository (Bioproject No. PRJNA1028823). Analysis of 16S rRNA sequences was performed using the Quantitative Insights into Microbial Ecology 2 (QIIME 2) pipeline (v. 2021.4) (Bolyen et al., 2019). The raw sequences (demultiplexed in fatsq files) were denoized, quality trimmed, and merged using the DADA2 software (Callahan et al., 2016) implemented in QIIME2 (Bolyen et al., 2019). The obtained amplicon sequence variants (ASVs) were aligned with q2-alignment of MAFFT (Katoh et al., 2002) and used to generate a phylogeny with q2-phylogeny of FastTree 2 (Price et al., 2010). Taxonomy was assigned to ASVs using a classify-sklearn naïve Bayes taxonomic classifier based on SILVA database (release 138) (Bokulich et al., 2018). Only the target sequence fragments were used for the classifier (i.e., the classifier was trained with primers 515F/806R) (Ren & Wu, 2016; Werner et al., 2012).

## 2.3 | Identification and removal of contaminants

The possible contaminants in the samples were statistically identified with the “Decontam” (Davis et al., 2018) package using the “prevalence” method. The method used compares the prevalence of each sequence feature in true samples to the prevalence in negative controls from the DNA extraction process to identify contaminants. Then, contaminants were removed from the dataset before downstream microbiome analysis (Davis et al., 2018).

## 2.4 | Microbial diversity, composition, and taxonomic differential relative abundance

To test the stability or variability of the microbiome over time, comparisons were made under three conditions: J-20, S-20, and M-21. To determine microbial diversity among the conditions, alpha and beta diversity metrics were calculated using q2-diversity plugin in QIIME 2 (Bolyen et al., 2019). Three alpha diversity metrics were explored using observed features (DeSantis et al., 2006) and Faith's phylogenetic diversity index (Faith, 1992) for richness, while evenness was explored with the Pielou's evenness index (Pielou, 1966). Differences in alpha-diversity metrics between groups were assessed with the Kruskal-Wallis test (p ≤ .05) using QIIME 2 (Bolyen et al., 2019). Beta-diversity was assessed with the Bray-Curtis dissimilarity index (Bray
was to evaluate the impact on clustering patterns, microbial composition and abundance of diverse taxa, community assembly, and network robustness over time. With this integrative approach, we aimed to uncover connections between *A. marginale* presence, temporal dynamics of the microbiome, and network structure, potentially identifying key microbial taxa to be used in anti-microbiota vaccines for the control of bovine anaplasmosis.

## 2 | MATERIALS AND METHODS

### 2.1 | Study design and tick samples

Tick samples collected from eight bovines on a farm in Mayabeque province, Cuba (Piloto-Sardiñas, Foucault-Simonin, et al., 2023) at three time points: July 2020 (J-20), September 2020 (S-20), and March 2021 (M-21) were included in this study. Engorged adult female ticks were manually collected from the same animals at different time points and morphologically identified as *R. microplus* using standardized taxonomic keys (Estrada-Peña et al., 2004; Piloto-Sardiñas, Foucault-Simonin, et al., 2023). Tick-borne pathogens (TBPs) were detected in individual tick samples through high-throughput real-time microfluidic PCR method (Piloto-Sardiñas, Foucault-Simonin, et al., 2023). This PCR method allows the detection of 27 bacterial species (belonging to the bacterial genera *Borrelia*, *Anaplasma*, *Ehrlichia*, *Rickettsia*, and *Mycoplasma*), 7 parasite species (such as *Babesia* and *Hepatozoon*), 5 bacterial genera, and 3 parasites taxa (Apicomplexa, *Theileria* and *Hepatozoon*) (Grech-Angelini et al., 2020; Michelet et al., 2014). The tested pathogens, target genes, and primer sequences used for amplification are shown in Table S1 (Gondard et al., 2020). Tick samples showing a low level of engorgement and single *A. marginale* infection (Piloto-Sardiñas, Foucault-Simonin, et al., 2023) were selected for microbiome sequencing.

Before DNA extraction, the collected ticks underwent a washing process, which involved two rounds of washing in miliQ sterile water and one round in 70% ethanol. It is worth noting that ethanol, rather than bleach, was used for washing to intentionally include both internal and external tick microbiome in our analysis, as we consider tick surface microbes to be part of the tick's microbiome. Following the washing process, the ticks were preserved in 70% ethanol and stored at -80°C until further processing. For the extraction of total DNA, the homogenization of whole ticks was performed on a MagNA Lyser instrument (Roche Molecular Diagnostics, Rotkreuz, Switzerland) at a speed of 5000rpm for 5 cycles of 60s each. Total DNA extraction was performed using the Wizard Genomic DNA Purification kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. The DNA samples were eluted in 60µL of DNA Rehydration Solution. The used of Colibri Microvolume Spectrophotometer (Titertek-Berthold, Pforzheim, Germany) allowed determining the quantitative and qualitative assessment of DNA extraction. Reagent extraction controls were set in DNA extraction process, using the same conditions as for the samples but using water as template. DNA amplification

was then performed on the extraction control in the same conditions as for any other sample.

### 2.2 | 16S rRNA amplicon sequencing and processing of raw sequences

A single lane of the Illumina MiSeq system was used to generate 251-base paired-end reads from variable region V4 of the 16S rRNA gene using barcoded universal primers (515F/806R) in ticks. The paired 16S rRNA raw sequences obtained from the J-20 (n=7), S-20 (n=7), and M-21 (n=8) samples were deposited in the SRA repository (Bioproject No. PRJNA1028823). Analysis of 16S rRNA sequences was performed using the Quantitative Insights into Microbial Ecology 2 (QIIME 2) pipeline (v. 2021.4) (Bolyen et al., 2019). The raw sequences (demultiplexed in fatsq files) were denoized, quality trimmed, and merged using the DADA2 software (Callahan et al., 2016) implemented in QIIME2 (Bolyen et al., 2019). The obtained amplicon sequence variants (ASVs) were aligned with q2-alignment of MAFFT (Katoh et al., 2002) and used to generate a phylogeny with q2-phylogeny of FastTree 2 (Price et al., 2010). Taxonomy was assigned to ASVs using a classify-sklearn naïve Bayes taxonomic classifier based on SILVA database (release 138) (Bokulich et al., 2018). Only the target sequence fragments were used for the classifier (i.e., the classifier was trained with primers 515F/806R) (Ren & Wu, 2016; Werner et al., 2012).

### 2.3 | Identification and removal of contaminants

The possible contaminants in the samples were statistically identified with the “Decontam” (Davis et al., 2018) package using the “prevalence” method. The method used compares the prevalence of each sequence feature in true samples to the prevalence in negative controls from the DNA extraction process to identify contaminants. Then, contaminants were removed from the dataset before downstream microbiome analysis (Davis et al., 2018).

### 2.4 | Microbial diversity, composition, and taxonomic differential relative abundance

To test the stability or variability of the microbiome over time, comparisons were made under three conditions: J-20, S-20, and M-21. To determine microbial diversity among the conditions, alpha and beta diversity metrics were calculated using q2-diversity plugin in QIIME 2 (Bolyen et al., 2019). Three alpha diversity metrics were explored using observed features (DeSantis et al., 2006) and Faith's phylogenetic diversity index (Faith, 1992) for richness, while evenness was explored with the Pielou's evenness index (Pielou, 1966). Differences in alpha-diversity metrics between groups were assessed with the Kruskal-Wallis test (p ≤ .05) using QIIME 2 (Bolyen et al., 2019). Beta-diversity was assessed with the Bray-Curtis dissimilarity index (Bray
& Curtis, 1957) with the PERMANOVA test ($p \le .01$) on QIIME 2. Beta dispersion was calculated using the betadisper function and the Vegan script implemented in RStudio (Oksanen et al., 2021), using an ANOVA test ($p \le .05$) as statistical analyses. Cluster analysis was performed with the Jaccard coefficient of similarity using Vegan (Oksanen et al., 2021) implemented in RStudio (RStudio Team, 2020). Unique and shared taxa among the three conditions were represented using Venn diagrams created with an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Differences in taxa relative abundance between the three conditions were tested using a Kruskal-Wallis test ($p \le .05$) and implemented using the ANOVA-Like Differential Expression (ALDEx2) package (Fernandes et al., 2013) on RStudio (RStudio Team, 2020). Only taxa with significant differences ($p \le .05$) were used for representation of the differential taxa relative abundance. Relative abundance was measured as centred log ratio (clr) transformation. The identified differentially abundant taxa were used to create a heatmap using the package "Heatplus" in RStudio (RStudio Team, 2020).

## 2.5 | Inference of bacterial co-occurrence networks

Co-occurrence networks were created for each dataset using the taxonomic profiles at family and genera level. The networks provide a graphical representation of the assembly of complex microbial communities within ecosystems. It allows us to analyze the associations that are established, as well as their nature. In the assembly, the nodes represent the taxa, while the edges represent the associations established between them. Analyses of significant positive (weight > 0.75) or negative (weight < -0.75) correlations were performed using the Sparse Correlations for Compositional data (SparCC) method (Friedman & Alm, 2012), implemented in RStudio (RStudio Team, 2020). Visualization and measurement of topological features (i.e., number of nodes and edges, network diameter, modularity, average degree, weighted degree, clustering coefficient and total count of basic undirected motifs of three fully connected vertices [triangles]) of the networks were performed using Gephi v0.10 (Bastian et al., 2009).

With the aim of identifying microbial taxa shared for the conditions, a Core Association Networks (CAN) were created for J-20/S20, S20/M21, and J-20/M21, using a software toolbox, anuran (a toolbox with null models for identification of nonrandom patterns in association networks) (Röttjers et al., 2021), and this version was tested in Python 3.6.

## 2.6 | Keystone taxa identification

Keystone taxa were identified within the community for each of the condition, based on three criteria, as previously reported (Mateos-Hernández et al., 2021): (i) ubiquitousness (microbial taxa present in all samples in an experimental group), (ii) eigenvector centrality

higher than 0.75, and (iii) high mean relative abundance (i.e., higher than that of the mean relative abundance of all taxa in an experimental group). Additionally, the common keystone taxa for J-20, S-20, and M-21 microbial community were identified.

## 2.7 | Local connectivity of *Anaplasma* in the microbial community

To explore the role of *Anaplasma* within the community, its direct relationship with the rest of the bacterial microbiome was determined. For this purpose, subnetworks were constructed where *Anaplasma* was visualized with its direct positive and negative associations. The analyses were carried out in Gephi v0.10 (Bastian et al., 2009), and the strength of the edges was presented with the SparCC weight.

## 2.8 | Analysis of centrality measures distribution in network nodes

The topology of the taxa in the network was analyzed with two connectivity types: (i) within-module connectivity ($Z_i$), which describes how the taxon is connected to others within its module and (ii) among-module connectivity ($P_i$), which describes the taxa connectivity with other taxa in different modules (Guimera & Nunes Amaral, 2005). The taxa are divided into four categories: (i) peripherals taxa ($Z_i \le 2.5$ and $P_i \le 0.62$), which contain taxa with few edges in and out of its module; (ii) connectors ($Z_i \le 2.5$ and $P_i > 0.62$), which contain taxa connected to other modules than its own; (iii) module hubs ($Z_i > 2.5$ and $P_i \le 0.62$), which contain taxa highly connected with members of their own module; and (iv) network hubs ($Z_i > 2.5$ and $P_i > 0.62$), which contain taxa highly connected with members within and among its module. For each taxon, $Z_i$ and $P_i$ values were calculated using only positive edges, with the R package "code-zi-pi-plot" described by (Cao et al., 2018) and (Guo et al., 2022) in Rstudio (R studio Team, 2020) and visualized with GraphPad Prism version 8.0.1 (GraphPad Software, San Diego, California USA).

## 2.9 | Differential network analysis and modules composition

With the aim of comparing the correlations between the same taxa in two different bacterial networks, a statistical network estimation analysis was performed using the network construction and comparison for microbiome (NetCoMi) method (Peschel et al., 2021) implemented in RStudio (RStudio Team, 2020). The comparison was carried out: (i) with *Anaplasma* (wA) and (ii) without *Anaplasma* (woA) (in-silico removal) for each time point, (iii) wA vs. woA in the same time point. To test for dissimilarities between the two networks [i.e., J-20 vs. S-20; J-20 (wA) vs. J-20 (woA)], the Jaccard index was calculated to test for dissimilarities between nodes in the two networks
& Curtis, 1957) with the PERMANOVA test ($p \le .01$) on QIIME 2. Beta dispersion was calculated using the betadisper function and the Vegan script implemented in RStudio (Oksanen et al., 2021), using an ANOVA test ($p \le .05$) as statistical analyses. Cluster analysis was performed with the Jaccard coefficient of similarity using Vegan (Oksanen et al., 2021) implemented in RStudio (RStudio Team, 2020). Unique and shared taxa among the three conditions were represented using Venn diagrams created with an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Differences in taxa relative abundance between the three conditions were tested using a Kruskal-Wallis test ($p \le .05$) and implemented using the ANOVA-Like Differential Expression (ALDEx2) package (Fernandes et al., 2013) on RStudio (RStudio Team, 2020). Only taxa with significant differences ($p \le .05$) were used for representation of the differential taxa relative abundance. Relative abundance was measured as centred log ratio (clr) transformation. The identified differentially abundant taxa were used to create a heatmap using the package "Heatplus" in RStudio (RStudio Team, 2020).

## 2.5 | Inference of bacterial co-occurrence networks

Co-occurrence networks were created for each dataset using the taxonomic profiles at family and genera level. The networks provide a graphical representation of the assembly of complex microbial communities within ecosystems. It allows us to analyze the associations that are established, as well as their nature. In the assembly, the nodes represent the taxa, while the edges represent the associations established between them. Analyses of significant positive (weight > 0.75) or negative (weight < -0.75) correlations were performed using the Sparse Correlations for Compositional data (SparCC) method (Friedman & Alm, 2012), implemented in RStudio (RStudio Team, 2020). Visualization and measurement of topological features (i.e., number of nodes and edges, network diameter, modularity, average degree, weighted degree, clustering coefficient and total count of basic undirected motifs of three fully connected vertices [triangles]) of the networks were performed using Gephi v0.10 (Bastian et al., 2009).

With the aim of identifying microbial taxa shared for the conditions, a Core Association Networks (CAN) were created for J-20/S20, S20/M21, and J-20/M21, using a software toolbox, anuran (a toolbox with null models for identification of nonrandom patterns in association networks) (Röttjers et al., 2021), and this version was tested in Python 3.6.

## 2.6 | Keystone taxa identification

Keystone taxa were identified within the community for each of the condition, based on three criteria, as previously reported (Mateos-Hernández et al., 2021): (i) ubiquitousness (microbial taxa present in all samples in an experimental group), (ii) eigenvector centrality

higher than 0.75, and (iii) high mean relative abundance (i.e., higher than that of the mean relative abundance of all taxa in an experimental group). Additionally, the common keystone taxa for J-20, S-20, and M-21 microbial community were identified.

## 2.7 | Local connectivity of *Anaplasma* in the microbial community

To explore the role of *Anaplasma* within the community, its direct relationship with the rest of the bacterial microbiome was determined. For this purpose, subnetworks were constructed where *Anaplasma* was visualized with its direct positive and negative associations. The analyses were carried out in Gephi v0.10 (Bastian et al., 2009), and the strength of the edges was presented with the SparCC weight.

## 2.8 | Analysis of centrality measures distribution in network nodes

The topology of the taxa in the network was analyzed with two connectivity types: (i) within-module connectivity ($Z_i$), which describes how the taxon is connected to others within its module and (ii) among-module connectivity ($P_i$), which describes the taxa connectivity with other taxa in different modules (Guimera & Nunes Amaral, 2005). The taxa are divided into four categories: (i) peripherals taxa ($Z_i \le 2.5$ and $P_i \le 0.62$), which contain taxa with few edges in and out of its module; (ii) connectors ($Z_i \le 2.5$ and $P_i > 0.62$), which contain taxa connected to other modules than its own; (iii) module hubs ($Z_i > 2.5$ and $P_i \le 0.62$), which contain taxa highly connected with members of their own module; and (iv) network hubs ($Z_i > 2.5$ and $P_i > 0.62$), which contain taxa highly connected with members within and among its module. For each taxon, $Z_i$ and $P_i$ values were calculated using only positive edges, with the R package "code-zi-pi-plot" described by (Cao et al., 2018) and (Guo et al., 2022) in Rstudio (R studio Team, 2020) and visualized with GraphPad Prism version 8.0.1 (GraphPad Software, San Diego, California USA).

## 2.9 | Differential network analysis and modules composition

With the aim of comparing the correlations between the same taxa in two different bacterial networks, a statistical network estimation analysis was performed using the network construction and comparison for microbiome (NetCoMi) method (Peschel et al., 2021) implemented in RStudio (RStudio Team, 2020). The comparison was carried out: (i) with *Anaplasma* (wA) and (ii) without *Anaplasma* (woA) (in-silico removal) for each time point, (iii) wA vs. woA in the same time point. To test for dissimilarities between the two networks [i.e., J-20 vs. S-20; J-20 (wA) vs. J-20 (woA)], the Jaccard index was calculated to test for dissimilarities between nodes in the two networks
for degree, betweenness centrality, closeness centrality and eigen-
vector centrality. The Jaccard index tests for the similarity between
sets of "most central nodes" of networks, which are defined as those
nodes with a centrality value above the empirical 75% quartile. This
index expresses the similarity of the sets of most central nodes as
well as the sets of hub taxa between the two networks. The Jaccard
index ranges from 0 (completely different sets) to 1 (sets equal).
The two p-values *p* (*J*≤*j*) and *p* (*J*≥*j*) for each Jaccard index are the
probability that the observed value of Jaccard's index is "less than or
equal" or "higher than or equal," respectively, to the Jaccard value
expected at random which is calculated taking into account the pre-
sent total number of taxa in both sets (Real & Vargas, 1996). The ARI
was calculated to test the dissimilarity of clustering in the networks.
The ARI values range from -1 to 1. Negative and positive ARI values
mean lower and higher than random clustering, respectively. An ARI
value of 1 corresponds to identical clustering and 0 to dissimilar clus-
tering. The p-value tests whether the calculated value is significantly
different from zero (Peschel et al., 2021).

To assess the potential direct and indirect consequences of
removing *Anaplasma* from the networks, we focused on two mod-
ules (M1 and M2) within the networks. M1 represented the mod-
ule containing *Anaplasma*, while M2 exhibited a higher number of
taxa and modularity value. This ensured the equivalence of modules
between wA and woA networks. Subsequently, subnetworks were
constructed for comparison (wA vs. woA) at the same time point, to
better understand the network dynamics.

## 2.10 | Network robustness analysis in nodes removal and addition

The robustness of the networks against disturbances due to removal and addition of nodes was determined. In the analysis of node removal, the proportion of eliminated nodes necessary to achieve a connectivity loss of 0.40 (40%) and 0.80 (80%) was recorded, after directed and random attacks. Two scenarios were evaluated: (i) robustness of the networks at each time point (J-20, S-20 and M-21) and (ii) robustness of the (wA-woA) networks at the same time point. The robustness of the networks was calculated using the Network Strengths and Weaknesses Analysis (NetSwan) package (Lhomme, 2015) in RStudio (RStudio Team, 2020).

The robustness of the networks at each time point (J-20, S-20, and M-21) and at the same time point (wA vs. woA) was analyzed for the addition of nodes, using the Network analysis and visualization package (Freitas et al., 2021). Nodes were incrementally added in sections ranging from 5 to 100, and network connectivity was measured based on the degree metric of the largest connected component (LCC) and average path length. A Wilcoxon signed-rank test was conducted to calculate *p*-values for LCC and average path length. The *p*-values were adjusted using the Benjamini–Hochberg (BH) method to control the false discovery rate. Additionally, bootstrapping was performed to obtain confidence intervals for the variables. Significance was determined at a threshold of *p* < .05.

3 | RESULTS

3.1 | Diversity, composition, and abundance of bacterial taxa in R. *microplus* microbiome over time

Diversity, composition, and abundance of bacterial taxa in the R. microplus microbiome were assessed over two consecutive years, 2020 and 2021, using 16S rRNA gene profiling after statistical identification and removal of DNA features identified as contaminants (Table S2). Differences in α-diversity were significant, with higher observed features in ticks collected in J-20 compared to those collected in M-21 (Kruskal-Wallis, *p* = .049, Figure 1a). Faith's phylogenetic diversity (Faith's PD) also differed between M-21 and J-20 (Kruskal-Wallis, *p* = .015, Figure 1b), as well as M-21 and S-20 (Kruskal-Wallis, *p* = .021, Figure 1b). However, no significant differences were observed between J-20 and S-20 for observed features (Kruskal-Wallis, *p* > .05, Figure 1a) or Faith's PD (Kruskal-Wallis, *p* > .05, Figure 1b) metrics. Evenness showed no significant changes over time (*p* > .05, Figure 1c).

Bray–Curtis index analysis revealed no significant differences
in microbiome composition between J-20 and S-20 (PERMANOVA,
*p* > .01), but both differed from M-21 (PERMANOVA, *p* = .001,
*F* = 2.652). Beta dispersion showed no significant within-group
variability (ANOVA test, *p* > .05, Figure 1d). Jaccard clustering con-
firmed two distinct clusters: one with all M-21 samples and two
with J-20 and S-20 samples (Figure 1e), aligning with Bray–Curtis
index findings. Compositional analysis identified 437 bacterial
taxa, with 67.7% shared across all samples (Figure 1f). Unique taxa
were found in M-21 (1.60%) and J-20 (0.92%), while S-20 had none
(Figure 1f, Table S3).

Differential relative abundance analysis identified significant
changes in 20 taxa across the three conditions (Figure 1g, Table S4).
Murdochiella, Neisseria, and Rickettsiales were more abundant
in J-20 and S-20, while Anaplasma, Cloacibacterium, Delftia, and
Frigoribacterium were higher in M-21. Anaplasma consistently showed
significantly higher mean relative abundance in M-21 (10.8±1.41)
compared to J-20 (3.01±0.78) and S-20 (4.36±0.93) (*p* < .05), align-
ing with previous PCR-confirmed findings (Piloto-Sardiñas, Foucault-
Simonin, et al., 2023). The increased relative abundance of Anaplasma,
coupled with reduced microbial diversity over time, suggest a poten-
tial interaction with R. microplus microbial communities.

## 3.2 | Dynamics of *Anaplasma* nesting in the microbial communities of *R. microplus*

The dynamics of *Anaplasma* nesting within the microbial communi-
ties of *R. microplus* were investigated using co-occurrence networks
to assess community assembly over time. Notably, J-20 displayed the
most total and connected nodes, while M-21 had the least, show-
casing topological variations (Figure 2a-c, Table 1). Additionally,
J-20 displayed the highest number of correlations with a balanced
positive-negative ratio (Figure 2a, Table 1). In contrast, S-20 and
for degree, betweenness centrality, closeness centrality and eigen-
vector centrality. The Jaccard index tests for the similarity between
sets of "most central nodes" of networks, which are defined as those
nodes with a centrality value above the empirical 75% quartile. This
index expresses the similarity of the sets of most central nodes as
well as the sets of hub taxa between the two networks. The Jaccard
index ranges from 0 (completely different sets) to 1 (sets equal).
The two p-values *p* (*J*≤*j*) and *p* (*J*≥*j*) for each Jaccard index are the
probability that the observed value of Jaccard's index is "less than or
equal" or "higher than or equal," respectively, to the Jaccard value
expected at random which is calculated taking into account the pre-
sent total number of taxa in both sets (Real & Vargas, 1996). The ARI
was calculated to test the dissimilarity of clustering in the networks.
The ARI values range from -1 to 1. Negative and positive ARI values
mean lower and higher than random clustering, respectively. An ARI
value of 1 corresponds to identical clustering and 0 to dissimilar clus-
tering. The p-value tests whether the calculated value is significantly
different from zero (Peschel et al., 2021).

To assess the potential direct and indirect consequences of
removing *Anaplasma* from the networks, we focused on two mod-
ules (M1 and M2) within the networks. M1 represented the mod-
ule containing *Anaplasma*, while M2 exhibited a higher number of
taxa and modularity value. This ensured the equivalence of modules
between wA and woA networks. Subsequently, subnetworks were
constructed for comparison (wA vs. woA) at the same time point, to
better understand the network dynamics.

## 2.10 | Network robustness analysis in nodes removal and addition

The robustness of the networks against disturbances due to removal and addition of nodes was determined. In the analysis of node removal, the proportion of eliminated nodes necessary to achieve a connectivity loss of 0.40 (40%) and 0.80 (80%) was recorded, after directed and random attacks. Two scenarios were evaluated: (i) robustness of the networks at each time point (J-20, S-20 and M-21) and (ii) robustness of the (wA-woA) networks at the same time point. The robustness of the networks was calculated using the Network Strengths and Weaknesses Analysis (NetSwan) package (Lhomme, 2015) in RStudio (RStudio Team, 2020).

The robustness of the networks at each time point (J-20, S-20, and M-21) and at the same time point (wA vs. woA) was analyzed for the addition of nodes, using the Network analysis and visualization package (Freitas et al., 2021). Nodes were incrementally added in sections ranging from 5 to 100, and network connectivity was measured based on the degree metric of the largest connected component (LCC) and average path length. A Wilcoxon signed-rank test was conducted to calculate *p*-values for LCC and average path length. The *p*-values were adjusted using the Benjamini–Hochberg (BH) method to control the false discovery rate. Additionally, bootstrapping was performed to obtain confidence intervals for the variables. Significance was determined at a threshold of *p* < .05.

3 | RESULTS

## 3.1 | Diversity, composition, and abundance of bacterial taxa in R. *microplus* microbiome over time

Diversity, composition, and abundance of bacterial taxa in the R. microplus microbiome were assessed over two consecutive years, 2020 and 2021, using 16S rRNA gene profiling after statistical identification and removal of DNA features identified as contaminants (Table S2). Differences in α-diversity were significant, with higher observed features in ticks collected in J-20 compared to those collected in M-21 (Kruskal-Wallis, *p* = .049, Figure 1a). Faith's phylogenetic diversity (Faith's PD) also differed between M-21 and J-20 (Kruskal-Wallis, *p* = .015, Figure 1b), as well as M-21 and S-20 (Kruskal-Wallis, *p* = .021, Figure 1b). However, no significant differences were observed between J-20 and S-20 for observed features (Kruskal-Wallis, *p* > .05, Figure 1a) or Faith's PD (Kruskal-Wallis, *p* > .05, Figure 1b) metrics. Evenness showed no significant changes over time (*p* > .05, Figure 1c).

Bray-Curtis index analysis revealed no significant differences
in microbiome composition between J-20 and S-20 (PERMANOVA,
*p* > .01), but both differed from M-21 (PERMANOVA, *p* = .001,
*F* = 2.652). Beta dispersion showed no significant within-group
variability (ANOVA test, *p* > .05, Figure 1d). Jaccard clustering con-
firmed two distinct clusters: one with all M-21 samples and two
with J-20 and S-20 samples (Figure 1e), aligning with Bray-Curtis
index findings. Compositional analysis identified 437 bacterial
taxa, with 67.7% shared across all samples (Figure 1f). Unique taxa
were found in M-21 (1.60%) and J-20 (0.92%), while S-20 had none
(Figure 1f, Table S3).

Differential relative abundance analysis identified significant
changes in 20 taxa across the three conditions (Figure 1g, Table S4).
Murdochiella, Neisseria, and Rickettsiales were more abundant
in J-20 and S-20, while Anaplasma, Cloacibacterium, Delftia, and
Frigoribacterium were higher in M-21. Anaplasma consistently showed
significantly higher mean relative abundance in M-21 (10.8±1.41)
compared to J-20 (3.01±0.78) and S-20 (4.36±0.93) (*p* < .05), align-
ing with previous PCR-confirmed findings (Piloto-Sardiñas, Foucault-
Simonin, et al., 2023). The increased relative abundance of Anaplasma,
coupled with reduced microbial diversity over time, suggest a poten-
tial interaction with R. microplus microbial communities.

## 3.2 | Dynamics of *Anaplasma* nesting in the microbial communities of *R. microplus*

The dynamics of *Anaplasma* nesting within the microbial communi-
ties of *R. microplus* were investigated using co-occurrence networks
to assess community assembly over time. Notably, J-20 displayed the
most total and connected nodes, while M-21 had the least, show-
casing topological variations (Figure 2a-c, Table 1). Additionally,
J-20 displayed the highest number of correlations with a balanced
positive-negative ratio (Figure 2a, Table 1). In contrast, S-20 and
<figure><img src="image_2.png" /><figcaption>FIGURE 1 Comparison of diversity of complex microbial communities within *Rhipicephalus microplus* over time. Comparison of alpha diversity between J-20, S-20 and M-21 (Kruskal-Wallis test, significant differences for *p* ≤ .05), (a) observed features, (b) Faith's phylogenetic diversity (PD), and (c) Pielou's evenness index. (d) Comparison of beta - diversity with Bray Curtis dissimilarity index between J-20, S-20, and M-21. Beta dispersion of three sets of samples (global comparison). Small circles, crosses and triangles represent samples, and ellipses represent centroid position for each group. This test use principal coordinate analysis (PCoA), it is used to explore and to visualize variability in a microbial community. ANOVA test was performed and showed that beta dispersion of the three sets of samples (three conditions) is not significantly different (*p* = .51). (e) Jaccard clusterisation of the tick samples collected in J-20, S-20, and M-21. The samples are represented by circles and the groups by colors (legend). (f) Venn Diagram displaying the comparison of taxa composition in ticks collected at the three sampling times. Common and unique taxa between the conditions are represented. (g) Comparison of relative abundance of complex microbial communities within *R. microplus* over time. 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 per group. Each line represents bacterial taxa with significant changes between the datasets. Color represent the clr value (range from -15 to 15).</figcaption></figure>

M-21 showed greater differences in positive-negative associa-
tions (Figure 2b,c, Table 1). Despite M-21 having fewer connected
nodes, it had a high proportion of positive associations (Figure 2c,
Table 1). S-20 and M-21 displayed lower modularity than J-20. The
three networks had similar diameter values between them (Table 1).
The total count of motifs over time was determined (Table 1). The
J-20 network presented the highest total number of motifs (13,322),
followed by S-20 (519), while M-21 (18) had the lowest (Table 1).

The high modularity values and a considerable number of motifs
within J-20 (Table 1) indicate a strong connection between the ver-
tices (nodes) compared to the S-20 and M-21 networks (Table S5).
Compositional analysis revealed that J-20 had the higher number of
unique nodes (69), followed by S-20 (51) and M-21 (12) (Figure 2d,
Table S6). The significant overlap in nodes between J-20 and S-20
compared to M-21 implies minimal taxa variability and minimal as-
sembly variation between J-20 and S-20 (Figure 2a-d).
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** Comparison of diversity of complex microbial communities within *Rhipicephalus microplus* over time. Comparison of alpha diversity between J-20, S-20 and M-21 (Kruskal-Wallis test, significant differences for *p* ≤ .05), (a) observed features, (b) Faith's phylogenetic diversity (PD), and (c) Pielou's evenness index. (d) Comparison of beta-diversity with Bray Curtis dissimilarity index between J-20, S-20, and M-21. Beta dispersion of three sets of samples (global comparison). Small circles, crosses and triangles represent samples, and ellipses represent centroid position for each group. This test uses principal coordinate analysis (PCoA), it is used to explore and to visualize variability in a microbial community. ANOVA test was performed and showed that beta dispersion of the three sets of samples (three conditions) is not significantly different (*p* = .51). (e) Jaccard clusterisation of the tick samples collected in J-20, S-20, and M-21. The samples are represented by circles and the groups by colors (legend). (f) Venn Diagram displaying the comparison of taxa composition in ticks collected at the three sampling times. Common and unique taxa between the conditions are represented. (g) Comparison of relative abundance of complex microbial communities within *R. microplus* over time. 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 per group. Each line represents bacterial taxa with significant changes between the datasets. Color represents the clr value (range from -15 to 15).</figcaption></figure>

M-21 showed greater differences in positive-negative associations (Figure 2b,c, Table 1). Despite M-21 having fewer connected nodes, it had a high proportion of positive associations (Figure 2c, Table 1). S-20 and M-21 displayed lower modularity than J-20. The three networks had similar diameter values between them (Table 1). The total count of motifs over time was determined (Table 1). The J-20 network presented the highest total number of motifs (13,322), followed by S-20 (519), while M-21 (18) had the lowest (Table 1).

The high modularity values and a considerable number of motifs within J-20 (Table 1) indicate a strong connection between the vertices (nodes) compared to the S-20 and M-21 networks (Table S5). Compositional analysis revealed that J-20 had the higher number of unique nodes (69), followed by S-20 (51) and M-21 (12) (Figure 2d, Table S6). The significant overlap in nodes between J-20 and S-20 compared to M-21 implies minimal taxa variability and minimal assembly variation between J-20 and S-20 (Figure 2a-d).
<figure><img src="image_3.png" /><figcaption>**FIGURE 2** Dynamics of *Anaplasma* nesting in the microbial communities of *Rhipicephalus microplus* over time. Global and local co-occurrence networks. Co-occurrence networks of (a) J-20, (b) S-20, and (c) M-21 networks. Node colors are based on modularity class metric and equal color means modules of co-occurring taxa. The size of the nodes is proportional to the eigenvector centrality of each taxon. The colors in the edges represent strong positive (blue) or negative (red) correlations (SparCC >0.75 or <-0.75). (d) Venn diagram displaying the comparison of networks composition. CAN for: (e) J-20/ S-20, (f) S-20/M-21, and (g) J-20/M21, the colors in the edges represent strong positive (blue) or negative (red) correlations and the nodes represent bacterial taxa. Sub-networks of local connectivity and indirect association of *Anaplasma* with keystone bacteria within the bacterial community in *R. microplus* over time: *Anaplasma*'s local connectivity within (h) J-20 and (i) S-20 networks, (j) *Anaplasma*'s indirect association with keystone bacteria in S-20 network and (k) *Anaplasma*'s local connectivity within M-21 network. Within-module and among-module connectivities, $Z_i-P_i$ plot of the individual genera from three groups: (l) J-20, (m) S-20, and (n) M-21. (o) Venn diagram displaying the comparison of the connectors between the S-20 and M-21 networks.</figcaption></figure>

The core association network (CAN) revealed 22 core associated nodes between J-20 and S-20 networks, supporting the observation of minimal variability in the assembly of the microbial community in these two conditions (Figure 2e). In contrast, four and three core associated nodes were found between M-21 and S-20 (Figure 2f), and J-20 networks (Figure 2g), respectively. Interestingly, *Anaplasma* was absent in the cores shared by J-20/S-20 and J-20/M-21 (Figure 2e,f), but present in the S-20/M-21 core, establishing a positive association with JS1 (Caldribacteriota) and Deferribacteraceae (Figure 2g). The presence of *Anaplasma* as a core associated node suggests that, alongside its increased abundance in the community, this pathogen may acquire importance in the community across time.

Analyzing the local connectivity of *Anaplasma* across the three specified time points revealed dynamic relationships. In the J-20 network, *Anaplasma* displayed numerous direct associations with nodes exhibiting high eigenvector centrality values (*Pajaroellobacter*, *KD4-96*, *Bacteroidales*, *Deferribacteraceae*) compared to other nodes in the network (Figure 2h). However, despite this, their mean clr values were lower than the mean relative abundance of all taxa in the experimental condition, disqualifying them as keystone taxa. In the S-20 network, *Anaplasma* mostly exhibited positive associations, except for a negative association with *Rickettsiella* (Figure 2i). Notably, in this condition, *Anaplasma* had higher eigenvector centrality compared to J-20, with a mean relative abundance value below the average relative abundance
<figure><img src="image_5.png" /></figure>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

FIGURE 2 Dynamics of *Anaplasma* nesting in the microbial communities of *Rhipicephalus microplus* over time. Global and local co-occurrence networks. Co-occurrence networks of (a) J-20, (b) S-20, and (c) M-21 networks. Node colors are based on modularity class metric and equal color means modules of co-occurring taxa. The size of the nodes is proportional to the eigenvector centrality of each taxon. The colors in the edges represent strong positive (blue) or negative (red) correlations (SparCC >0.75 or <-0.75). (d) Venn diagram displaying the comparison of networks composition. CAN for: (e) J-20/ S-20, (f) S-20/M-21, and (g) J-20/M21, the colors in the edges represent strong positive (blue) or negative (red) correlations and the nodes represent bacterial taxa. Sub-networks of local connectivity and indirect association of *Anaplasma* with keystone bacteria within the bacterial community in *R. microplus* over time: *Anaplasma*'s local connectivity within (h) J-20 and (i) S-20 networks, (j) *Anaplasma*'s indirect association with keystone bacteria in S-20 network and (k) *Anaplasma*'s local connectivity within M-21 network. Within-module and among-module connectivities, $Z_i-P_i$ plot of the individual genera from three groups: (l) J-20, (m) S-20, and (n) M-21. (o) Venn diagram displaying the comparison of the connectors between the S-20 and M-21 networks.

The core association network (CAN) revealed 22 core associated nodes between J-20 and S-20 networks, supporting the observation of minimal variability in the assembly of the microbial community in these two conditions (Figure 2e). In contrast, four and three core associated nodes were found between M-21 and S-20 (Figure 2f), and J-20 networks (Figure 2g), respectively. Interestingly, *Anaplasma* was absent in the cores shared by J-20/S-20 and J-20/M-21 (Figure 2e,f), but present in the S-20/M-21 core, establishing a positive association with JS1 (Caldribacteriota) and Deferribacteraceae (Figure 2g). The presence of *Anaplasma* as a core associated node suggests that, alongside its increased abundance in the community, this pathogen may acquire importance in the community across time.

Analyzing the local connectivity of *Anaplasma* across the three specified time points revealed dynamic relationships. In the J-20 network, *Anaplasma* displayed numerous direct associations with nodes exhibiting high eigenvector centrality values (*Pajaroellobacter*, *KD4-96*, *Bacteroidales*, *Deferribacteraceae*) compared to other nodes in the network (Figure 2h). However, despite this, their mean clr values were lower than the mean relative abundance of all taxa in the experimental condition, disqualifying them as keystone taxa. In the S-20 network, *Anaplasma* mostly exhibited positive associations, except for a negative association with *Rickettsiella* (Figure 2i). Notably, in this condition, *Anaplasma* had higher eigenvector centrality compared to J-20, with a mean relative abundance value below the average relative abundance
<table><caption>TABLE 1 Topological features of J-20, S-20, and M-21 networks with (wA) and without (woA) *Anaplasma* for each time.</caption><thead><tr><th>Topological features</th><th>J-20 (wA)</th><th>S-20 (wA)</th><th>M-21 (wA)</th><th>J-20 (woA)</th><th>S-20 (woA)</th><th>M-21 (woA)</th></tr></thead><tbody><tr><td>Total nodes</td><td>401</td><td>389</td><td>369</td><td>400</td><td>388</td><td>368</td></tr><tr><td>Connected nodes</td><td>240</td><td>215</td><td>90</td><td>231</td><td>195</td><td>91</td></tr><tr><td>Edges</td><td>975</td><td>369</td><td>76</td><td>867</td><td>375</td><td>72</td></tr><tr><td>Positive correlations</td><td>514 (52.7%)</td><td>232 (62.9%)</td><td>48 (60.5%)</td><td>493 (56.9%)</td><td>222 (59.2%)</td><td>44 (61.1%)</td></tr><tr><td>Negative correlations</td><td>461 (47.3%)</td><td>137 (37.1%)</td><td>28 (36.8%)</td><td>374 (43.1%)</td><td>153 (40.8%)</td><td>28 (38.9%)</td></tr><tr><td>Modularity</td><td>6.43</td><td>2.09</td><td>2.01</td><td>3.24</td><td>2.61</td><td>2.46</td></tr><tr><td>Network diameter</td><td>12</td><td>15</td><td>10</td><td>12</td><td>14</td><td>14</td></tr><tr><td>Average degree</td><td>8.13</td><td>3.43</td><td>1.69</td><td>7.5</td><td>3.85</td><td>1.58</td></tr><tr><td>Weighted degree</td><td>0.40</td><td>0.71</td><td>0.36</td><td>0.86</td><td>0.57</td><td>0.27</td></tr><tr><td>Clustering coefficient</td><td>0.52</td><td>0.39</td><td>0.29</td><td>0.46</td><td>0.38</td><td>0.21</td></tr><tr><td>Number of motifs</td><td>13,322</td><td>519</td><td>18</td><td>9549</td><td>561</td><td>12</td></tr></tbody></table>

<table><caption>TABLE 2 Keystone taxa of the bacterial communities within *R. microplus* by condition and in common.</caption><thead><tr><th>Condition</th><th>Keystone taxa by condition and in common</th></tr></thead><tbody><tr><td>July 2020</td><td>TRA3-20<br>Flavonifractor<br>Rhizobiaceae uncultured<br>Guggenheimella<br>Hydrocarboniphaga<br>Ga0077536<br>Ideonella<br>Nakamurella<br>Acidibacter<br>Staphylococcaceae<br>Parvibacter<br>Aquicella</td></tr><tr><td>September 2020</td><td>Ga0077536<br>Murdochielia<br>Acetobacteriaceae<br>Solimonas</td></tr><tr><td>March 2021</td><td>Quadrisphaera<br>Acidibacter<br>Anaplasma<br>Cellulomonadaceae</td></tr><tr><td>July 2020–September 2020</td><td>Ga0077536</td></tr><tr><td>July 2020–March 2021</td><td>Acidibacter</td></tr></tbody></table>

of all taxa in S-20. This node indirectly linked with three keystone taxa (Ga0077536, *Murdochielia*, *Solimonas*) (Table 2), through a direct association with *Haliangium* (Figure 2j). Notably, in the M-21 network, *Anaplasma* appeared as a keystone taxon (Table 2), with positive association module with seven nodes (Figure 2k). In the M-21 network, *Rickettsiella* was no longer part of the network, possibly due to co-exclusion caused by a negative interaction observed with *Anaplasma* in S20 (Figure 2i). Most direct associations of *Anaplasma* at a given time point were either displaced or did not co-occur with the taxon in the rest of the conditions, except for *Deferribacteraceae* and *JS1* (Caldribacteriota) in J-20 and M-21 (Table S7).

The evaluation of within-module ($Z_i$) and among-module ($P_i$) connectivity revealed similar distributions of nodes across the three conditions, with most nodes as peripheral (low $Z_i$ and $P_i$) and no network hubs identified (high $Z_i$ and $P_i$; Figure 2l-n). However, the presence of module hubs differed among the three groups, with nine, two and none taxa considered as module hubs in M-21, S-20, and J-20 networks, respectively (Figure 2l-n). The analysis revealed 24 and 14 unique connectors in S-20 and M-21 networks, respectively. Propionibacteriaceae and Lachnospiraceae UCG-002 promote coherence of both networks (Figure 2o, Table S8). *Anaplasma*'s position within the networks varied; it was peripheral in J-20 and S-20 networks (Figure 2l,m), but contributed to module coherence in the M-21 network (Figure 2n).

*Anaplasma* exhibited varying associations, evolving from diverse connections in J-20 to higher centrality and positive associations in S-20. Remarkably, in M-21, *Anaplasma* emerged as a keystone taxon, forming an independent cluster nested in the network establishing positives associations with other taxa. The absence of certain taxa and the potential co-exclusion of *Rickettsiella* may underscore the distinctive impact of *Anaplasma* on the community's dynamics over time.

### 3.3 | Influence of *Anaplasma* on the assembly, and hierarchy of *R. microplus* microbiome over time

To investigate *Anaplasma*'s impact on community assembly, network topology after its removal (woA) was analyzed and compared to the network with *Anaplasma* (wA). The removal of *Anaplasma* led to the loss of nodes in both the J-20 (woA) and S-20 (woA) networks, resulting in the depletion of both positive and negative correlations (Table 1, Figure S1A,B). Notably, within S-20 (woA) network, the associations changed with 16 new negative correlations appearing and 10 positive correlations being lost, suggesting that *Anaplasma* exerts an influence on interactions within the microbial community (Table 1, Figure S1B). In M-21 (woA), only one new connection formed within the network, and four positive correlations were lost (Table 1, Figure S1C). The removal of *Anaplasma* resulted in a decrease in the total number of motifs across all three networks, with the modularity
<table><caption>TABLE 1 Topological features of J-20, S-20, and M-21 networks with (wA) and without (woA) *Anaplasma* for each time.</caption><thead><tr><th>Topological features</th><th>J-20 (wA)</th><th>S-20 (wA)</th><th>M-21 (wA)</th><th>J-20 (woA)</th><th>S-20 (woA)</th><th>M-21 (woA)</th></tr></thead><tbody><tr><td>Total nodes</td><td>401</td><td>389</td><td>369</td><td>400</td><td>388</td><td>368</td></tr><tr><td>Connected nodes</td><td>240</td><td>215</td><td>90</td><td>231</td><td>195</td><td>91</td></tr><tr><td>Edges</td><td>975</td><td>369</td><td>76</td><td>867</td><td>375</td><td>72</td></tr><tr><td>Positive correlations</td><td>514 (52.7%)</td><td>232 (62.9%)</td><td>48 (60.5%)</td><td>493 (56.9%)</td><td>222 (59.2%)</td><td>44 (61.1%)</td></tr><tr><td>Negative correlations</td><td>461 (47.3%)</td><td>137 (37.1%)</td><td>28 (36.8%)</td><td>374 (43.1%)</td><td>153 (40.8%)</td><td>28 (38.9%)</td></tr><tr><td>Modularity</td><td>6.43</td><td>2.09</td><td>2.01</td><td>3.24</td><td>2.61</td><td>2.46</td></tr><tr><td>Network diameter</td><td>12</td><td>15</td><td>10</td><td>12</td><td>14</td><td>14</td></tr><tr><td>Average degree</td><td>8.13</td><td>3.43</td><td>1.69</td><td>7.5</td><td>3.85</td><td>1.58</td></tr><tr><td>Weighted degree</td><td>0.40</td><td>0.71</td><td>0.36</td><td>0.86</td><td>0.57</td><td>0.27</td></tr><tr><td>Clustering coefficient</td><td>0.52</td><td>0.39</td><td>0.29</td><td>0.46</td><td>0.38</td><td>0.21</td></tr><tr><td>Number of motifs</td><td>13,322</td><td>519</td><td>18</td><td>9549</td><td>561</td><td>12</td></tr></tbody></table>

<table><caption>TABLE 2 Keystone taxa of the bacterial communities within *R. microplus* by condition and in common.</caption><thead><tr><th>Condition</th><th>Keystone taxa by condition and in common</th></tr></thead><tbody><tr><td>July 2020</td><td>TRA3-20<br>Flavonifractor<br>Rhizobiaceae uncultured<br>Guggenheimella<br>Hydrocarboniphaga<br>Ga0077536<br>Ideonella<br>Nakamurella<br>Acidibacter<br>Staphylococcaceae<br>Parvibacter<br>Aquicella</td></tr><tr><td>September 2020</td><td>Ga0077536<br>Murdochielia<br>Acetobacteriaceae<br>Solimonas</td></tr><tr><td>March 2021</td><td>Quadrisphaera<br>Acidibacter<br>Anaplasma<br>Cellulomonadaceae</td></tr><tr><td>July 2020–September 2020</td><td>Ga0077536</td></tr><tr><td>July 2020–March 2021</td><td>Acidibacter</td></tr></tbody></table>

of all taxa in S-20. This node indirectly linked with three keystone taxa (Ga0077536, *Murdochielia*, *Solimonas*) (Table 2), through a direct association with *Haliangium* (Figure 2j). Notably, in the M-21 network, *Anaplasma* appeared as a keystone taxon (Table 2), with positive association module with seven nodes (Figure 2k). In the M-21 network, *Rickettsiella* was no longer part of the network, possibly due to co-exclusion caused by a negative interaction observed with *Anaplasma* in S20 (Figure 2i). Most direct associations of *Anaplasma* at a given time point were either displaced or did not co-occur with the taxon in the rest of the conditions, except for *Deferribacteraceae* and *JS1* (Caldribacteriota) in J-20 and M-21 (Table S7).

The evaluation of within-module ($Z_i$) and among-module ($P_i$) connectivity revealed similar distributions of nodes across the three conditions, with most nodes as peripheral (low $Z_i$ and $P_i$) and no network hubs identified (high $Z_i$ and $P_i$; Figure 2l-n). However, the presence of module hubs differed among the three groups, with nine, two and none taxa considered as module hubs in M-21, S-20, and J-20 networks, respectively (Figure 2l-n). The analysis revealed 24 and 14 unique connectors in S-20 and M-21 networks, respectively. Propionibacteriaceae and Lachnospiraceae UCG-002 promote coherence of both networks (Figure 2o, Table S8). *Anaplasma*'s position within the networks varied; it was peripheral in J-20 and S-20 networks (Figure 2l,m), but contributed to module coherence in the M-21 network (Figure 2n).

*Anaplasma* exhibited varying associations, evolving from diverse connections in J-20 to higher centrality and positive associations in S-20. Remarkably, in M-21, *Anaplasma* emerged as a keystone taxon, forming an independent cluster nested in the network establishing positives associations with other taxa. The absence of certain taxa and the potential co-exclusion of *Rickettsiella* may underscore the distinctive impact of *Anaplasma* on the community's dynamics over time.

### 3.3 | Influence of *Anaplasma* on the assembly, and hierarchy of *R. microplus* microbiome over time

To investigate *Anaplasma*'s impact on community assembly, network topology after its removal (woA) was analyzed and compared to the network with *Anaplasma* (wA). The removal of *Anaplasma* led to the loss of nodes in both the J-20 (woA) and S-20 (woA) networks, resulting in the depletion of both positive and negative correlations (Table 1, Figure S1A,B). Notably, within S-20 (woA) network, the associations changed with 16 new negative correlations appearing and 10 positive correlations being lost, suggesting that *Anaplasma* exerts an influence on interactions within the microbial community (Table 1, Figure S1B). In M-21 (woA), only one new connection formed within the network, and four positive correlations were lost (Table 1, Figure S1C). The removal of *Anaplasma* resulted in a decrease in the total number of motifs across all three networks, with the modularity
value in the J-20 (woA) network halving (Table 1). However, despite the removal of *Anaplasma*, the J-20 (woA) network retained a substantial number of motifs compared to the S-20 (woA) and M-21 (woA) networks, suggesting that association patterns among taxa with shared functions evolve over time (Table 1, Table S5). Overall, the removal of *Anaplasma* did not lead to significant changes in other topological features (Table 1). Regarding the composition, 53 nodes were shared by the three networks (Figure S1D, Table S9). J-20 (woA) network had the most unique nodes, followed by S-20 (woA) and M-21 (woA) (Figure S1D, Table S9).

Jaccard index comparison was used as a test to evaluate the local centrality measures in networks. The analysis considered dynamic interaction patterns within (wA and woA) networks at each time point and evaluated the impact of *Anaplasma*'s removal (wA vs. woA) at the same time point. Jaccard index values of the comparisons of central-ity measures, were lower than expected by random in the wA and woA network comparison for each time point ($p$ (≤ Jacc) <.05). With the exception, in both cases, of betweenness centrality in the J/S-20 comparison which has a random distribution (Tables S10 and S11). On the other hand, when the Jaccard index values were compared, the centrality measures between the (wA vs. woA) networks, in the same time point, were higher than expected by random (Table S12). The comparison of node clustering for wA and woA networks, in each time point, showed the higher ARI value for J-20 compared to S-20 networks, followed by S-20/M-21. The J-20/M-21 comparisons showed low similarities in clustering (Table 3). When both, the dynamic interaction patterns and impact of *Anaplasma*'s removal were analyzed, the high dissimilarity between J-20 and M-21 is in correspondence with what was observed visually and by the topological features values. When comparing wA vs. woA networks at the same time point, we found a higher ARI value, close to 1, indicating strong similarities and suggesting that *Anaplasma*'s removal did not affect clustering (Table 3).

The impact of removing *Anaplasma* was evaluated on two modules: M1, where *Anaplasma* was initially present, and M2, a module with a higher number of taxa and modularity value indirectly affected by *Anaplasma*'s removal. Micropepsaceae and *Craurococcus-Caldovatus* were the only common taxa found in the comparison between the J-20 (wA vs. woA) networks in modules 1 (J-20 M1) and 2 (J-20 M2), respectively (Figure 3a,b). After *Anaplasma*'s removal, both the composition

of the modules and the direct associations of the shared taxa changed drastically (Figure 3b,c, Table S13). In the S-20 network, *Anaplasma* is contained in the M1 in indirect association through *Haliagium* with keystone taxa found in M2 (Figure 3d). After *Anaplasma*'s removal, the composition of the M1 changes drastically, no taxon is shared (Figure 3e,f). Keystone taxa's presence in M2 (wA vs. woA) contributes to the similarity in the composition of the module, sharing a total of 23 taxa (Figure 3d-f, Table S13). As a result of the large number of shared nodes in M2 (wA vs. woA), due to the stability in the composition conferred by the presence of keystone taxa, the majority of the connected nodes represented in the differential sub-network correspond to M1 (wA and woA) (Figure 3f). After the analysis of both modules in the M-21 network, no taxon was shared, which may be due to the independent clustering formed by *Anaplasma* and its role as keystone taxa (Figure 3g-i, Table S13).

## 3.4 | Influence of *Anaplasma* on network robustness

A comprehensive analysis of network stability and capacity to withstand disturbances, such as node removal and addition, was performed over different months, both in the presence and absence of *Anaplasma* at the same time point. The robustness of networks with *Anaplasma* against directed and random attacks was assessed by comparing the fraction of nodes necessary for a connectivity loss of 40% and 80% over time.

Results indicated that J-20 (wA) and S-20 (wA) network connectivity remained similar after directed attacks in betweenness (Figure 4a), cascading (Figure S1E), degree (Figure S1F), as well as random attack (Figure S1G), removing equal fractions of nodes caused losses of connectivity of 40% and 80% (Table S14). This behavior corresponded with the random distribution for betweenness centrality observed in the comparison of local centrality measures (Table S10). In contrast, the M-21 network was severely affected during directed and random attacks due to its lower number of nodes and motifs compared to J-20 and S-20 (Figure 4a, Figure S1E-G, Table S14).

As demonstrated, *Anaplasma* influenced community assembly by slightly modifying topological features and drastically

<table><caption>TABLE 3 Network clustering comparisons.</caption><thead><tr><th>Conditions</th><th>Network comparisons</th><th>Adjusted Rand index (ARI)</th><th>p-Value</th></tr></thead><tbody><tr><td rowspan="3">With <i>Anaplasma</i> (wA)</td><td>J-20 (wA) vs. S-20 (wA)</td><td>0.25</td><td>0***</td></tr><tr><td>S-20 (wA) vs. M-21 (wA)</td><td>0.14</td><td>0***</td></tr><tr><td>J-20 (wA) vs. M-21 (wA)</td><td>0.09</td><td>0***</td></tr><tr><td rowspan="3">Without <i>Anaplasma</i> (woA)</td><td>J-20 (woA) vs. S-20 (woA)</td><td>0.26</td><td>0***</td></tr><tr><td>S-20 (woA) vs. M-21 (woA)</td><td>0.18</td><td>0***</td></tr><tr><td>J-20 (woA) vs. M-21 (woA)</td><td>0.09</td><td>0***</td></tr><tr><td rowspan="3">wA vs. woA</td><td>J-20 (wA vs. woA)</td><td>0.82</td><td>0***</td></tr><tr><td>S-20 (wA vs. woA)</td><td>0.81</td><td>0***</td></tr><tr><td>M-21 (wA vs. woA)</td><td>0.79</td><td>0***</td></tr></tbody></table>
value in the J-20 (woA) network halving (Table 1). However, despite the removal of *Anaplasma*, the J-20 (woA) network retained a substantial number of motifs compared to the S-20 (woA) and M-21 (woA) networks, suggesting that association patterns among taxa with shared functions evolve over time (Table 1, Table S5). Overall, the removal of *Anaplasma* did not lead to significant changes in other topological features (Table 1). Regarding the composition, 53 nodes were shared by the three networks (Figure S1D, Table S9). J-20 (woA) network had the most unique nodes, followed by S-20 (woA) and M-21 (woA) (Figure S1D, Table S9).

Jaccard index comparison was used as a test to evaluate the local centrality measures in networks. The analysis considered dynamic interaction patterns within (wA and woA) networks at each time point and evaluated the impact of *Anaplasma*'s removal (wA vs. woA) at the same time point. Jaccard index values of the comparisons of centrality measures, were lower than expected by random in the wA and woA network comparison for each time point ($p$ (≤ Jacc) <.05). With the exception, in both cases, of betweenness centrality in the J/S-20 comparison which has a random distribution (Tables S10 and S11). On the other hand, when the Jaccard index values were compared, the centrality measures between the (wA vs. woA) networks, in the same time point, were higher than expected by random (Table S12). The comparison of node clustering for wA and woA networks, in each time point, showed the higher ARI value for J-20 compared to S-20 networks, followed by S-20/M-21. The J-20/M-21 comparisons showed low similarities in clustering (Table 3). When both, the dynamic interaction patterns and impact of *Anaplasma*'s removal were analyzed, the high dissimilarity between J-20 and M-21 is in correspondence with what was observed visually and by the topological features values. When comparing wA vs. woA networks at the same time point, we found a higher ARI value, close to 1, indicating strong similarities and suggesting that *Anaplasma*'s removal did not affect clustering (Table 3).

The impact of removing *Anaplasma* was evaluated on two modules: M1, where *Anaplasma* was initially present, and M2, a module with a higher number of taxa and modularity value indirectly affected by *Anaplasma*'s removal. Micropepsaceae and *Craurococcus-Caldovatus* were the only common taxa found in the comparison between the J-20 (wA vs. woA) networks in modules 1 (J-20 M1) and 2 (J-20 M2), respectively (Figure 3a,b). After *Anaplasma*'s removal, both the composition

of the modules and the direct associations of the shared taxa changed drastically (Figure 3b,c, Table S13). In the S-20 network, *Anaplasma* is contained in the M1 in indirect association through *Haliagium* with keystone taxa found in M2 (Figure 3d). After *Anaplasma*'s removal, the composition of the M1 changes drastically, no taxon is shared (Figure 3e,f). Keystone taxa's presence in M2 (wA vs. woA) contributes to the similarity in the composition of the module, sharing a total of 23 taxa (Figure 3d-f, Table S13). As a result of the large number of shared nodes in M2 (wA vs. woA), due to the stability in the composition conferred by the presence of keystone taxa, the majority of the connected nodes represented in the differential sub-network correspond to M1 (wA and woA) (Figure 3f). After the analysis of both modules in the M-21 network, no taxon was shared, which may be due to the independent clustering formed by *Anaplasma* and its role as keystone taxa (Figure 3g-i, Table S13).

## 3.4 | Influence of *Anaplasma* on network robustness

A comprehensive analysis of network stability and capacity to withstand disturbances, such as node removal and addition, was performed over different months, both in the presence and absence of *Anaplasma* at the same time point. The robustness of networks with *Anaplasma* against directed and random attacks was assessed by comparing the fraction of nodes necessary for a connectivity loss of 40% and 80% over time.

Results indicated that J-20 (wA) and S-20 (wA) network connectivity remained similar after directed attacks in betweenness (Figure 4a), cascading (Figure S1E), degree (Figure S1F), as well as random attack (Figure S1G), removing equal fractions of nodes caused losses of connectivity of 40% and 80% (Table S14). This behavior corresponded with the random distribution for betweenness centrality observed in the comparison of local centrality measures (Table S10). In contrast, the M-21 network was severely affected during directed and random attacks due to its lower number of nodes and motifs compared to J-20 and S-20 (Figure 4a, Figure S1E-G, Table S14).

As demonstrated, *Anaplasma* influenced community assembly by slightly modifying topological features and drastically

<table><caption>TABLE 3 Network clustering comparisons.</caption><thead><tr><th>Conditions</th><th>Network comparisons</th><th>Adjusted Rand index (ARI)</th><th>p-Value</th></tr></thead><tbody><tr><td rowspan="3">With <i>Anaplasma</i> (wA)</td><td>J-20 (wA) vs. S-20 (wA)</td><td>0.25</td><td>0***</td></tr><tr><td>S-20 (wA) vs. M-21 (wA)</td><td>0.14</td><td>0***</td></tr><tr><td>J-20 (wA) vs. M-21 (wA)</td><td>0.09</td><td>0***</td></tr><tr><td rowspan="3">Without <i>Anaplasma</i> (woA)</td><td>J-20 (woA) vs. S-20 (woA)</td><td>0.26</td><td>0***</td></tr><tr><td>S-20 (woA) vs. M-21 (woA)</td><td>0.18</td><td>0***</td></tr><tr><td>J-20 (woA) vs. M-21 (woA)</td><td>0.09</td><td>0***</td></tr><tr><td rowspan="3">wA vs. woA</td><td>J-20 (wA vs. woA)</td><td>0.82</td><td>0***</td></tr><tr><td>S-20 (wA vs. woA)</td><td>0.81</td><td>0***</td></tr><tr><td>M-21 (wA vs. woA)</td><td>0.79</td><td>0***</td></tr></tbody></table>
J-20

(a)

Sub-networks
wA

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

(b)

Sub-networks
woA

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

(c)

J-20 M1

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

J-20 M2

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

S-20

(d)

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

(e)

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

(f)

S-20 M1

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

S-20 M2

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

M-21

(g)

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

(h)

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

(i)

M-21 M1

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

M-21 M2

<figure><img src="image_15.png" /><figcaption>FIGURE 3 Co-occurrence networks of the main modules (M1 and M2). Sub-networks that contain M1 and M2 from global co-occurrence networks in *Anaplasma*'s presence (wA) and removal (woA). Venn diagram displaying the comparison of module composition (M1-M2) in *Anaplasma*'s presence and removal: (a) J-20 wA, (b) J-20 woA, (c) J-20 (M1M2), (d) S-20 wA, (e) S-20 woA, (f) S-20 (M1M2), (g) M-21 wA, (h) M-21 woA, and (i) M-21 (M1M2).</figcaption></figure>

changing the composition, especially in modules where it was present. *Anaplasma*'s removal did not alter the robustness of J-20 woA and S-20 woA networks against directed (cascading and degree) and random attacks (Figure S1H–M, Table S14). However, the M-21 woA network experienced a rapid loss of connectivity (80%) compared to M-21 wA at a lower fraction of removed nodes, against both directed and random attacks (Figure S1N–P, Table S14). The most notable robustness behavior occurred when comparing networks in the *Anaplasma*'s presence and removal after a directed attack in

betweenness (Figure 4b–d). In J-20 and S-20, *Anaplasma*'s removal did not significantly change the proportion of nodes required for a loss of connectivity (Figure 4b,c, Table S14). Conversely, in M-21, *Anaplasma*'s removal conferred instability to the networks, making them more susceptible to disturbances (Figure 4b,c, Table S14).

As a result of the addition of nodes, the values of LCC Size and average path length were compared between the networks at each time point (J-20, S-20 and M-21) and at the same time point (wA vs. woA) (Figure 4e–I, Table S15). In general, J-20 and S-20 networks
J-20

(a)

Sub-networks
wA

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

(b)

Sub-networks
woA

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

(c)

J-20 M1

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

J-20 M2

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

S-20

(d)

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

(e)

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

(f)

S-20 M1

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

S-20 M2

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

M-21

(g)

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

(h)

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

(i)

M-21 M1

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

M-21 M2

<figure><img src="image_31.png" /><figcaption>FIGURE 3 Co-occurrence networks of the main modules (M1 and M2). Sub-networks that contain M1 and M2 from global co-occurrence networks in *Anaplasma*'s presence (wA) and removal (woA). Venn diagram displaying the comparison of module composition (M1-M2) in *Anaplasma*'s presence and removal: (a) J-20 wA, (b) J-20 woA, (c) J-20 (M1M2), (d) S-20 wA, (e) S-20 woA, (f) S-20 (M1M2), (g) M-21 wA, (h) M-21 woA, and (i) M-21 (M1M2).</figcaption></figure>

changing the composition, especially in modules where it was present. *Anaplasma*'s removal did not alter the robustness of J-20 woA and S-20 woA networks against directed (cascading and degree) and random attacks (Figure S1H–M, Table S14). However, the M-21 woA network experienced a rapid loss of connectivity (80%) compared to M-21 wA at a lower fraction of removed nodes, against both directed and random attacks (Figure S1N–P, Table S14). The most notable robustness behavior occurred when comparing networks in the *Anaplasma*'s presence and removal after a directed attack in

betweenness (Figure 4b–d). In J-20 and S-20, *Anaplasma*'s removal did not significantly change the proportion of nodes required for a loss of connectivity (Figure 4b,c, Table S14). Conversely, in M-21, *Anaplasma*'s removal conferred instability to the networks, making them more susceptible to disturbances (Figure 4b,c, Table S14).

As a result of the addition of nodes, the values of LCC Size and average path length were compared between the networks at each time point (J-20, S-20 and M-21) and at the same time point (wA vs. woA) (Figure 4e–I, Table S15). In general, J-20 and S-20 networks
<figure><img src="image_16.png" /><figcaption>FIGURE 4 Robustness comparison between J-20, S-20, and M-21 networks after removal and addition of nodes in *Anaplasma*'s presence and removal. (a) Connectivity loss measured after directed attack in *Anaplasma*'s presence between networks for each time point, removing first the nodes with the highest betweenness centrality (J-20 wA/S-20 wA/M-21 wA). Connectivity loss measured after directed attack in *Anaplasma*'s presence and removal between networks in the same point time, removing first the nodes with the highest betweenness centrality: (b) J-20 (wA vs. woA), (c) S-20 (wA vs. woA), (d) M-21 (wA vs. woA). Largest connected component (LCC) values are represented and compared in *Anaplasma*'s presence between networks for each time point: (e) J-20 (wA)/S-20 (wA)/M-21 (wA). LCC values are represented and compared in *Anaplasma*'s presence and removal between networks in the same point time: (f) J-20 (wA vs. woA), (g) S-20 (wA vs. woA), and (h) M-21 (wA vs. woA). Average path length (APL) values are represented and compared in *Anaplasma*'s presence between networks for each time point: (i) J-20 (wA)/S-20 (wA)/M-21 (wA). APL values are represented and compared in *Anaplasma*'s presence and removal between networks in the same point time: (j) J-20 (wA vs. woA), (k) S-20 (wA vs. woA), and (l) M-21 (wA vs. woA).</figcaption></figure>

maintained similar and larger LCC size values compared to M-21 network (Figure 4e, Table S15). The wA networks maintained larger values of LCC size compared to the woA networks in J-20 (Figure 4f, Table S15) and S-20 (Figure 4g, Table S15). In contrast, the M-21

(wA vs. woA) networks the LCC size values were more similar with overlap at some points (Figure 4h, Table S15). However, in all cases, the value of LCC Size was increased as greater number of nodes were added. With respect to the average path length J-20 and S-20
<figure><img src="image_32.png" /><figcaption>**FIGURE 4** Robustness comparison between J-20, S-20, and M-21 networks after removal and addition of nodes in *Anaplasma*'s presence and removal. (a) Connectivity loss measured after directed attack in *Anaplasma*'s presence between networks for each time point, removing first the nodes with the highest betweenness centrality (J-20 wA/S-20 wA/M-21 wA). Connectivity loss measured after directed attack in *Anaplasma*'s presence and removal between networks in the same point time, removing first the nodes with the highest betweenness centrality: (b) J-20 (wA vs. woA), (c) S-20 (wA vs. woA), (d) M-21 (wA vs. woA). Largest connected component (LCC) values are represented and compared in *Anaplasma*'s presence between networks for each time point: (e) J-20 (wA)/S-20 (wA)/M-21 (wA). LCC values are represented and compared in *Anaplasma*'s presence and removal between networks in the same point time: (f) J-20 (wA vs. woA), (g) S-20 (wA vs. woA), and (h) M-21 (wA vs. woA). Average path length (APL) values are represented and compared in *Anaplasma*'s presence between networks for each time point: (i) J-20 (wA)/S-20 (wA)/M-21 (wA). APL values are represented and compared in *Anaplasma*'s presence and removal between networks in the same point time: (j) J-20 (wA vs. woA), (k) S-20 (wA vs. woA), and (l) M-21 (wA vs. woA).</figcaption></figure>

maintained similar and larger LCC size values compared to M-21 network (Figure 4e, Table S15). The wA networks maintained larger values of LCC size compared to the woA networks in J-20 (Figure 4f, Table S15) and S-20 (Figure 4g, Table S15). In contrast, the M-21

(wA vs. woA) networks the LCC size values were more similar with overlap at some points (Figure 4h, Table S15). However, in all cases, the value of LCC Size was increased as greater number of nodes were added. With respect to the average path length J-20 and S-20
networks maintained similar values and behavior compared to M-21 (Figure 4i, Table S15). Examining the (wA vs. woA) APL comparison, the similarity and overlapping of the values at some points between the networks in the presence and absence of *Anaplasma* was evident in J-20 and M-21, while in S-20 they were different and approached as the number of nodes added increased (Figure 4j–l, Table S15). This predictive behavior, in the conditions where *Anaplasma* is important in the community either due to abundance, direct association with key taxon or its own role as a keystone, suggests that its presence confers stability in the assembly of the microbial community in the face of a possible loss or addition of taxa.

4 | DISCUSSION

The dynamic nature of microbial communities within ticks is a topic of considerable interest due to its potential implications for tick-borne diseases. In this study, we aimed to elucidate the stable and variable elements in the *R. microplus* microbiome over time and explore the interactions between microbial communities and specific pathogens, with a focus on *Anaplasma* nesting within these communities over time. The tick microbiome is considered a fluctuating microecosystem influenced by internal factors, such as interactions among pathogenic and nonpathogenic microorganisms, and responses to external perturbations (Aguilar-Díaz et al., 2021; Cabezas-Cruz et al., 2018; Chicana et al., 2019; Swei & Kwan, 2017; Wu-Chuang et al., 2021). To analyze the interplay between pathogens and the temporal dynamics of microbial communities: the covariates (engorgement level, tick species, stages and host) were reduced, individual tick samples positive for single *A. marginale* infection were selected, and extraction process controls were used.

Observations revealed differences in α-diversity between ticks collected in different years and a significantly higher relative abundance of *Anaplasma* in ticks from M-21 compared to J-20 and S-20, indicating substantial temporal variation in the *R. microplus* microbiome. Changes in microbial diversity demonstrated that the tick microbiome is a dynamic system with varying patterns of species abundance, suggesting interference in the acquisition and colonization process in *R. microplus* due to *Anaplasma*'s relative abundance variability. High microbial diversity within ecosystems makes communities more resistant to pathogen colonization, as trophic interactions prevent dominance by a single pathogen (Wei et al., 2015). Although differences in the abundance of bacterial taxa in the presence of *Anaplasma* imply its impact on the *R. microplus* microbiome, further exploration is necessary to understand the co-occurrence patterns it establishes with the bacterial community over time.

The microbial community assembly in *R. microplus* reveals a dynamic interplay among microorganisms, forming intricate and ever-changing microbial consortia. The configuration of microbial communities in ecosystems is likely influenced by various biotic interactions, such as commensalism, mutualism, and parasitism, among the composing microorganisms (Freilich et al., 2011; Nemergut et al., 2013; Stolyar et al., 2007). Comparative analyses of bacterial

assembly in hard ticks (Maitre et al., 2023) and soft ticks (Piloto-Sardiñas, Cano-Argüelles, et al., 2023) indicate a higher frequency of positive interactions between bacterial taxa. Our co-occurrence pattern analysis in *R. microplus* underscores a predominant occurrence of cooperative interactions across all three conditions. However, a notable decline in associations among community members was evident over time.

The nature of microbial cooperation or competition hinges on factors such as metabolic diversity (Stolyar et al., 2007), genotypic and phenotypic variations between species (Ackermann, 2015; Mitri & Richard Foster, 2013), environmental carrying capacity (Freilich et al., 2011) and the presence, role and microbial load of pathogenic microorganisms (Adegoke et al., 2020; Maitre et al., 2022). Consequently, the observed reduction in associations across the three conditions may be attributed to the loss of bacterial taxa or the increasing influence of *Anaplasma* within the community over time.

Over time, assembly patterns in the microbial community have evolved, marked by the emergence of *Anaplasma* as a key taxon. *Anaplasma* forms an independent cluster within the M-21 network, while there is a notable reduction in bacterial taxa, with *Rickettsiella* being one of the taxa displaced from the community. Previous investigations in *I. ricinus* ticks indicated no preferential or antagonistic association between *Rickettsiella* and *Borrelia* species (*B. burgdorferi* and *B. miyamotoi*) (Garcia-Vozmediano et al., 2022). In contrast, our network analysis suggests a potential displacement of *Rickettsiella* by *Anaplasma*, possibly due to competition or co-exclusion interactions established between the two in the preceding condition. The variations in co-occurrence and nesting patterns of *Anaplasma* at the three time points, along with its presumed impact on the reduction of bacterial taxa, indicate a critical role played by this pathogen in shaping interactions within the microbial community and influencing the dynamics of the *R. microplus* microbiome.

*In silico* removal of *Anaplasma* reveals modified interaction patterns, motif numbers, module compositions, and stability against disturbances across time points. The most significant impacts are observed when *Anaplasma* establishes direct interaction with keystone taxa or is considered one. Recent evidence indicates that bacterial pathogens like *R. helvetica* (Maitre et al., 2022), and *B. afzelii* (Wu-Chuang et al., 2023) modulates the composition and assembly of the bacterial community in *I. ricinus*, while *A. phagocytophilum* modifies the *I. scapularis* microbiome thus facilitating its colonization (Abraham et al., 2017).

Studies also show that the removal of *Rickettsia* affects conserved patterns of community assembly in *H. marginatum* and *R. bursa* ticks, suggesting it acts as a community assembly driver (Maitre et al., 2023). Conversely, in *R. microplus*, infection by the protozoan *Theileria* sp. leads to a significant reduction in the bacterial community, termed “pathogen-induced dysbiosis” by the authors (Adegoke et al., 2020). Keystone taxa, in general, influence the composition, structure, assembly, and functioning of microbial communities (Banerjee et al., 2018; Modlmeier et al., 2014). Those sustaining and stabilizing a microbiome associated with pathological states are referred to as “keystone pathogens” (Hajishengallis
networks maintained similar values and behavior compared to M-21 (Figure 4i, Table S15). Examining the (wA vs. woA) APL comparison, the similarity and overlapping of the values at some points between the networks in the presence and absence of *Anaplasma* was evident in J-20 and M-21, while in S-20 they were different and approached as the number of nodes added increased (Figure 4j–l, Table S15). This predictive behavior, in the conditions where *Anaplasma* is important in the community either due to abundance, direct association with key taxon or its own role as a keystone, suggests that its presence confers stability in the assembly of the microbial community in the face of a possible loss or addition of taxa.

4 | DISCUSSION

The dynamic nature of microbial communities within ticks is a topic of considerable interest due to its potential implications for tick-borne diseases. In this study, we aimed to elucidate the stable and variable elements in the *R. microplus* microbiome over time and explore the interactions between microbial communities and specific pathogens, with a focus on *Anaplasma* nesting within these communities over time. The tick microbiome is considered a fluctuating microecosystem influenced by internal factors, such as interactions among pathogenic and nonpathogenic microorganisms, and responses to external perturbations (Aguilar-Díaz et al., 2021; Cabezas-Cruz et al., 2018; Chicana et al., 2019; Swei & Kwan, 2017; Wu-Chuang et al., 2021). To analyze the interplay between pathogens and the temporal dynamics of microbial communities: the covariates (engorgement level, tick species, stages and host) were reduced, individual tick samples positive for single *A. marginale* infection were selected, and extraction process controls were used.

Observations revealed differences in α-diversity between ticks collected in different years and a significantly higher relative abundance of *Anaplasma* in ticks from M-21 compared to J-20 and S-20, indicating substantial temporal variation in the *R. microplus* microbiome. Changes in microbial diversity demonstrated that the tick microbiome is a dynamic system with varying patterns of species abundance, suggesting interference in the acquisition and colonization process in *R. microplus* due to *Anaplasma*'s relative abundance variability. High microbial diversity within ecosystems makes communities more resistant to pathogen colonization, as trophic interactions prevent dominance by a single pathogen (Wei et al., 2015). Although differences in the abundance of bacterial taxa in the presence of *Anaplasma* imply its impact on the *R. microplus* microbiome, further exploration is necessary to understand the co-occurrence patterns it establishes with the bacterial community over time.

The microbial community assembly in *R. microplus* reveals a dynamic interplay among microorganisms, forming intricate and ever-changing microbial consortia. The configuration of microbial communities in ecosystems is likely influenced by various biotic interactions, such as commensalism, mutualism, and parasitism, among the composing microorganisms (Freilich et al., 2011; Nemergut et al., 2013; Stolyar et al., 2007). Comparative analyses of bacterial

assembly in hard ticks (Maitre et al., 2023) and soft ticks (Piloto-Sardiñas, Cano-Argüelles, et al., 2023) indicate a higher frequency of positive interactions between bacterial taxa. Our co-occurrence pattern analysis in *R. microplus* underscores a predominant occurrence of cooperative interactions across all three conditions. However, a notable decline in associations among community members was evident over time.

The nature of microbial cooperation or competition hinges on factors such as metabolic diversity (Stolyar et al., 2007), genotypic and phenotypic variations between species (Ackermann, 2015; Mitri & Richard Foster, 2013), environmental carrying capacity (Freilich et al., 2011) and the presence, role and microbial load of pathogenic microorganisms (Adegoke et al., 2020; Maitre et al., 2022). Consequently, the observed reduction in associations across the three conditions may be attributed to the loss of bacterial taxa or the increasing influence of *Anaplasma* within the community over time.

Over time, assembly patterns in the microbial community have evolved, marked by the emergence of *Anaplasma* as a key taxon. *Anaplasma* forms an independent cluster within the M-21 network, while there is a notable reduction in bacterial taxa, with *Rickettsiella* being one of the taxa displaced from the community. Previous investigations in *I. ricinus* ticks indicated no preferential or antagonistic association between *Rickettsiella* and *Borrelia* species (*B. burgdorferi* and *B. miyamotoi*) (Garcia-Vozmediano et al., 2022). In contrast, our network analysis suggests a potential displacement of *Rickettsiella* by *Anaplasma*, possibly due to competition or co-exclusion interactions established between the two in the preceding condition. The variations in co-occurrence and nesting patterns of *Anaplasma* at the three time points, along with its presumed impact on the reduction of bacterial taxa, indicate a critical role played by this pathogen in shaping interactions within the microbial community and influencing the dynamics of the *R. microplus* microbiome.

*In silico* removal of *Anaplasma* reveals modified interaction patterns, motif numbers, module compositions, and stability against disturbances across time points. The most significant impacts are observed when *Anaplasma* establishes direct interaction with keystone taxa or is considered one. Recent evidence indicates that bacterial pathogens like *R. helvetica* (Maitre et al., 2022), and *B. afzelii* (Wu-Chuang et al., 2023) modulates the composition and assembly of the bacterial community in *I. ricinus*, while *A. phagocytophilum* modifies the *I. scapularis* microbiome thus facilitating its colonization (Abraham et al., 2017).

Studies also show that the removal of *Rickettsia* affects conserved patterns of community assembly in *H. marginatum* and *R. bursa* ticks, suggesting it acts as a community assembly driver (Maitre et al., 2023). Conversely, in *R. microplus*, infection by the protozoan *Theileria* sp. leads to a significant reduction in the bacterial community, termed “pathogen-induced dysbiosis” by the authors (Adegoke et al., 2020). Keystone taxa, in general, influence the composition, structure, assembly, and functioning of microbial communities (Banerjee et al., 2018; Modlmeier et al., 2014). Those sustaining and stabilizing a microbiome associated with pathological states are referred to as “keystone pathogens” (Hajishengallis
et al., 2012). When these pathogens dominate the community by suppressing other microbes, they induce an alteration in the microbiome (Hajishengallis et al., 2012).

The changes in microbial community assembly induced by *Anaplasma* and its role as a keystone taxon suggest a “keystone pathogen-induced dysbiosis” in *R. microplus*. Although further studies are required to validate this hypothesis, we propose that the effects induced by *Anaplasma* on the community dynamics likely create a favorable environment facilitating its colonization within *R. microplus*.

In summary, our findings indicate that the microbiome of *R. microplus* is a dynamic system. The results highlight the significant influence of *A. marginale* on microbial communities, suggesting that its high relative abundance, widespread presence, and increasing importance shape the dynamics of microbial interactions. This, in turn, may potentially alter the physiology and vector capacity of *R. microplus*. Analyzing microbiome dynamics is crucial for understanding the cause-and-effect relationships in responses to both biotic and abiotic factors. Recognizing the temporal dimension in tick-microbiome interactions becomes fundamental for strategies aiming to manipulate bacterial communities to modify tick physiology and vector capacity. Nevertheless, further research is required to uncover the mechanistic basis of these effects and their broader implications for the transmission dynamics of tick-borne pathogens, the ecology of tick microbiome, and the development of effective control strategies.

## AUTHOR CONTRIBUTIONS

**Elianne Piloto-Sardiñas:** Conceptualization (equal); formal analysis (lead); investigation (lead); validation (lead); visualization (lead); writing – original draft (lead); writing – review and editing (equal). **Lianet Abuin-Denis:** Formal analysis (supporting); writing – review and editing (equal). **Apolline Maitre:** Formal analysis (supporting); supervision (supporting); writing – review and editing (equal). **Angélique Foucault-Simonin:** Investigation (supporting); writing – review and editing (equal). **Belkis Corona-González:** Conceptualization (equal); investigation (supporting); resources (equal); supervision (equal); writing – review and editing (equal). **Cristian Díaz-Corona:** Writing – review and editing (equal). **Lisset Roblejo-Arias:** Writing – review and editing (equal). **Lourdes Mateos-Hernández:** Formal analysis (lead); writing – review and editing (equal). **Roxana Marrero-Perera:** Investigation (equal); writing – review and editing (supporting). **Dasiel Obregon:** Software (lead); writing – review and editing (equal). **Karolína Svobodová:** Formal analysis (equal); methodology (equal); writing – review and editing (equal). **Alejandra Wu-Chuang:** Supervision (lead). **Alejandro Cabezas-Cruz:** Conceptualization (equal); resources (equal); supervision (equal); writing – original draft (supporting); writing – review and editing (lead).

## ACKNOWLEDGEMENTS

This work was funded by a grant to from Piloto-Sardiñas E. in the frame of FSPI project (Fonds de solidarité pour les projects innovants)

VectoCaribe from the embassy of France in Cuba. UMR BIPAR is supported by the French Government's Investissement d'Avenir program, Laboratoire d'Excellence 'Integrative Biology of Emerging Infectious Diseases' (grant no. ANR-10-LABX-62-IBEID). Apolline Maitre is supported by the 'Collectivité de Corse', grant: 'Formations superieures' (SGCE - RAPPORT No 0300). AI technologies were used to improve spelling, grammar, and general editing of the text. AI technologies were not used to produce scientific insights, or drawing scientific conclusions.

## CONFLICT OF INTEREST STATEMENT

The authors declare no competing interests.

## DATA AVAILABILITY STATEMENT

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

## BENEFIT-SHARING STATEMENT

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

## ORCID

Dasiel Obregon ID https://orcid.org/0000-0002-5786-1114

Alejandro Cabezas-Cruz ID https://orcid.

org/0000-0002-8660-730X

## REFERENCES

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

Ackermann, M. (2015). A functional perspective on phenotypic heterogeneity in microorganisms. *Nature Reviews Microbiology*, 13(8), 497-508.

Adegoke, A., Kumar, D., Bobo, C., Rashid, M. I., Durrani, A. Z., Sajid, M. S., & Karim, S. (2020). Tick-borne pathogens shape the native microbiome within tick vectors. *Microorganisms*, 8(9), 1299.

Aguilar-Díaz, H., Quiroz-Castañeda, R. E., Cobaxin-Cárdenas, M., Salinas-Estrella, E., & Amaro-Estrada, I. (2021). Advances in the study of the tick cattle microbiota and the influence on vectorial capacity. *Frontiers in Veterinary Science*, 8, 710352.

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

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

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

Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., & Al-Ghalith, G. A. (2019). Reproducible, interactive, scalable
et al., 2012). When these pathogens dominate the community by suppressing other microbes, they induce an alteration in the microbiome (Hajishengallis et al., 2012).

The changes in microbial community assembly induced by *Anaplasma* and its role as a keystone taxon suggest a “keystone pathogen-induced dysbiosis” in *R. microplus*. Although further studies are required to validate this hypothesis, we propose that the effects induced by *Anaplasma* on the community dynamics likely create a favorable environment facilitating its colonization within *R. microplus*.

In summary, our findings indicate that the microbiome of *R. microplus* is a dynamic system. The results highlight the significant influence of *A. marginale* on microbial communities, suggesting that its high relative abundance, widespread presence, and increasing importance shape the dynamics of microbial interactions. This, in turn, may potentially alter the physiology and vector capacity of *R. microplus*. Analyzing microbiome dynamics is crucial for understanding the cause-and-effect relationships in responses to both biotic and abiotic factors. Recognizing the temporal dimension in tick-microbiome interactions becomes fundamental for strategies aiming to manipulate bacterial communities to modify tick physiology and vector capacity. Nevertheless, further research is required to uncover the mechanistic basis of these effects and their broader implications for the transmission dynamics of tick-borne pathogens, the ecology of tick microbiome, and the development of effective control strategies.

## AUTHOR CONTRIBUTIONS

**Elianne Piloto-Sardiñas:** Conceptualization (equal); formal analysis (lead); investigation (lead); validation (lead); visualization (lead); writing – original draft (lead); writing – review and editing (equal). **Lianet Abuin-Denis:** Formal analysis (supporting); writing – review and editing (equal). **Apolline Maitre:** Formal analysis (supporting); supervision (supporting); writing – review and editing (equal). **Angélique Foucault-Simonin:** Investigation (supporting); writing – review and editing (equal). **Belkis Corona-González:** Conceptualization (equal); investigation (supporting); resources (equal); supervision (equal); writing – review and editing (equal). **Cristian Díaz-Corona:** Writing – review and editing (equal). **Lisset Roblejo-Arias:** Writing – review and editing (equal). **Lourdes Mateos-Hernández:** Formal analysis (lead); writing – review and editing (equal). **Roxana Marrero-Perera:** Investigation (equal); writing – review and editing (supporting). **Dasiel Obregon:** Software (lead); writing – review and editing (equal). **Karolína Svobodová:** Formal analysis (equal); methodology (equal); writing – review and editing (equal). **Alejandra Wu-Chuang:** Supervision (lead). **Alejandro Cabezas-Cruz:** Conceptualization (equal); resources (equal); supervision (equal); writing – original draft (supporting); writing – review and editing (lead).

## ACKNOWLEDGEMENTS

This work was funded by a grant to from Piloto-Sardiñas E. in the frame of FSPI project (Fonds de solidarité pour les projects innovants)

VectoCaribe from the embassy of France in Cuba. UMR BIPAR is supported by the French Government's Investissement d'Avenir program, Laboratoire d'Excellence 'Integrative Biology of Emerging Infectious Diseases' (grant no. ANR-10-LABX-62-IBEID). Apolline Maitre is supported by the 'Collectivité de Corse', grant: 'Formations superieures' (SGCE - RAPPORT No 0300). AI technologies were used to improve spelling, grammar, and general editing of the text. AI technologies were not used to produce scientific insights, or drawing scientific conclusions.

## CONFLICT OF INTEREST STATEMENT

The authors declare no competing interests.

## DATA AVAILABILITY STATEMENT

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

## BENEFIT-SHARING STATEMENT

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

## ORCID

Dasiel Obregon ID https://orcid.org/0000-0002-5786-1114

Alejandro Cabezas-Cruz ID https://orcid.

org/0000-0002-8660-730X

## REFERENCES

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

Ackermann, M. (2015). A functional perspective on phenotypic heterogeneity in microorganisms. *Nature Reviews Microbiology*, 13(8), 497-508.

Adegoke, A., Kumar, D., Bobo, C., Rashid, M. I., Durrani, A. Z., Sajid, M. S., & Karim, S. (2020). Tick-borne pathogens shape the native microbiome within tick vectors. *Microorganisms*, 8(9), 1299.

Aguilar-Díaz, H., Quiroz-Castañeda, R. E., Cobaxin-Cárdenas, M., Salinas-Estrella, E., & Amaro-Estrada, I. (2021). Advances in the study of the tick cattle microbiota and the influence on vectorial capacity. *Frontiers in Veterinary Science*, 8, 710352.

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

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

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

Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., & Al-Ghalith, G. A. (2019). Reproducible, interactive, scalable
and extensible microbiome data science using QIIME 2. *Nature Biotechnology*, 37(8), 852–857.

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

Cabezas-Cruz, A., Pollet, T., Estrada-Peña, A., Allain, E., Bonnet, S. I., & Moutailler, S. (2018). Handling the microbial complexity associated to ticks. In M. Abubakar & P. K. Perera (Eds.), *Ticks and tick-borne pathogens* (Vol. 1, pp. 1–36). IntechOpen Limited.

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

Cao, X., Zhao, D., Xu, H., Huang, R., Zeng, J., & Yu, Z. (2018). Heterogeneity of interactions of microbial communities in regions of Taihu Lake with different nutrient loadings: A network analysis. *Scientific Reports*, 8(1), 8890.

Chicana, B., Couper, L. I., Kwan, J. Y., Tahiraj, E., & Swei, A. (2019). Comparative microbiome profiles of sympatric tick species from the far-western United States. *Insects*, 10(10), 353.

Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A., & Callahan, B. J. (2018). Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. *Microbiome*, 6, 1–14.

De La Fourniere, S., Guillemi, E. C., Paletta, 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.

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

Estrada-Peña, A., Bouattour, A., Camicas, J. L., & Walker, A. R. (2004). *Ticks of domestic animals in the Mediterranean region* (p. 131). University of Zaragoza.

Faith, D. P. (1992). Conservation evaluation and phylogenetic diversity. *Biological Conservation*, 61, 1–10.

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

Freilich, S., Zarecki, R., Eilam, O., Segal, E. S., Henry, C. S., Kupiec, M., Gophna, U., Sharan, R., & Ruppin, E. (2011). Competitive and co-operative metabolic interactions in bacterial communities. *Nature Communications*, 2(1), 589.

Freitas, S., Yang, D., Kumar, S., Tong, H., & Chau, D. H. (2021). Evaluating graph vulnerability and robustness using tiger. In *proceedings of the 30th ACM international conference on Information & Knowledge Management* (pp. 4495–4503).

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

Garcia-Vozmediano, A., Tomassone, L., Fonville, M., Bertolotti, L., Heylen, D., Fabri, N. D., Medlock, J. M., Nijhof, A. M., Hansford, K. M., Sprong, H., & Krawczyk, A. I. (2022). The genetic diversity of *Rickettsiella* symbionts in *Ixodes ricinus* throughout Europe. *Microbial Ecology*, 84, 1–14.

Gondard, M., Delannoy, S., Pinarello, V., Aprelon, R., Devillers, E., Galon, C., Pradel, J., Vayssier-Taussat, M., Albina, E., & Moutailler, S. (2020). Upscaling the surveillance of tick-borne pathogens in the French Caribbean islands. *Pathogens*, 9(3), 176.

Grech-Angelini, S., Stachurski, F., Vayssier-Taussat, M., Devillers, E., Casabianca, F., Lancelot, R., Uilenberg, G., & Moutailler, S. (2020). Tick-borne pathogens in ticks (Acari: Ixodidae) collected from various domestic and wild hosts in Corsica (France), a Mediterranean Island environment. *Transboundary and Emerging Diseases*, 67(2), 745–757.

Guimera, R., & Nunes Amaral, L. A. (2005). Functional cartography of complex metabolic networks. *Nature*, 433(7028), 895–900.

Guo, B., Zhang, L., Sun, H., Gao, M., Yu, N., Zhang, Q., Mou, A., & Liu, Y. (2022). Microbial co-occurrence network topological properties link with reactor parameters and reveal importance of low-abundance genera. *npj Biofilms and Microbiomes*, 8(1), 3.

Hajishengallis, G., Darveau, R. P., & Curtis, M. A. (2012). The keystone-pathogen hypothesis. *Nature Reviews Microbiology*, 10(10), 717–725.

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

Krawczyk, A. I., Röttjers, L., Fonville, M., Takumi, K., Takken, W., Faust, K., & Sprong, H. (2022). Quantitative microbial population study reveals geographical differences in bacterial symbionts of *Ixodes ricinus*. *Microbiome*, 10(1), 120.

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

Lejal, E., Marsot, M., Chalvet-Monfray, K., Cosson, J. F., Moutailler, S., Vayssier-Taussat, M., & Pollet, T. (2019). A three-year assessment of *Ixodes ricinus*-borne pathogens in a French peri-urban forest. *Parasites & Vectors*, 12(1), 1–14.

Lhomme, S. (2015). Analyse spatiale de la structure des réseaux techniques dans un contexte de risques. *Cybergeo: European Journal of Geography* [En ligne], Systèmes, Modélisation, Géostatistiques. http://journals.openedition.org/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. *Scientific Reports*, 12(1), 11464.

Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Piloto-Sardiñas, E., Foucault-Simonin, A., Cicculli, V., Moutailler, S., Paoli, J. C., Falchi, A., Obregón, D., & Cabezas-Cruz, A. (2023). Rickettsial pathogens drive microbiota assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks. *Molecular Ecology*, 32(16), 4660–4676.

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

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

Mitri, S., & Richard Foster, K. (2013). The genotypic view of social interactions in microbial communities. *Annual Review of Genetics*, 47, 247–273.

Modlmeier, A. P., Keiser, C. N., Watters, J. V., Sih, A., & Pruitt, J. N. (2014). The keystone individual concept: An ecological and evolutionary overview. *Animal Behaviour*, 89, 53–62.

Narasimhan, S., Rajeevan, N., Liu, L., Zhao, Y. O., Heisig, J., Pan, J., Eppler-Epstein, R., DePonte, K., Fish, D., & Fikrig, E. (2014). Gut microbiota of the tick vector *Ixodes scapularis* modulate colonization of the Lyme disease spirochete. *Cell Host & Microbe*, 15(1), 58–71.

Narasimhan, S., Schuijt, T. J., Abraham, N. M., Rajeevan, N., Coumou, J., Graham, M., Robson, A., Wu, M. J., Daffre, S., Hovius, J. W., & Fikrig, E. (2017). Modulation of the tick gut milieu by a secreted tick protein favors borrelia burgdorferi colonization. *Nature Communications*, 8(1), 184.
and extensible microbiome data science using QIIME 2. *Nature Biotechnology*, 37(8), 852–857.

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

Cabezas-Cruz, A., Pollet, T., Estrada-Peña, A., Allain, E., Bonnet, S. I., & Moutailler, S. (2018). Handling the microbial complexity associated to ticks. In M. Abubakar & P. K. Perera (Eds.), *Ticks and tick-borne pathogens* (Vol. 1, pp. 1–36). IntechOpen Limited.

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

Cao, X., Zhao, D., Xu, H., Huang, R., Zeng, J., & Yu, Z. (2018). Heterogeneity of interactions of microbial communities in regions of Taihu Lake with different nutrient loadings: A network analysis. *Scientific Reports*, 8(1), 8890.

Chicana, B., Couper, L. I., Kwan, J. Y., Tahiraj, E., & Swei, A. (2019). Comparative microbiome profiles of sympatric tick species from the far-western United States. *Insects*, 10(10), 353.

Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A., & Callahan, B. J. (2018). Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. *Microbiome*, 6, 1–14.

De La Fourniere, S., Guillemi, E. C., Paletta, 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.

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

Estrada-Peña, A., Bouattour, A., Camicas, J. L., & Walker, A. R. (2004). *Ticks of domestic animals in the Mediterranean region* (p. 131). University of Zaragoza.

Faith, D. P. (1992). Conservation evaluation and phylogenetic diversity. *Biological Conservation*, 61, 1–10.

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

Freilich, S., Zarecki, R., Eilam, O., Segal, E. S., Henry, C. S., Kupiec, M., Gophna, U., Sharan, R., & Ruppin, E. (2011). Competitive and co-operative metabolic interactions in bacterial communities. *Nature Communications*, 2(1), 589.

Freitas, S., Yang, D., Kumar, S., Tong, H., & Chau, D. H. (2021). Evaluating graph vulnerability and robustness using tiger. In *proceedings of the 30th ACM international conference on Information & Knowledge Management* (pp. 4495–4503).

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

Garcia-Vozmediano, A., Tomassone, L., Fonville, M., Bertolotti, L., Heylen, D., Fabri, N. D., Medlock, J. M., Nijhof, A. M., Hansford, K. M., Sprong, H., & Krawczyk, A. I. (2022). The genetic diversity of *Rickettsiella* symbionts in *Ixodes ricinus* throughout Europe. *Microbial Ecology*, 84, 1–14.

Gondard, M., Delannoy, S., Pinarello, V., Aprelon, R., Devillers, E., Galon, C., Pradel, J., Vayssier-Taussat, M., Albina, E., & Moutailler, S. (2020). Upscaling the surveillance of tick-borne pathogens in the French Caribbean islands. *Pathogens*, 9(3), 176.

Grech-Angelini, S., Stachurski, F., Vayssier-Taussat, M., Devillers, E., Casabianca, F., Lancelot, R., Uilenberg, G., & Moutailler, S. (2020). Tick-borne pathogens in ticks (Acari: Ixodidae) collected from various domestic and wild hosts in Corsica (France), a Mediterranean Island environment. *Transboundary and Emerging Diseases*, 67(2), 745–757.

Guimera, R., & Nunes Amaral, L. A. (2005). Functional cartography of complex metabolic networks. *Nature*, 433(7028), 895–900.

Guo, B., Zhang, L., Sun, H., Gao, M., Yu, N., Zhang, Q., Mou, A., & Liu, Y. (2022). Microbial co-occurrence network topological properties link with reactor parameters and reveal importance of low-abundance genera. *npj Biofilms and Microbiomes*, 8(1), 3.

Hajishengallis, G., Darveau, R. P., & Curtis, M. A. (2012). The keystone-pathogen hypothesis. *Nature Reviews Microbiology*, 10(10), 717–725.

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

Krawczyk, A. I., Röttjers, L., Fonville, M., Takumi, K., Takken, W., Faust, K., & Sprong, H. (2022). Quantitative microbial population study reveals geographical differences in bacterial symbionts of *Ixodes ricinus*. *Microbiome*, 10(1), 120.

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

Lejal, E., Marsot, M., Chalvet-Monfray, K., Cosson, J. F., Moutailler, S., Vayssier-Taussat, M., & Pollet, T. (2019). A three-year assessment of *Ixodes ricinus*-borne pathogens in a French peri-urban forest. *Parasites & Vectors*, 12(1), 1–14.

Lhomme, S. (2015). Analyse spatiale de la structure des réseaux techniques dans un contexte de risques. *Cybergeo: European Journal of Geography* [En ligne], Systèmes, Modélisation, Géostatistiques. http://journals.openedition.org/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. *Scientific Reports*, 12(1), 11464.

Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Piloto-Sardiñas, E., Foucault-Simonin, A., Cicculli, V., Moutailler, S., Paoli, J. C., Falchi, A., Obregón, D., & Cabezas-Cruz, A. (2023). Rickettsial pathogens drive microbiota assembly in *Hyalomma marginatum* and *Rhipicephalus bursa* ticks. *Molecular Ecology*, 32(16), 4660–4676.

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

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

Mitri, S., & Richard Foster, K. (2013). The genotypic view of social interactions in microbial communities. *Annual Review of Genetics*, 47, 247–273.

Modlmeier, A. P., Keiser, C. N., Watters, J. V., Sih, A., & Pruitt, J. N. (2014). The keystone individual concept: An ecological and evolutionary overview. *Animal Behaviour*, 89, 53–62.

Narasimhan, S., Rajeevan, N., Liu, L., Zhao, Y. O., Heisig, J., Pan, J., Eppler-Epstein, R., DePonte, K., Fish, D., & Fikrig, E. (2014). Gut microbiota of the tick vector *Ixodes scapularis* modulate colonization of the Lyme disease spirochete. *Cell Host & Microbe*, 15(1), 58–71.

Narasimhan, S., Schuijt, T. J., Abraham, N. M., Rajeevan, N., Coumou, J., Graham, M., Robson, A., Wu, M. J., Daffre, S., Hovius, J. W., & Fikrig, E. (2017). Modulation of the tick gut milieu by a secreted tick protein favors borrelia burgdorferi colonization. *Nature Communications*, 8(1), 184.
Nemergut, D. R., Schmidt, S. K., Fukami, T., O'Neill, S. P., Bilinski, T. M., Stanish, L. F., Knelman, J. E., Darcy, J. L., Lynch, R. C., Wickey, P., & Ferrenberg, S. (2013). Patterns and processes of microbial community assembly. *Microbiology and Molecular Biology Reviews*, **77**(3), 342–356.

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

Ostfeld, R. S., & Keesing, F. (2023). Does experimental reduction of blacklegged tick (*Ixodes scapularis*) abundance reduce Lyme disease incidence? *Pathogens*, **12**(5), 714.

Pereira, D. F. S., Ribeiro, H. S., Gonçalves, A. A. M., da Silva, A. V., Lair, D. F., de Oliveira, D. S., Boas, D. F. V., Conrado, I. D. S. S., Leite, J. C., Barata, L. M., Reis, P. C. C., Mariano, R. M. D. S., Santos, T. A. P., Coutinho, D. C. O., Gontijo, N. F., Araujo, R. N., Galdino, A. S., Paes, P. R. O., Melo, M. M., ... Giunchetti, R. C. (2022). *Rhipicephalus microplus*: An overview of vaccine antigens against the cattle tick. *Ticks and Tick-borne Diseases*, **13**(1), 101828.

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

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

Piloto-Sardiñas, E., Cano-Argüelles, A. L., Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Corduneanu, A., Obregón, D., Oleaga, A., Pérez-Sánchez, R., & Cabezas-Cruz, A. (2023). Comparison of salivary gland and midgut microbiome in the soft ticks *Ornithodoros erraticus* and *Ornithodoros moubata*. *Frontiers in Microbiology*, **14**, 1173609.

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ąc, 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.

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

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

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

Rodríguez, S. D., Ortiz, M. A. G., Ocampo, R. J., & Vega y Murguía, C. A. (2009). Molecular epidemiology of bovine anaplasmosis with a particular focus in Mexico. *Infection, Genetics and Evolution*, **9**(6), 1092–1101.

Röttjers, L., Vandeputte, D., Raes, J., & Faust, K. (2021). Null-model-based network comparison reveals core associations. *ISME Communications*, **1**(1), 36.

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

Salinas-Estrella, E., Amaro-Estrada, I., Cobaxin-Cárdenas, M. E., Preciado de la Torre, J. F., & Rodríguez, S. D. (2022). Bovine anaplasmosis: Will there ever be an almighty effective vaccine? *Frontiers in Veterinary Science*, **9**, 946545.

Stolyar, S., Van Dien, S., Hillesland, K. L., Pinel, N., Lie, T. J., Leigh, J. A., & Stahl, D. A. (2007). Metabolic modeling of a mutualistic microbial community. *Molecular Systems Biology*, **3**(1), 92.

Swei, A., & Kwan, J. Y. (2017). Tick microbiome and pathogen acquisition altered by host blood meal. *The ISME Journal*, **11**(3), 813–816.

Takumi, K., Sprong, H., & Hofmeester, T. R. (2019). Impact of vertebrate communities on *Ixodes ricinus*-borne disease risk in forest areas. *Parasites & Vectors*, **12**(1), 1–12.

Tonk-Rügen, M., Zając, Z., & Cabezas-Cruz, A. (2023). Can tick microbiome explain nonlinear relationship between tick abundance and Lyme disease incidence? *Pathogens*, **12**(10), 1229.

Wei, Z., Yang, T., Friman, V. P., Xu, Y., Shen, Q., & Jousset, A. (2015). Trophic network architecture of root-associated bacterial communities determines pathogen invasion and plant health. *Nature communications*, **6**, 8413. https://doi.org/10.1038/ncomms9413

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

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

Wu-Chuang, A., Mateos-Hernandez, L., Maitre, A., Rego, R. O., Šíma, R., Porcelli, S., Rakotobe, S., Foucault-Simonin, A., Moutailler, S., Palinauskas, V., Aželytė, J., Šimo, L., Obregón, D., & Cabezas-Cruz, A. (2023). Microbiota perturbation by anti-microbiota vaccine reduces the colonization of *Borrelia afzelii* in *Ixodes ricinus*. *Microbiome*, **11**(1), 151.

Zając, Z., Obregón, D., Foucault-Simonin, A., Wu-Chuang, A., Moutailler, S., Galon, C., Kulisz, J., Woźniak, A., Bartosik, K., & Cabezas-Cruz, A. (2023). Disparate dynamics of pathogen prevalence in *Ixodes ricinus* and *Dermacentor reticulatus* ticks occurring sympatrically in diverse habitats. *Scientific Reports*, **13**(1), 10645.

Zivkovic, Z., Esteves, E., Almazán, C., Daffre, S., Nijhof, A. M., Kocan, K. M., Jongejan, F., & De la Fuente, J. (2010). Differential expression of genes in salivary glands of male *Rhipicephalus (Boophilus) microplus* in response to infection with *Anaplasma marginale*. *BMC Genomics*, **11**, 1–12.

## SUPPORTING INFORMATION

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

**How to cite this article:** 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., Obregón, D., Svobodová, K., Wu-Chuang, A., & Cabezas-Cruz, A. (2024). Dynamic nesting of *Anaplasma marginale* in the microbial communities of *Rhipicephalus microplus*. *Ecology and Evolution*, **14**, e11228. https://doi.org/10.1002/ece3.11228
Nemergut, D. R., Schmidt, S. K., Fukami, T., O'Neill, S. P., Bilinski, T. M., Stanish, L. F., Knelman, J. E., Darcy, J. L., Lynch, R. C., Wickey, P., & Ferrenberg, S. (2013). Patterns and processes of microbial community assembly. *Microbiology and Molecular Biology Reviews*, *77*(3), 342–356.

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

Ostfeld, R. S., & Keesing, F. (2023). Does experimental reduction of blacklegged tick (*Ixodes scapularis*) abundance reduce Lyme disease incidence? *Pathogens*, *12*(5), 714.

Pereira, D. F. S., Ribeiro, H. S., Gonçalves, A. A. M., da Silva, A. V., Lair, D. F., de Oliveira, D. S., Boas, D. F. V., Conrado, I. D. S. S., Leite, J. C., Barata, L. M., Reis, P. C. C., Mariano, R. M. D. S., Santos, T. A. P., Coutinho, D. C. O., Gontijo, N. F., Araujo, R. N., Galdino, A. S., Paes, P. R. O., Melo, M. M., ... Giunchetti, R. C. (2022). *Rhipicephalus microplus*: An overview of vaccine antigens against the cattle tick. *Ticks and Tick-borne Diseases*, *13*(1), 101828.

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

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

Piloto-Sardiñas, E., Cano-Argüelles, A. L., Maitre, A., Wu-Chuang, A., Mateos-Hernández, L., Corduneanu, A., Obregón, D., Oleaga, A., Pérez-Sánchez, R., & Cabezas-Cruz, A. (2023). Comparison of salivary gland and midgut microbiome in the soft ticks *Ornithodoros erraticus* and *Ornithodoros moubata*. *Frontiers in Microbiology*, *14*, 1173609.

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ąc, 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.

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

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

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

Rodríguez, S. D., Ortiz, M. A. G., Ocampo, R. J., & Vega y Murguía, C. A. (2009). Molecular epidemiology of bovine anaplasmosis with a particular focus in Mexico. *Infection, Genetics and Evolution*, *9*(6), 1092–1101.

Röttjers, L., Vandeputte, D., Raes, J., & Faust, K. (2021). Null-model-based network comparison reveals core associations. *ISME Communications*, *1*(1), 36.

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

Salinas-Estrella, E., Amaro-Estrada, I., Cobaxin-Cárdenas, M. E., Preciado de la Torre, J. F., & Rodríguez, S. D. (2022). Bovine anaplasmosis: Will there ever be an almighty effective vaccine? *Frontiers in Veterinary Science*, *9*, 946545.

Stolyar, S., Van Dien, S., Hillesland, K. L., Pinel, N., Lie, T. J., Leigh, J. A., & Stahl, D. A. (2007). Metabolic modeling of a mutualistic microbial community. *Molecular Systems Biology*, *3*(1), 92.

Swei, A., & Kwan, J. Y. (2017). Tick microbiome and pathogen acquisition altered by host blood meal. *The ISME Journal*, *11*(3), 813–816.

Takumi, K., Sprong, H., & Hofmeester, T. R. (2019). Impact of vertebrate communities on *Ixodes ricinus*-borne disease risk in forest areas. *Parasites & Vectors*, *12*(1), 1–12.

Tonk-Rügen, M., Zając, Z., & Cabezas-Cruz, A. (2023). Can tick microbiome explain nonlinear relationship between tick abundance and Lyme disease incidence? *Pathogens*, *12*(10), 1229.

Wei, Z., Yang, T., Friman, V. P., Xu, Y., Shen, Q., & Jousset, A. (2015). Trophic network architecture of root-associated bacterial communities determines pathogen invasion and plant health. *Nature communications*, *6*, 8413. https://doi.org/10.1038/ncomms9413

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

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

Wu-Chuang, A., Mateos-Hernandez, L., Maitre, A., Rego, R. O., Šíma, R., Porcelli, S., Rakotobe, S., Foucault-Simonin, A., Moutailler, S., Palinauskas, V., Aželytė, J., Šimo, L., Obregón, D., & Cabezas-Cruz, A. (2023). Microbiota perturbation by anti-microbiota vaccine reduces the colonization of *Borrelia afzelii* in *Ixodes ricinus*. *Microbiome*, *11*(1), 151.

Zając, Z., Obregón, D., Foucault-Simonin, A., Wu-Chuang, A., Moutailler, S., Galon, C., Kulisz, J., Woźniak, A., Bartosik, K., & Cabezas-Cruz, A. (2023). Disparate dynamics of pathogen prevalence in *Ixodes ricinus* and *Dermacentor reticulatus* ticks occurring sympatrically in diverse habitats. *Scientific Reports*, *13*(1), 10645.

Zivkovic, Z., Esteves, E., Almazán, C., Daffre, S., Nijhof, A. M., Kocan, K. M., Jongejan, F., & De la Fuente, J. (2010). Differential expression of genes in salivary glands of male *Rhipicephalus (Boophilus) microplus* in response to infection with *Anaplasma marginale*. *BMC Genomics*, *11*, 1–12.

## SUPPORTING INFORMATION

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

**How to cite this article:** 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., Obregón, D., Svobodová, K., Wu-Chuang, A., & Cabezas-Cruz, A. (2024). Dynamic nesting of *Anaplasma marginale* in the microbial communities of *Rhipicephalus microplus*. *Ecology and Evolution*, *14*, e11228. https://doi.org/10.1002/ece3.11228