Disruption of bacterial interactions and community assembly in Babesia-infected Haemaphysalis longicornis following antibiotic treatment

Abstract

Background

A previous study highlighted the role of antibiotic-induced dysbiosis in the tick microbiota, facilitating the transstadial transmission of Babesia microti from nymph to adult in Haemaphysalis longicornis. This study builds on previous findings by analyzing sequence data from an earlier study to investigate bacterial interactions that could be linked to enhanced transstadial transmission of Babesia in ticks. The study employed antibiotic-treated (AT) and control-treated (CT) Haemaphysalis longicornis ticks to investigate shifts in microbial community assembly. Network analysis techniques were utilized to assess bacterial interactions, comparing network centrality measures between AT and CT groups, alongside studying network robustness and connectivity loss. Additionally, functional profiling was conducted to evaluate metabolic diversity in response to antibiotic treatment.

Results

The analysis revealed notable changes in microbial community assembly in response to antibiotic treatment. Antibiotic-treated (AT) ticks displayed a greater number of connected nodes but fewer correlations compared to control-treated (CT) ticks, indicating a less interactive yet more connected microbial community. Network centrality measures such as degree, betweenness, closeness, and eigenvector centrality, differed significantly between AT and CT groups, suggesting alterations in local network dynamics due to antibiotic intervention. Coxiella and Acinetobacter exhibited disrupted connectivity and roles, with the former showing reduced interactions in AT group and the latter displaying a loss of connected nodes, emphasizing their crucial roles in microbial network stability. Robustness tests against node removal showed decreased stability in AT networks, particularly under directed attacks, confirming a susceptibility of the microbial community to disturbances. Functional profile analysis further indicated a higher diversity and richness in metabolic capabilities in the AT group, reflecting potential shifts in microbial metabolism as a consequence of antimicrobial treatment.

Conclusions

Our findings support that bacterial interaction traits boosting the transstadial transmission of Babesia could be associated with reduced colonization resistance. The disrupted microbial interactions and decreased network robustness in AT ticks suggest critical vulnerabilities that could be targeted for managing tick-borne diseases.

Background

Ticks are obligate hematophagous arthropods that play a significant role in the transmission of various infectious agents, such as bacteria (e.g., Borrelia and Anaplasma), viruses (e.g., tick-borne encephalitis virus, TBEV), and protozoan parasites (e.g., Babesia and Theileria) [1]. The Asian long-horned tick, Haemaphysalis longicornis, has been associated with more than 30 human and animal pathogens, raising medical and veterinary concerns [2]. Originating from Eastern Asia, this tick has the capability to rapidly spread into new areas, presenting an emerging disease threat (e.g., recent invasion of USA) [2]. Notably, H. longicornis has been identified as a vector for Babesia [3]. Among Babesia pathogenic species, B. microti has gained significant recognition as the primary etiological agent responsible for babesiosis in humans, particularly in the USA [4, 5]. This hemoparasite can infect small rodents like mice and voles, serving as the primary reservoir hosts for B. microti.

The typical mode of Babesia transmission involves transovarial transmission through adult tick to their offspring [6, 7]. Transovarial transmission enhances species diversification by facilitating host switching to other vertebrate hosts [8]. In addition to transovarial transmission, Babesia pathogens exhibit transstadial transmission [9]. However, transovarial transmission of B. microti has been found absent in ticks belonging to the genus Ixodes, including Ixodes ricinus [10] as well as in ticks outside the Ixodes genus, such as Rhipicephalus haemaphysalis [11]. Consequently, this pathogen relies exclusively on transstadial transmission after acquisition from an infected host [10]. Furthermore, Gray et al. [10] experimentally demonstrated that while transstadial transmission occurs, the parasite does not persist after molting.

Combination therapies, typically involving an antiprotozoal agent and an antibiotic, were recommended and applied for the treatment of human babesiosis [12]. Antibiotics ingested with the tick's blood meal have the potential to disrupt the tick microbiota, leading to tick gut dysbiosis [5]. Ticks harbor a range of pathogenic microorganisms alongside endosymbionts and commensals, emphasizing the impact of these microbes on tick fitness and pathogen transmission [13, 14]. Microbiota dysbiosis may either reduce [15] or increase [16] the tick's susceptibility to tick-borne pathogens (TBPs), potentially shaping tick-borne diseases (TBDs) ecology. In the context of B. microti, gut microbiota has been found to play a role in facilitating the transstadial transmission of this pathogen from nymphs to adult H. longicornis ticks [16]. Furthermore, in their study, Wei et al. [16] found that antibiotics administered to mouse hosts altered the microbiota of adult ticks, with significant differences in the abundance of Coxiella and Acinetobacter between the antibiotic-treated group (AT) and the untreated control group (CT). Coxiella was the most abundant genus in CT adults, whereas Acinetobacter dominated in AT adults, indicating a shift in the microbiota composition due to antibiotics [16]. Similar evidence has been observed in Plasmodium falciparum parasites, where antibiotics have been shown to increase parasite colonization in mosquitoes [17], an effect mediated by the vector microbiota [17]. This highlights the significance of tick-microbiota interactions for vector competence, including the modulation of tick vector capacity by influencing pathogen colonization of tick tissues [18]. Overall, the evidence suggests a potential role of colonization resistance in the context of TBPs [19] where the microbiota within ticks may resist the colonization or transmission of TBPs like B. microti [16], a protection disrupted by antibiotic-induced dysbiosis [20].

Colonization resistance is the phenomenon where established microbial communities prevent the invasion and establishment of new, often pathogenic, species [21, 22, 242523]. This study aims to build upon the findings of Wei et al. [16] by identifying bacterial interaction traits associated with reduced colonization resistance and enhanced transstadial transmission of Babesia in ticks. While microbial diversity traits linked to colonization resistance have been extensively studied in vectors like ticks [16] and mosquitoes [17], employing a network-based approach to examine the impact of antibiotics on tick microbiota offers a comprehensive method for evaluating factors affecting microbial community interactions, structure, and functionality. Microbes with co-occurrence patterns are known to be influenced by metabolic interactions and competition for resources [26, 27], offering the potential to capture crucial community characteristics that may not be revealed in analyses based solely on microbial diversity or abundance [27, 28].

Network approach enables the detailed mapping of community interactions [29], potentially revealing how antibiotics disrupt normal microbial relationships essential for pathogens such as Babesia. For example, it can identify how altered interactions of key taxa like Coxiella and Acinetobacter may compromise the community's ability to resist Babesia colonization. By assessing changes in metrics such as modularity and centrality, network analysis provides quantitative evidence of the effects of disturbance on microbial community organization [30]. For instance, increased modularity and altered centrality in antibiotic-treated ticks may indicate a less compartmentalized with strong separation between communities and potentially less resilient microbial network, facilitating Babesia establishment and proliferation.

Furthermore, networks allow the evaluation of emerging system properties such as the robustness of microbial communities [30], which allows deciphering how changes impact community structure and function, particularly in response to the loss [29, 31] and addition [31] of specific nodes or bacterial taxa. Ultimately, this approach not only enhances our comprehension of how antibiotics indirectly influence Babesia dynamics by altering host microbiota but may also contribute to predictive modelling of disease transmission, offering crucial insights for the effective management of TBDs. Therefore, this study aims to explore the repercussions of antibiotic-induced dysbiosis on microbial interactions within Haemaphysalis longicornis ticks. Specifically, it targets the identification of bacterial interaction traits associated with colonization resistance and heightened transstadial transmission of B. microti. Employing a network-based approach, the study intends to comprehensively map and analyze the structural, functional, and interactive aspects of tick microbiota under antibiotic exposure. The research holds a scientific merit by elucidating the impact of antibiotics on tick microbiota and subsequent pathogen transmission dynamics. It advances our understanding on the impacts of antibiotics on tick microbiota and subsequent B. microti transmission dynamics.

Methods

Original 16S rRNA datasets

In this study, we used available 16S rRNA amplicon sequence datasets generated by Wei et al. [16]. Sequencing was performed on an Illumina MiSeq system and the resultant data has been deposited in the National Center for Biotechnology Information (NCBI)’s GenBank, under Sequence Read Archive (SRA) accession numbers SRP322057 and SRP323180. The raw data was collected as part of a study evaluating the effects of antibiotics on the microbiota of H. longicornis ticks and their implications for B. microti transmission, comparing the microbiota of nymphal ticks fed on Babesia microti-infected mice treated with antibiotics (AT) with the control-treated group (CT).

Analysis of 16S rRNA amplicon sequences

To facilitate a comprehensive exploration of the tick microbiota's composition and function, we retrieved 16S rRNA sequences from the Sequence Read Archive (SRA) repository within the Quantitative Insights Into Microbial Ecology (QIIME2) 2023.7 environment. Briefly, the q-fondue script, following the methodology outlined by Bolyen et al. [32], facilitated the download process. The sequence data (demultiplexed fastq files) underwent denoising, quality trimming, merging, chimera removal, and filtering using DADA2 software [33] implemented within the QIIME2 [32]. Taxonomic assignment to Amplicon Sequence Variants (ASVs) was accomplished using the Classify-Sklearn Naive Bayes method based on the 16S rRNA SILVA database v.138 [34]. The resulting taxonomic table was collapsed at the genus level, and subsequently employed for network analysis and pathway prediction.

Bacterial co-occurrence networks analysis

Co-occurrence networks were constructed for the CT and AT groups using taxonomic information at the genus level to investigate the effect of antibiotic treatment on community assembly. The prevalence of Coxiella and Acinetobacter within the microbiota in CT and AT groups respectively [16] suggested a possible influence of both taxa on community assembly and stability in response to antibiotic treatment. In addition, global co-occurrence networks were constructed for each condition (CT and AT) after the individual removal of Coxiella and Acinetobacter to assess their specific impacts. In the graphical visualization of microbial community assemblies, bacterial taxa are represented by nodes and the significant interactions between taxa are represented by edges. Sparse Correlations for Compositional data (SparCC) method [35] implemented in the SpiecEasi R package [36] analysis was used to identify significant positive (weight > 0.75) or negative (weight < -0.75) interactions between taxa, Gephi software 0.10.0 [37] was employed to visualize and analyze network features [i.e. number of nodes and edges, network diameter, modularity, average degree, weighted degree, clustering coefficient and centrality]. Additionally, the Core Association Network (CAN) analysis identified common microbial associations between CT and AT networks, conducted using the Anuran toolbox [38] within the Anaconda Python environment (Anaconda Software Distribution, 2023).

Comparative network analysis

The network comparisons between the same taxa in two different bacterial networks were conducted using the package Network Construction and Comparison for Microbiome Data (NetCoMi) [39] in R v.4.0.3 (R Core Team, 2023) [40], and performed using RStudio (RStudio Team, 2020) [41]. By comparing networks: (i) between CT and AT groups, (ii) between CT and AT groups after individual removal of Coxiella (woC) and Acinetobacter (woA), and (iii) within the same CT or AT group in the presence and after the individual removal of both taxa, the study aims to understand how these taxa influence network structure. The Jaccard index was computed to assess dissimilarities in local centrality measures (degree, betweenness centrality, closeness centrality, and eigenvector centrality) between the two networks [i.e., CT vs. AT; CT vs. CT (woC); CT (woC) vs. AT (woC)], providing insights into the impact of antibiotic treatment and taxa presence and absence on network structure. This index evaluates the similarity between sets of 'most central nodes' in networks, defined as nodes with centrality values above the empirical 75% quartile. It ranges from 0 (totally different) to 1 (unique). The associated p-values, P(J ≤ j) and P(J ≥ j), indicate the probability that the observed Jaccard index value is 'less than or equal to' or ‘higher than or equal’ to' the Jaccard value expected at random, considering the total number of taxa in both sets [42]. The ARI was calculated to test the dissimilarity of clustering in the networks. 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 clustering. The p-value test if the calculated value is significantly different from zero [39].

Local connectivity of Acinetobacter and Coxiella in the microbial community

The effect of antibiotic treatment on specific taxa was analyzed by determining the direct associations of Acinetobacter and Coxiella with the rest of the bacterial microbiota in the CT and AT groups. This analysis aimed to understand how these taxa influenced the overall microbial community structure. Sub-networks representing the positive and negative associations of local connectivity were constructed in Gephi 0.10.0 [37], with connections between microbes quantified using SparCC (SparCC > 0.75 or < -0.75) as implemented in the SpiecEasi R package [36], enabling the identification of key microbial associations influenced by antibiotic treatment.

Centrality measures and module dynamics in networks of CT and AT microbiota

To understand the networks structure of the CT and AT groups, the centrality measures distribution among different taxa, including Coxiella and Acinetobacter, were analyzed. The function of each taxon within the network was assessed using measures of within-module (Zi) and among-module (Pi) connectivity [43]. The measures allowed the taxa to be categorized into four different roles according to their connectivity: (i) network hubs: taxa that serve as central connectors both within their module and across the entire network, with high connectivity (Zi > 2.5 and Pi > 0.62), (ii) module hubs: highly connected taxa within their own module but not significantly connected to other modules (Zi > 2.5 and Pi ≤ 0.62), (iii) connectors: taxa that primarily link different modules together, indicating their role in connecting disparate parts of the network (Zi ≤ 2.5 and Pi > 0.62), and (iv) peripheral taxa: taxa with limited connections within the module and minimal interactions with other modules (Zi ≤ 2.5 and Pi ≤ 0.62). The Zi and Pi values were calculated using the "code-zi-pi-plot" package [44, 45] in R (R Core Team, 2023) [40], considering only positive interactions. Visualization was done using GraphPad Prism version 8.0.1 (GraphPad Software, San Diego, California, USA), enabling the analysis of taxon positions within the networks. This analysis helps identify major key players influencing community assembly in response to antibiotic treatment.

Identification of keystone taxa

To pinpoint microbial taxa that play crucial roles in maintaining the stability and function of the tick microbiota, keystone taxa were identified based on three criteria as described by Mateos-Hernandez et al. [46]: (i) ubiquitousness (microbial taxa present across all the samples of an experimental group), (ii) eigenvector centrality higher than 0.75, and (iii) abundance higher than the mean clr value (i.e., higher than that of the mean relative abundance of all taxa in an experimental group).

Network robustness analysis using node removal and addition

The stability of the CT and AT networks was evaluated and compared to determine if antibiotic treatment impacts the stability and resilience of the tick microbiota network in response to node disturbances. The robustness of the networks was compared under two types of disturbances: node removal and addition. The comparisons were made between: (i) CT and AT groups, (ii) CT and AT groups after individual removal of Coxiella and Acinetobacter, and (iii) within the same CT or AT group in the presence and after the individual removal of each taxon. To evaluate network resistance to node removal, an attack tolerance test was conducted using Network Strengths and Weaknesses Analysis (NetSwan) [47] in R v.4.0.3 (R Core Team, 2023) [40], performed using the RStudio (RStudio Team, 2020) [41]. The networks underwent both random and directed attacks. For directed node removal, three scenarios were employed: (i) a directed attack removing nodes in decreasing order of their betweenness centrality (BNC) value, (ii) a cascading attack recalculating BNC values after each node removal, and (iii) a degree centrality removal prioritizing nodes with the highest degree centrality values.

Conversely, a node addition analysis was performed following the approach outlined by Freitas et al. [48] in R v.4.3.1 (R Core Team, 2023), and performed using the RStudio (RStudio Team, 2020) [41]. In this analysis, new nodes were randomly connected to the existing network, and resultant changes were quantified by evaluating the size of the Largest Connected Component (LCC) and the Average Path Length (APL). To enhance the precision of the network's robustness assessment, multiple simulations were conducted with varying sets of nodes, introducing 500, 700, and 1000 nodes. The outcomes were graphically represented using GraphPad Prism 9.0.2 (GraphPad Software Inc., San Diego, CA, USA). For all comparisons made between the control (CT, CT (woC) and CT (woA)) and antibiotic-treated (AT, AT (woC) and AT (woA)) groups, a delta value was calculated. For node removal analysis, it was the difference in the fraction of nodes needed to achieve a connectivity loss of 0.80. For the addition of nodes, delta values for LCC size and APL were calculated by subtracting values in the CT network from those in the AT network after adding 100 and 1000 nodes.

Functional profile prediction

A step-forward analysis was performed to predict microbial functional traits, specifically enzymatic pathways, utilizing PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) standalone version [49]. Gene catalogues, including Kyoto Encyclopedia of Genes and Genomes (KEGG), Orthologs (KO), Enzyme Classification numbers (EC), and Cluster of Orthologous Genes (COGs) [50], along with the MetaCyc database [51], were employed to annotate major pathway categories and facilitate mapping. Following the output table, the taxa's contribution to predicted metabolic pathways was investigated. To ensure robust statistical analysis, various methods were employed. Initially, alpha diversity was assessed using observed features [52] and Pielou’s evenness metrics [53] via the q2-diversity method in QIIME2 plugin. Subsequently, differences in pathway frequency were evaluated using the DESeq2 package [54] in R v.4.0.3 (R Core Team, 2023) [40], enabling the identification of statistically significant alterations in pathway abundance between the CT and AT groups. This analysis resulted in a Volcano plot with Benjamini correlation, providing a visual representation of the significance and magnitude of pathway abundance changes. Analyses were performed using the RStudio Integrated Development Environment (IDE) v.2023.03.0-daily + 82.pro2 (RStudio Team, 2020) [41].

Results

Changes in microbial community assembly and node centrality in response to antibiotic treatment

The potential impact of antibiotics on the assembly of microbial communities in ticks’ microbiota and the infection of Babesia was investigated, specifically examining changes beyond bacterial composition, richness, and relative abundance. To explore this, microbial co-occurrence patterns through co-occurrence networks were analyzed. It was observed that AT exhibited a greater number of connected nodes compared to CT, showing topological variations between both groups (Table 1). In contrast, CT displayed the highest number of correlations compared to AT (Fig. 1 A-B, Table 1). However, a similar balance between positive and negative correlations was observed in both groups (Table 1), suggesting that antibiotic treatment reduces the level of interaction within the community but does not change its nature. Additionally, CT displayed lower modularity and diameter values than AT (Table 1). The core association network (CAN) revealed 71 core associated nodes between CT and AT networks (Fig. 1C) demonstrating variability in community configuration in these two conditions. NetCoMi was used to test dissimilarities between local network centrality measures of the CT and AT networks, Jaccard index was calculated for degree, betweenness centrality, closeness centrality and eigenvector centrality (Jacc = 0, lowest similarity and Jacc = 1, highest similarity). All these measures between the two networks were found to be lower than expected by random (P (≤ Jacc) < 0.05, Table 2).

Table 1 Topological features of taxonomic co-occurrence networks

Table 2 Jaccard index for comparison between CT and AT network

Fig. 1

Antibiotic treatment effects on microbiota diversity and community assembly in Babesia-infected ticks. Global co-occurrence networks of (A) CT and (B) AT networks. C Core Association Network between CT/AT networks. Coxiella’s local connectivity in (D) CT and (E) AT networks. F Direction of associations of common direct neighbor to Coxiella between CT and AT groups. Acinetobacter’s local connectivity in (G) CT and (H) AT networks. I Direction of associations of common direct neighbor to Acinetobacter between CT and AT groups. Within-module and among-module connectivity, Zi-Pi plot of the individual genera from (J) CT and (K) AT groups. Only nodes with at least one significant correlation are represented. Node colors are based on modularity class metric and equal color means modules of co-occurring taxa. The size of the nodes is proportional to the eigenvector centrality of each taxon. The colors in the edges represent strong positive (green) or negative (red) correlations (SparCC > 0.75 or <  − 0.75)

The most abundant genera in the CT adults were Ammoniphilus and Coxiella. However, the gut microbiota in AT adults was dominated by Acinetobacter [16]. To understand the role of Coxiella and Acinetobacter in the community assembly, the local connectivity in the CT and AT groups was visually inspected. Both networks displayed two principal modules with negative and positive co-occurring interactions between their nodes. In the CT group's microbial community of infected ticks with Babesia, Coxiella occupied a central position, displaying numerous interactions with other genera (Fig. 1D, Supplementary Table S1). However, visual inspection of the networks showed that the treatment caused a shift in the bacterial community assembly patterns with a notable reduction in Coxiella's co-occurrence network with other taxa (Fig. 1E, Supplementary Table S1). Although the local connectivity sub-networks between CT and AT groups were unique in their connections, Coxiella maintained common direct association with four taxa in both networks (Fig. 1F). Similarly, Acinetobacter in CT showed high centrality and extensive connections with various taxa (Fig. 1G, Supplementary Table S1). Antibiotics disrupted Acinetobacter's interactions, reducing its connectivity and affecting the microbial community co-occurrence network resulting in a loss of a few connected nodes (Fig. 1H, Supplementary Table S1). However, Acinetobacter maintained a common direct association with eight taxa in both groups (Fig. 1I). Markedly, Coxiella and Acinetobacter presented a direct association of co-exclusion, which was maintained in the AT group (Fig. D-I, Supplementary Table S1).

When analyzing the distribution of connections (Zi and Pi connectivity), a similar pattern was observed in both CT (Fig. 1J) and AT (Fig. 1K) networks. All taxa, including Coxiella and Acinetobacter, were classified as peripheral, indicating that they do not function as central hubs within their respective networks. Despite the absence of central nodes within CT and AT networks, eight taxa met the criteria to be considered keystone taxa for the CT and six for the AT group microbiota (Table 3). In general, the topological variations, the changed clustering patterns of Coxiella and Acinetobacter and the absence of other central hubs demonstrate notable differences in community assembly between both conditions and a considerable susceptibility to antibiotic treatment. This behavior suggests a possible influence of both taxa on community assembly as support for the destabilizing effect that antimicrobial treatment can cause.

Table 3 Keystone taxa of the bacterial communities of CT and AT groups

Influence of Coxiella and Acinetobacter on the assembly and hierarchy of the microbiota in response to antibiotic treatment

To investigate the impact of Coxiella and Acinetobacter on community assembly in response to antibiotic treatment, the topology of the network was analyzed after individual removal of both taxa (woC and woA), comparing the CT and AT groups (Table 1). Removal of Coxiella resulted in a loss of interactions in both CT (woC) and AT (woC) networks (Fig. 2 A-B, Table 1) compared to conditions with the taxon present (Table 1), suggesting potential stability conferred by the bacteria in its endosymbiont condition. Similar to those observed in the presence of Coxiella, there was a balance between positive and negative interactions in the CT (woC) and AT (woC) networks, with cooperation being slightly greater (Fig. 2 A-B, Table 1). Moreover, after removing Acinetobacter, a decrease in the number of interactions was observed in both CT (woA) and AT (woA) (Fig. 2 C-D, Table 1) compared to conditions with the taxon present. However, a disproportion between positive and negative associations was evident in AT (woA) compared to CT (woA) (Fig. 2 C-D, Table 1), indicating increased competition post-taxa removal. Furthermore, the decrease in interactions of cooperation and co-exclusion (Table 1), upon comparing networks of the same condition after individual bacteria removal suggests a reconfiguration of the remaining taxa.

Fig. 2

Comparative analysis of microbial co-occurrence networks in CT and AT groups: Impact of Coxiella and Acinetobacter removal. Global co-occurrence networks after Coxiella’s removal of (C) CT(woC) and (D) AT(woC). Global co-occurrence networks after Acinetobacter’s removal of (E) CT(woA) and (F) AT(woA). The colors in the edges represent positive (green) or negative (red) correlations. In global networks, the 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 impact of individual removal of Coxiella (woC) and Acinetobacter (woA) on the interaction patterns between CT and AT networks (Table 4) were compared, as well as within the same condition before and after in-silico manipulation of the taxa (Supplementary Table S2 and S3). Jaccard index values were lower than expected by random (P (≤ Jacc) < 0.05) in woC (CT vs. AT) and woA (CT vs. AT) comparison, respectively (Table 4). Moreover, the centrality measures were higher than expected by random (P (≥ Jacc) < 0.05) in CT vs. CT (woC) (Supplementary Table S2) and CT vs. CT (woA) (Supplementary Table S3) comparison, respectively. Jaccard index values of the AT vs. AT (woC) and AT vs. AT (woA) comparisons, were higher than expected by random (P (≥ Jacc) < 0.05) except for betweenness centrality and closeness centrality which had a random distribution. Comparing node clustering for woC (CT vs. AT) and woA (CT vs. AT) revealed low ARI values, indicating the antibiotic treatment strongly influences microbial community clustering patterns (Table 5). In contrast, ARI values were close to 1 within the same condition (CT and AT) after individual taxon removal (Table 5). High ARI values indicated strong similarities and suggested that despite the topological variations (Table 1) the impact of in-silico manipulation of Coxiella and Acinetobacter on clustering patterns under the same conditions is moderated.

Table 4 Jaccard index for comparison between CT and AT groups after Coxiella and Acinetobacter ‘s removal

Table 5 Network clustering comparisons

Changes in network robustness in response to antibiotic treatment

The robustness of the CT and AT co-occurrence networks was compared to determine if the effect of antibiotic treatment compromises stability against node removal and addition. Assessment of connectivity loss showed that directed attacks (degree, cascading, and betweenness methods) (Fig. 3A-B and Supplementary Fig. S1A) had more significant impacts on both CT and AT networks compared to the random method (Supplementary Fig. S1B, Supplementary Table S4). During node removal, delta values were calculated by subtracting the fraction of nodes necessary to achieve a connectivity loss of 0.8 of CT minus AT network. Degree-directed node removal yielded a negative delta value, indicating AT's greater robustness compared to CT (Fig. 3A, Supplementary Table S4). Conversely, cascading and betweenness node removal showed positive values, suggesting antibiotic treatment reduces network stability against both types of directed attacks (Fig. 3B, Supplementary Fig. S1A, Supplementary Table S4). The delta value equal to 0 confirms that, regardless of the possible community-destabilizing effect of antibiotic treatment, the robustness of the CT and AT networks is not affected by random attacks (Supplementary Fig. S1B, Supplementary Table S4). Comparing robustness against node removal between CT and AT networks after individual Coxiella (Fig. 3C-D, Supplementary Fig. S1C-D) and Acinetobacter (Fig. 3E-F, Supplementary Fig. S1E-F) removal showed that only in the absence of Coxiella and after directed attack in degree and betweenness, the AT network was more robust than the CT network (Fig. 3C and Supplementary Fig. S1C). Additionally, positive delta values in the absence of Acinetobacter and after directed attack in degree and cascading (Fig. 3E-F) highlight the susceptibility of the AT group to external perturbations and suggest a possible stabilizing potential exerted by the taxon when present (Supplementary Table S4).

Fig. 3

Robustness comparison after removal and addition of nodes between the CT and AT groups. Connectivity loss measured after directed (degree and cascading) attack in CT and AT networks: (A) CT/AT (degree), and (B) CT/AT (cascading), Connectivity loss measured against directed (degree and cascading) attack between the CT and AT networks after Coxiella’ removal: (C) CT(woC)/AT(woC) (degree), and (D) CT(woC)/AT(woC) (cascading). Connectivity loss measured against directed (degree and cascading) attack between the CT and AT networks after Acinetobacter’s removal: (E) CT(woA)/AT(woA) (degree), and (F) CT(woA)/AT(woA) (cascading). Robustness comparison between CT and AT networks after the addition of nodes. The largest connected component (LCC) and average path length (APL) values are represented and compared between CT and AT networks: (G) CT/AT (LCC), and (H) CT/AT (APL). Robustness comparison between CT and AT networks after Coxiella’s removal: (I) CT(woC)/AT(woC) (LCC), and (J) CT(woC)/AT(woC) (APL). Robustness comparison between CT and AT networks after Acinetobacter’s removal: (K) CT(woA)/AT(woA) (LCC), and (L) CT(woA)/ AT(woA) (APL)

For node addition, 1000 nodes were added, and two key network properties were quantified: LCC (Fig. 3G, Supplementary Table S5) and APL (Fig. 3H, Supplementary Table S5). Applying the nodes addition strategy increased the LCC for the AT group compared to the CT group (Fig. 3G, Supplementary Table S5). Specifically, The AT group's average robustness increased at the 750th node addition, followed by a similar rise in the CT group at the 800th node addition. Both groups showed a rapid increase in robustness by the 900th node added. Regarding the APL test, initially, the AT group exhibited slightly higher values than the CT group until the addition of the 800th node, after which the CT network rapidly increased (Fig. 3H, Supplementary Table S5). Ultimately, APL values for both groups converged to approximately equal values upon the addition of 1000th node. However, the observed overlap and increase in APL values in both CT and AT networks suggest that the addition of nodes causes a loss of functional connectivity within the networks independently of the effect of antibiotic treatment (Fig. 3H, Supplementary Table S5). The higher LCC and APL values of AT compared to CT were highlighted by negative delta values after adding the 100th and 1000th nodes (Fig. 3G-H, Supplementary Table S5). This persisted after Coxiella (Fig. 3I-J) and Acinetobacter (Fig. 3K-L) removal, with higher LCC and APL values in AT networks Fig. 3I-L, Supplementary Table S5), except for LCC values after adding the 1000th nodes in Coxiella-removed networks, showing a positive delta value (Fig. 3I, Supplementary Table S5).

Influence of the removal of Coxiella and Acinetobacter on network robustness within the CT and AT groups

To demonstrate the direct influence of both bacteria within the same condition (CT or AT) after node removal, the robustness of the networks was analyzed and compared after individual removal of Coxiella (woC) (Fig. 4A-D, Supplementary Fig. S2A-D) and Acinetobacter (woA) (Fig. 4E-H, Supplementary Fig. S2E-H). The removal of Coxiella did not modify the robustness of the networks within the same condition against directed (Fig. 4A-D) and random attacks (Supplementary Fig. S2A-D). In contrast, the removal of Acinetobacter conferred robustness and stability to the AT (woA) network during a connectivity loss of 0.8 after the directed attack in degree (Fig. 4F) and, in addition, to CT (woA) network in betweenness (Fig. 4G). This behavior highlights the important role of Acinetobacter within the community conditioned by its high abundance and possible potential for antibiotic resistance.

Fig. 4

Robustness comparison after node removal within groups (CT and AT). Connectivity loss measured after directed attack (degree and betweenness) in Coxiella’s presence (wC) and removal (woC): (A) CT (wC vs. woC) and (B) AT (wC vs. woC) in degree, (C) CT (wC vs. woC), and (D) AT (wC vs. woC) in betweenness. Connectivity loss measured after directed attack (degree and betweenness) in Acinetobacter’s presence (wA) and removal (woA) between the same group: (E) CT (wA vs. woA), and (F) AT (wA vs. woA) in degree, (G) CT (wA vs. woA), and (H) AT (wA vs. woA) in betweenness

During the addition of nodes, the LCC values was higher in the CT network following Coxiella removal (Fig. 5A). In contrast, an overlap was observed between the AT and AT (woC) networks (Fig. 5B), suggesting that the absence of the taxon did not influence the response to the addition of nodes. On the other hand, the presence of Acinetobacter conferred stability in both CT and AT conditions (Fig. 5C-D). The robustness conferred by Acinetobacter may be due to its potential for resistance to antibiotic treatment and therefore its possible action in maintaining stability. Interestingly, only the LCC value measured in the AT network after the addition of 1000 nodes demonstrated greater robustness in the absence of the taxon (Fig. 5D). The overlap and increase in APL values within the same condition for the CT and AT groups in the presence and after the individual removal of each taxon continue to suggest loss of functional connectivity and increase in competition within the network with the addition of nodes (Fig. 5E-H).

Fig. 5

Robustness comparison against addition of nodes after taxa’s removal between the same groups (CT and AT). LCC values are represented and compared in Coxiella’s presence and removal between the same groups: (A) CT (wC vs. woC), and (B) AT (wC vs. woC). LCC values are represented and compared in Acinetobacter’s presence and removal between the same groups: (C) CT (wA vs. woA), and (D) AT (wA vs. woA). APL values are represented and compared in Coxiella’s presence and removal between the same groups: (E) CT (wC vs. woC), and (F) AT (wC vs. woC). APL values are represented and compared in Acinetobacter’s presence and removal between the same groups: (G) CT (wA vs. woA), and (H) AT (wA vs. woA)

Changes in predicted functional profiles in response to antibiotic treatment

The alterations in microbial community composition and structure were assessed to determine whether they affected the inferred functional profile of tick microbiota. A thorough analysis was conducted, comparing the composition, diversity, and relative abundance of metabolic pathways in the microbiota of H. longicornis ticks from both the AT and CT groups. The investigation revealed a greater richness in functional profiles within the AT group compared to the CT group. Specifically, diversity metrics such as observed features (Fig. 6A) and Pielou’s evenness index (Fig. 6B) were higher for the AT group compared to the CT group. Differences were also observed in the relative abundance of pathways between the functional profiles of AT and CT groups (Fig. 6C, Supplementary Table S6). Analysis revealed both unique and shared predicted metabolic pathways within the microbiota of both groups. Specifically, five pathways were identified as unique to the CT tick microbiota, while the AT tick microbiota exhibited eleven unique pathways within which the 2-heptyl-3-hydroxy-4(1H)-quinolone biosynthesis pathway is found (Fig. 6D, Supplementary Table S7 and 8). Furthermore, both microbiota shared 398 pathways (Fig. 6D, Supplementary Table S7 and 8), predominantly associated with biosynthesis, as documented on MetaCyc database [51].

Fig. 6

Predicted functional profile analysis. Comparative analysis of predicted functional profiles between CT (green) and AT groups (red). A Observed features, and (B) Pielou’s evenness index were employed to assess pathway alpha diversity. C Differential abundance of predicted metabolic pathways. Volcano plot showing the differential abundance of pathways in CT and AT microbiota. The pathways with significant differences (Wald test, p < 0.05) between the groups are represented by pink dots. D Venn diagram illustrates shared and unique predicted bacterial pathways between the CT and AT microbiota

Discussion

Our findings provide support for the hypothesis that antibiotic treatment (AT) disrupts the microbial community assembly and network interactions within H. longicornis ticks. This disruption affects key taxa such as Coxiella and Acinetobacter, and complement those of Wei et al. [16] in understanding how antibiotics disrupt colonization resistance within tick microbiota, thereby facilitating the transstadial transmission of B. microti. Colonization resistance is a critical ecological function provided by the native microbiota as it prevents the establishment and proliferation of pathogens within the host [55]. This resistance is mediated by several mechanisms, including competition for resources [26, 27] and modulation of the host's immune response [56, 57]. The network-based traits associated with AT and their influence on reduced colonization resistance and enhanced Babesia transstadial transmission can be summarized: 1) the AT group had more connected nodes but less interactive microbial community compared to the CT group. 2) Despite having more connections, the AT group had fewer correlations and greater modularity compared to the CT group. 3) There were significant differences in network centrality measures between AT and CT groups and in the connectivity and network roles of key taxa such as Coxiella and Acinetobacter. 4) The robustness tests (i.e., node removal strategies) demonstrated that the AT networks were less stable against specific attacks compared to CT networks. And 5) The AT group displayed a greater richness in functional profiles.

Enhanced network connectivity and modularity in AT microbiota

Antibiotic treatment modifies the community assembly as a consequence of the reduction of microorganisms due to its bacteriostatic or bactericidal effect [58, 59]. Consequently, interactions such as co-occurrence (positive correlation) or co-exclusion (negative correlation) may change [60, 61]. We found topological variations indicating distinct community assembly between the CT and AT groups, with a decrease in the number of interactions between community members. However, no considerable changes were evident in the nature of the interactions, with a balance between cooperation and co-exclusion existing in both groups. This change in topology may reflect a shift towards a microbial community where fewer dominant interactions exist, potentially reducing competitive exclusion and colonization resistance. The interactions established between microorganisms within communities can alter the susceptibility of their members to antibiotic treatment [57]. The nature and number of interactions established can determine whether the resistance of a microbial species to a certain antibiotic confers a protective effect [57, 62].

Altered centrality measures

Analyzing antibiotic effects on node centrality and network strength in Babesia-infected ticks' co-occurrence networks revealed significant insights into microbial interaction dynamics under perturbations. Specifically, differences in local network centrality measures (degree, betweenness centrality, closeness centrality, and eigenvector centrality) between CT and AT groups highlights distinct microbial network patterns. Interestingly, the AT group showed greater adaptability to disruptions with additional nodes, while both groups exhibit similar overall robustness. Furthermore, taxa such as Coxiella and Acinetobacter showed disrupted connectivity in the AT group. Coxiella exhibited reduced interactions, whereas Acinetobacter experiences a loss of connected nodes, underscoring their critical roles in maintaining microbial network stability. This provides insights by demonstrating differences in network centrality measures and the roles of key taxa like Coxiella and Acinetobacter between both groups which deepen our understanding of how antibiotic treatment impacts tick microbiota dynamics.

Impact of key taxa including Coxiella and Acinetobacter

In their initial study, Wei et al. [16] observed that antibiotic treatment led to a significant alteration in the diversity and abundance of the tick gut microbiota. They noted shifts towards genera like Acinetobacter, known to thrive under antibiotic pressure [16], and a decreased in Coxiella, which typically acts as an endosymbiont in ticks [63], providing protection against tick-borne pathogens (TBPs) [55, 64]. This alteration suggests a disruption of the balanced microbial ecosystem that typically supports colonization resistance to Babesia infection. Our study further demonstrated that these shifts are not just alterations in abundance but also critical changes in the connectivity and network roles of key taxa such as Coxiella and Acinetobacter. The reduced interaction within the microbial community and the loss of connected nodes indicates a weakened network structure, reducing the microbiota’s ability to collectively resist new colonization by pathogens, in this case, Babesia species. Symbiont-mediated protection, predicted through ecological modelling, is increasingly observed in natural insect populations, serving as a potent mechanism to maintain symbiont prevalence [55]. Ticks harbor various nutritional symbionts, including Coxiella-like endosymbionts (CLE) [65], Francisella-like endosymbionts (FLE) [66] and Rickettsia [67]. Additionally, ticks can occasionally harbor other symbionts, such as Wolbachia [68] and Arsenophonus [69], which are less commonly found in ticks. Experimental evidence demonstrates vertically transmitted symbiont’s protective role against pathogens or predators. Notably, studies show that antibiotics targeting Coxiella-LE and Rickettsia-LE reduce their densities [70], potentially affecting host fitness. This is exemplified in H. longicornis, where tetracycline treatment reduced Coxiella-LE levels, affecting the tick's fitness [70]. Additionally, the obligate Coxiella-like symbiont has been shown to manipulate the reproduction of its host, Amblyomma americanum [71].

Despite varying assembly patterns, Acinetobacter contributed to the stability of the AT group following node removal. Recently, this genus has been identified in the microbiome of field collected H. longicornis larvae, nymphs, and adults, with a higher abundance in larvae [72]. Additionally, it was detected in the midgut microbiome of fed specimens collected in China, demonstrating that Acinetobacter can stably colonize the midgut of H. longicornis [73]. Known for its broad antibiotic resistance [74, 75], this genus monitors antibiotic resistance genes and underscores ticks' potential as pathogen reservoirs [70]. For instance, the pathogen Acinetobacter baumannii shows high antimicrobial resistance development and acquisition of new resistance determinants [70, 76, 77]. Its response to antibiotics under iron limitation and oxidative stress is of interest, given its ability to express resistance to a wide range of antibiotics used in human medicine [78, 79].

Furthermore, the analysis of keystone taxa identification in both the CT and AT groups revealed distinct compositions between the two groups. Ammoniphilus, Noviherbaspirillum, and others were prominent in the CT group, whereas Methyloceanibacter, Lysobacter, and others were notable in the AT group. This suggests a shift in the microbial ecosystem following antibiotic intervention. These findings are consistent with prior research highlighting the significant role of specific bacterial taxa in tick microbiota and their potential susceptibility to antibiotics [18] [46]. These disruptions in microbial community assembly and node centrality likely contribute to the reduced colonization resistance observed in the AT group, as the antibiotics compromise the integrity and stability of the microbial network, making it more susceptible to colonization by pathogens like Babesia. Integrating community structure with functional dynamics represents a fundamental pursuit in microbial ecology [80], necessitating the exploration of microbial co-occurrence patterns and the identification of keystone taxa essential for ecosystem processes. The intricate interplay elucidated in this study among antibiotic treatment, keystone taxa, and microbial community dynamics provides profound insights beyond diversity studies. These findings establish a scientific basis for potential interventions targeting the modulation of tick microbiota to mitigate tick-borne diseases [81].

Robustness against network perturbations

The concept of robustness, which refers to the resistance of a network, can be elucidated through percolation theory [82], offering insights into information flow among network nodes [30]. Prior studies have employed in silico node removal to evaluate microorganism influence on plant microbiota properties [83]. This approach validated as a tool for predicting ecosystem behavior [30], suggests network robustness as indicative of microbial community resilience in diverse animal species, from arthropods [84, 46] to mammals [85]. This approach aims to induce infection-refractory states in ticks or other vectors [86, 87], ultimately reducing or blocking vector-borne pathogen transmission [88]. In this study, percolation theory was applied to evaluate network robustness by assessing loss in connectivity through degree, cascading, betweenness, and random attacks [30]. The results showed node removal had a more pronounced effect on both CT and AT networks compared to the random method, especially in the betweenness method. These findings suggest that while antibiotics administration can be discerned in the robustness test, its impact on bacterial community assembly appears limited compared to other disturbance factors [84]. This observation is consistent with previous research demonstrating a significant reduction in network robustness following directed attacks in the microbiota of mice exposed to antibiotics and fed a high-fat diet, as compared to untreated mice on the same diet [85]. From the result of this study, it can be proposed that specific antibiotics may target keystone taxa, potentially leading to ecosystem collapse [85]. The mechanisms by which antibiotics modulate the taxonomic and functional profiles of the tick microbiota are yet to be determined.

Functional implications

The analysis of predicted metabolic profiling of the microbiota was incorporated based on the novel bioinformatics tool PICRUSt2 [49]. Researching the effects of antibiotics on the functional characteristics of the microbiota of H. longicornis ticks has yielded some interesting results. In addition to changes in assembly patterns and community robustness, the functional profiles of microorganisms are affected as a consequence of antibiotic treatment. The AT had greater functional richness and evenness, implying broader metabolic capabilities compared to the CT group. This suggests that antibiotics may influence the metabolic diversity of the microbiota. Despite the existence of shared metabolic pathways, unique metabolic pathways were evident in both groups, denoting metabolic differences. Changes in the physicochemical properties of the ecosystem and the biological products that certain species produce may contribute to antibiotic resistance in other species [57]. Among the unique metabolic pathways in the AT group, the 2-heptyl-3-hydroxy-4(1H)-quinolone biosynthesis pathway was found. This molecule, known as the Pseudomonas quinolone signal (PQS), is important for its role in the quorum-sensing system that regulates biofilm formation, secondary metabolite production, pigment and virulence factor production, motility, and membrane vesicle formation [89]. Previous studies have revealed that exposure to subinhibitory concentrations of antibiotics induces the production of PQS in Pseudomonas aeruginosa [90]. Additionally, PQS has been shown to influence not only P. aeruginosa populations but also other bacterial species, regulating microbial communities within a specific ecosystem [91]. Therefore, in our study, antibiotic treatment may be triggering these cell-to-cell mechanisms to adapt and survive, facilitating the overgrowth of Pseudomonas species and other antibiotic-resistant bacteria. This could reduce the growth of beneficial endosymbionts species in the tick microbiome, which play an important role in tick metabolism and reproductive fitness [71, 92]. A disrupted microbiome could reduce tick’s resistance to colonization of pathogens and other antibiotic-resistant bacteria, increasing the risk of transmitting diseases to the host. Furthermore, the proliferation of these bacteria species can potentially turn ticks into reservoirs and spreaders of antibiotic-resistant genes [93]. This underscores the importance of understanding the role of commensal organisms, such as Coxiella-like endosymbionts and Acinetobacter as significant metagenomic biomarkers [94]. For instance, the clinical significance of variances in antimicrobial susceptibility profiles among distinct genomic clusters of Acinetobacter has been elucidated [95], while the escalating recognition of efflux systems in facilitating multidrug resistance in Acinetobacter further underscores their critical role [96, 97].

Study limitations

This study has four major limitations: 1) Working with a small sample size may not capture the full spectrum of microbial diversity in both the AT and CT groups. Rare or less abundant microbial taxa might be underrepresented or missed, leading to an incomplete picture of the microbiome [98] and making the findings less robust for drawing definitive conclusions about the diversity and composition of the tick microbiome between both groups [99]. 2) Bacterial community structures vary significantly across tick life stages [100]. In some cases, the interaction pattern (co-occurrence or co-exclusion) changed according to the development stage [101]. For instance, higher bacterial diversity was observed in nymphal stages compared to adult stages of Ixodes ricinus ticks, likely due to distinct host-selection behaviors between immature and mature ticks [100]. These variations can significantly impact how antibiotics affect the tick microbiota. 3) Our study utilized metagenomic data to characterize the functional profiles of the tick microbiota, focusing on potential genes and pathways present, but did not directly assess their active expression or functional significance under antibiotic exposure. Using richness and evenness metrics provide valuable insights into the diversity and distribution of functional pathways but may not fully capture functional redundancy, where different microbial species perform similar roles. This redundancy can mask changes in community composition without altering richness or evenness. This contrasts with meta-transcriptomics, which can provide direct insights into actively expressed genes by studying transcriptional regulation, metabolite dynamics, and protein signaling within tick microbial communities under antibiotic treatment. 4) Antibiotics targeting bacterial communities can disrupt microbial diversity and abundance, impacting microbiota-host interactions and potentially compromising the tick's ability to control pathogens [102]. These changes can also influence ecological niches, community stability, and interactions within the environment. Previous research has highlighted specific tick genes like longicin which is defensin-like protein of H. longicornis exerting anti-microbial and anti-fungal activity [103] and TROSPA, serum amyloid A, and calreticulin, which are implicated in vector-pathogen interactions during Babesia infection [104]. Further investigation into these biomarkers and proteins would enhance our understanding of how tick microbiota responds to antibiotics during Babesia infection, providing a better understanding of the tick immune response under antibiotics exposure and bring insights for optimizing treatment strategies against tick-borne diseases.

Conclusions

The current study provided a comprehensive exploration of the H. longicornis microbiota's response to antibiotic treatment, particularly in the context of B. microti transstadial transmission. The investigation into co-occurrence networks and keystone taxa showed the vulnerability of the tick microbiota to antibiotic-induced perturbations, observing notable differences in network centrality measures (degree, betweenness, closeness, and eigenvector centrality) and distinct keystone taxa compositions. Specifically, antibiotic treatment altered the clustering pattern of key taxa such as Coxiella and Acinetobacter, leading to a less robust microbial network and increased susceptibility to disturbance. Remarkably, the study uncovered novel insights into the functional consequences of antibiotic treatment, revealing increased functional richness and evenness in the antibiotic-treated group, implying broader metabolic capabilities and potential shifts in the importance of specific functions. These results suggest that antibiotic treatment can disrupt the microbial balance within ticks, decreasing their resistance to pathogen colonization. These findings carry broader ecological implications, emphasizing the need to consider functional aspects in understanding antibiotic-mediated reduction of colonization resistance and its implications for Babesia transstadial transmission. Our results offer a more detailed comprehension of tick microbiota dynamics under antibiotic treatment, revealing insights that were not evident through traditional assessments of bacterial diversity or abundance alone.

The intricate interplay between antibiotic treatment, microbial community dynamics, and functional profiles underscore the complexity of the tick microbiota, offering avenues for further research to manipulate these microbial communities for effective control of tick-borne diseases. In this context, anti-microbiota vaccines have been designed to modulate the tick microbiome by targeting essential microbial taxa [46, 105]. The use of co-occurrence networks and centrality measures to identify key taxa provides a robust method for developing these vaccines [105]. Anti-microbiota vaccines have been shown to effectively modulate the tick microbiome, impacting tick performance and pathogen colonization, thus supporting the development of this strategy for controlling tick-borne pathogens [106, 46, 107]. Additionally, alternative strategies in other vectors, such as paratransgenesis and phage therapy, have been explored to target specific taxa, which have demonstrated success in reducing pathogen load and vector competence [108, 109].