Resistance of tick microbiome to biological disturbance

doi:10.21203/rs.2.17733/v1

Download PDF

Research

Resistance of tick microbiome to biological disturbance

https://doi.org/10.21203/rs.2.17733/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background : Ixodes scapularis ticks harbor microbial communities including pathogenic and non-pathogenic microbes. Pathogen infection increases the expression of several tick gut proteins which disturb the tick gut microbiota and impact bacterial biofilm formation. Anaplasma phagocytophilum induces ticks to express I. scapularis IAFGP, a protein with antimicrobial activity while Borrelia burgdorferi induces the expression of PIXR. Here, we tested the resistance of I. scapularis microbiome to A. phagocytophilum infection, antimicrobial peptide IAFGP, and anti-tick immunity specific to PIXR. Results : We demonstrate that A. phagocytophilum infection and IAFGP affect the taxonomic composition and taxa co-occurrence networks but had no effect on the functional traits of tick microbiome. In contrast, anti-tick immunity disturbed the taxonomic composition and the functional profile of tick microbiome, by increasing both taxonomic and pathways diversity. Mechanistically, we show that anti-tick immunity increases the representation and importance of polysaccharide biosynthesis pathways involved in biofilm formation while these pathways are under-represented in the microbiome of ticks infected by A. phagocytophilum or exposed to IAFGP. Conclusions : These analyses revealed that tick microbiota is highly sensitive to anti-tick immunity, while it is less sensitive to pathogen infection and antimicrobial peptides. Results suggest that biofilm formation is a defensive response of tick microbiome to anti-tick immunity.

General Microbiology

Ixodes scapularis

microbiota

disturbance

functional redundancy

core microbiome

Ixodes scapularis is an important vector of Borrelia burgdorferi and Anaplasma phagocytophilum, causal agents of Lyme disease and human granulocytic anaplasmosis, respectively (1). Borrelia burgdorferi is an extracellular bacterium and A. phagocytophilum is an obligate intracellular bacterium. However, both bacteria are maintained in a mammal-tick infectious cycle that involves intimate interactions with ticks during the colonization of its gut and subsequent migration to the salivary gland before transmission to the vertebrate host (1). Ixodes scapularis harbors a diverse group of native microbes ranging from viruses to bacteria (2). These microbial communities influence the ability of B. burgdorferi and A. phagocytophilum to colonize and persist within the vector (3–5). Microbial biofilms (6), generated by diverse bacterial species, are essential for bacterial colonization and successful symbiotic relationships in arthropod hosts (7–9). Borrelia burgdorferi and A. phagocytophilum colonization in I. scapularis is associated with changes in the taxonomic composition of the microbiota resulting in diffuse biofilms (3,4). Tick gene expression is affected by B. burgdorferi and A. phagocytophilum infection in I. scapularis gut. A secreted protein with a Reeler domain (PIXR) and an antifreeze glycoprotein (IAFGP) with antimicrobial properties are upregulated by B. burgdorferi and A. phagocytophilum, respectively (3,4,6). The IAFGP and the anti-PIXR immunity also alter the tick microbiota and result in diffuse and dense biofilms, respectively (3,4).

Based on these evidences, we asked to what degree exposure to pathogens, antimicrobial peptides or anti-tick immunity would impact the resistance of tick microbiome. We hypothesized that these disturbing factors induce taxonomic and functional changes that impact the resistance of tick microbiota and may explain the properties of biofilms α and β (Figure 1). To test this hypothesis, published 16S rRNA gene amplicon sequence datasets obtained from tick microbiota upon disturbance by A. phagocytophilum (3), antimicrobial peptide IAFGP (3) and anti-tick immunity specific to PIXR (4) were used for annotation of taxonomic profiles and prediction of functional traits of the microbiome using the state-of-the-art metagenomics tool PICRUSt2 (7). The taxonomic and functional profiles were quantified and compared in terms of abundance, diversity and composition, parameters used to measure the resistance of microbial communities to disturbance (8,9). The taxonomic and functional structure of the microbial communities was then quantified and compared using co-occurrence networks. Resistance to taxa extinction and connectedness loss was tested on the taxonomic networks by systematic removal of taxa in random and directed manners. Large effect on resistance was considered when all taxonomic and functional parameters of abundance, diversity, composition and resistance to taxa extinction were affected by disturbance. Here we show that anti-tick immunity has a large impact on the tick microbiome, enhancing both taxonomic and functional diversity. Microbial community response to anti-tick immunity produced the over-representation of pathways involved in biofilm formation. In contrast, tick microbiota is less sensitive to A. phagocytophilum and antimicrobial peptide, since they alter the taxonomic composition but not the pathway profile of the microbiome. We conclude that anti-tick immunity is a major disruptor of tick microbiota favoring the abundance of strong biofilm formers. Results demonstrate that communities of bacteria are functionally redundant, suggesting a mechanism by ticks selecting the adequate microbiome fulfilling a core set of functions.

Effect of biological disturbance on taxonomic and functional profiles of tick microbiome

Tick microbiota showed different taxonomic composition in the six 16S datasets (Additional File 1 Figure S1a). Major differences were found between the tick stages (i.e. nymph and larvae). A Venn diagram analysis shows that each experimental group was associated to a set of specific bacterial genera that was not shared by the others (Figure 2a). A small taxonomic core of 61 bacterial genera (total: 821, 7.4%) was shared by at least one individual tick of each of the experimental groups (Figure 2a). The presence of a highly reduced taxonomic core was confirmed by taxa counting after progressive sampling, since 20% and 100% of all samples (n=98) shared only 140 (total: 821, 17.1%) and 10 taxa (1.2%), respectively (Figure 2b). The only 10 taxa present in 100% of the samples were included in Additional File 2 Table S1. The existence of a reduced taxonomic core was confirmed by fuzzy logic analysis which showed that less than 50 bacterial taxa have a probability of membership higher than 50% (Figure 2c).

Community diversity indexes showed that anti-tick immunity produced a significant increase in phylogenetic diversity (Kruskal-Wallis test, p <0.001) and species evenness (Kruskal-Wallis test, p <0.001) (Additional File 1 Figure S2a). However, A. phagocytophilum infection or the antimicrobial peptide did not modify the alpha diversity indexes (Additional File 1 Figure S2b,c). Disturbing factors also changed the abundance of several taxa (Additional File 1 Figure S3). Pairwise comparisons between datasets within the same experiment showed that anti-tick immunity, A. phagocytophilum infection and the antimicrobial peptide changed dramatically the abundance of specific bacterial genera (Additional File 1 Figure S4). We observed that the abundance of 63 (total 78, 80.8%), 18 (total 46, 39.0%) and 22 (total 51, 43.1%) bacteria genera increased in response to anti-tick immunity, A. phagocytophilum infection, and antimicrobial peptide, respectively. These results suggest that biofilm α (Figure 1) may be associated with an increase in the abundance of biofilm formers, while biofilms β1 and β2 (Figure 1) may be associated with a reduction in the number of biofilm formers. In agreement with this idea, the bacterial genera with higher abundance in response to anti-tick immunity include the strong biofilm formers Mycobacterium (33), Tepidimonas (34), Rothia (35) and Leuconostoc (36), whereas A. phagocytophilum infection and antimicrobial peptide reduced the presence of the biofilm formers Gracilibacteria (37) and Enterococcus (3), respectively.

Based on the changes observed in the taxonomic composition and taxa abundance after disturbance, we expected a similar trend in the functional traits encoded in the tick microbiomes. Surprisingly, the Venn diagram analysis showed that all experimental groups shared the majority (i.e. 381) of the metabolic pathways (total: 437, 87.2%) and almost none of the disturbing factors was associated to a set of specific metabolic pathways (Figure 3a). Progressive sampling of metabolic pathways (Figure 3b) and fuzzy logic analysis (Figure 3c) showed the persistence of a broad functional core despite disturbance. Metabolic pathways included in the functional core across all the samples are available in Additional File 3 Table S2.

We next asked whether this functional redundancy would be observed at the bacterial community level. To evaluate this idea, we built networks of co-occurring bacteria (described in detail in the next section) to test the functional profiles of bacterial communities in tick microbiota. Network modules are formed by bacteria that co-occur more often with each other than with the other bacteria in the network and they can be considered as bacterial communities (38). The three largest modules (>75% of nodes) were selected in each network. Subsequently, we predicted and compared their pathway profiles. The results show that the network-derived modules of co-occurring bacteria in each network are similar in pathway content and abundance (Additional File 1 Figure S5). These findings suggest that bacterial communities of tick microbiota are functional units containing a highly redundant set of metabolic pathways.

Disturbance reduces the tolerance of co-ocurrence networks to taxa extinction

Bacteria co-occurrence networks were used to quantify changes in topology and resistance to taxa extinction. Visual inspection of the taxonomic networks revealed that each disturbing factor changed the network topology (Figure 4) a result supported by their numerical features (Table 1). The topological changes observed in the co-occurrence networks of tick microbiota after disturbance indicate that anti-tick immunity, pathogen infection and antimicrobial peptide reshape the co-occurrence of microbial communities in the tick gut. Disturbance affects the connections between modules, explained by the average clustering coefficient. The anti-tick immunity induced modules with higher clustering while A. phagocytophilum and the antimicrobial peptide reduced dramatically this property. Anti-tick immunity produced a network with two large modules and higher number of nodes compared with the control (Additional File 1 Figure S6). Anaplasma phagocytophilum completely altered the network by increasing the number of modules and decreasing the connection among them (Additional File 1 Figure S7), properties that were similar in the antimicrobial peptide network (Additional File 1 Figure S8). All the networks of disturbed microbiota had higher diameters (increased laxity) than the controls.

To test the effect of disturbing factors on network resistance, the networks of co-occurring taxa were subjected to two different attack strategies: (i) random removal of nodes and (ii) directed removal of nodes starting from those with the highest centrality. In each case, we assessed the loss of connections and secondary extinctions in the network. Bacterial networks resulting from ticks exposed to anti-tick immunity (Additional File 1 Figure S9a), A. phagocytophilum (Additional File 1 Figure S9b) and anti-microbial peptide (Additional File 1 Figure S9c) were less tolerant to both strategies of taxa removal. Every disturbance had an effect on the connectivity of the networks, reaching a 50% of disconnected taxa when 12%, 20%, or 22% of taxa were removed for anti-tick immunity, A. phagocytophilum infection and anti-microbial peptide, respectively (against 30%, 24% and 23% removal in control groups, respectively). The anti-tick immunity induced a larger disconnection of the network than other disturbances to tick microbiota.

Anti-tick immunity increases the representation of biofilm formation pathways in tick microbiome

We then measured the impact of disturbance on the metabolic traits of tick microbiome. We first compared the alpha and beta diversities of metabolic pathways found in the microbiomes. Anti-tick immunity induced an increase in the pathway richness (Kruskal-Wallis test, p <0.001) and evenness (Kruskal-Wallis test, p <0.001) whereas A. phagocytophilum and the antimicrobial peptide had no significant effect on these alpha-diversity metrics (Additional File 1 Figure S10). In agreement with this, the principal coordinate analysis (PCoA) of beta diversity showed that only anti-tick immunity produced a significant differentiation (PERMANOVA test, p <0.001) in the functional profiles of tick microbiome (Additional File 1 Figure S10). Anaplasma phagocytophilum and the antimicrobial peptide produced minor and no change, respectively, on the functional diversity of tick microbiome (Additional File 1 Figure S10).

Significant differences in the relative abundance of the metabolic pathways in response to the disturbing factors were observed only in ticks exposed to anti-tick immunity (Figure 5a). Despite the changes in the abundance of bacterial genera in response to A. phagocytophilum and antimicrobial peptide these disturbing factors produced minor (Figure 5b) and no change (Figure 5c), respectively, on the abundance of metabolic pathways of the tick microbiome. A detailed comparison of metabolic pathways abundance across all the samples (Gneiss test) resulted in a dendrogram heatmap showing that five clusters of functional categories explained the major differences among groups (Additional File 1 Figure S11). The metabolic pathway composition of each cluster is available in Additional File 4 Table S3. To evaluate the relative importance of each of these pathways, we built networks using pathway co-occurrence and used the WD as a proxy of pathway importance. Functional network of tick microbiome exposed to anti-tick immunity showed higher connectivity than those of A. phagocytophilum-infected and antimicrobial peptide-treated ticks (Table 2). Disturbance of tick microbiome changed the weighted degree (WD) of metabolic pathways in the functional networksAmong the pathways with at least 2fold change (log2 =1) in WD, 87 (total 105, 82.9%) were over-represented in response to anti-tick immunity, while among the pathways with significant changes in response to A. phagocytophilum and antimicrobial peptide, 79 (total 87, 90.8%) and 74 (total 93, 79.6%) were under-represented, respectively.

The pathways with the greatest differences between groups (log2 ≥/≤ 2) are shown in Additional File 1 Figure S12. Among the pathways with WD log2 ≥/≤ 2 in response to anti-tick immunity, three were related with polysaccharide biosynthesis, an essential process in biofilm formation. These pathways were colanic acid biosynthesis (39), peptidoglycan biosynthesis (40), and O-antigen biosynthesis (41). With few exceptions, all the molecular components of these polysaccharide biosynthesis pathways were identified (Additional File 1 Figure S13a). The relative abundance of these pathways was also found to increase significantly in response to anti-tick immunity (Additional File 1 Figure S13b). Other pathways also involved in biofilm formation, including polysaccharide degradation (39) and siderophore biosynthesis (42,43), were also identified. Based on the comparison of WD values for the three group of pathways (i.e. polysaccharide biosynthesis, polysaccharide degradation, and siderophore biosynthesis), a model of the functional profiles of biofilms α, β1 and β2 was proposed (Figure 6).

In this study we tested the resistance of I. scapularis microbiome to three disturbing factors, namely anti-tick immunity, pathogen infection and an antimicrobial peptide. To this end, we used 16S gene sequences available from previous publications (3,4). These studies were originally aimed to elucidate the effects of anti-tick immunity specific to PIXR, A. phagocytophilum infection, and antimicrobial peptide P1, on tick gut microbiota composition and bacterial biofilm formation (3,4). Using recent advances on bioinformatics tools and benchmarks (i.e. Dada2 (13), and the classify‐sklearn naïve Bayes taxonomy classifier (19)), our taxonomic analysis pipeline improved the accuracy in taxonomic classification of the sequences. This methodology allowed us to resolve each amplicon sequence variants (ASVs) instead of operational taxonomic units (OTUs). In addition, we incorporated the analysis of predicted metabolic profiling of the microbiomes, based on the novel bioinformatics tool PICRUSt2 (7). Resistance could be defined as the degree to which the community withstands change in the face of disturbance (8). To assess resistance of tick microbiota, we measured to what extent the taxonomic or functional profiles of the tick microbiome remained unchanged under different disturbing factors (9). Furthermore, we explored the structure of microbiota by applying network analysis, a well-established methodology for detecting interactions within the microbiota (44,45). A further test of resistance was performed on the taxonomic co-occurrence networks by quantifying the loss in connectivity and secondary extinctions due to taxa removal.

The complex community of microbes living in the tick guts (2) plays an important role in pathogen colonization (3,4) and potentially on tick fitness. However, despite its expected importance defining the composition of tick microbiota has remained an elusive task. In this regard, a central question is whether a ‘core’ microbiota consisting of bacterial groups common to all ticks exists (46). The existence of a core microbiota has been previously addressed in terms of taxonomic composition (47,48). Here we showed that the taxonomic core of I. scapularis is highly reduced whereas the functional core is an important element of tick microbiome. We found evidence of a functional core in I. scapularis defined as 300 pathways present in 100% of the samples analyzed. However, I. scapularis core microbiome was not defined by taxonomic associations as only 10 taxa were shared by 100% of the samples. This suggests that microbial communities in ticks shared a set of metabolic pathways, regardless of the taxonomic identity of the microbes in tick microbiota. The existence of a functional core instead of a taxonomic core underlay the existence of functional redundancy (when multiple taxa contribute with the same metabolic function) of tick microbiome.

This study demonstrated that functional redundancy is an important property of tick microbiome, since we also demonstrated that a core functional microbiome occurs at the level of bacterial communities within the microbiome. Selected network modules can be considered as communities of bacteria with the highest frequency of co-ocurrence (29). The presence of different microbial communities in ticks suggest the existence of microbe–microbe interactions shaping the composition of tick microbiota. Little is known, however, about the driving forces behind these interactions and their effects on the dynamics of microbial communities. One hypothesis suggests that tick genotype and ecological factors influence tick symbionts (49) and/or strongly interconnected taxa to recruit other microbes to form a holobiont, the assemblage of different species (e.g. host and bacteria) that form an ecological unit (50). Our results showed that network modules had similar metabolic traits which suggest that different bacterial communities within ticks shared a functional core. A possible implication of such property of tick microbiota is to keep a functional stability independently of taxonomic changes caused by disturbance. We propose that the functional redundancy that emanates from the modular structure ensures that key functions are maintained and reinforce the criterion that a functional core microbiome is resilient to biological disturbance.

Ixodes scapularis microbiota was altered by disturbing factors including anti-tick immunity, pathogen infection and antimicrobial peptide. Different mechanisms thus accounted for tick microbiota modulation (3,4). The rewiring of taxonomic networks after disturbance revealed the sensibility of tick microbiota to disturbance. Moreover, the taxonomic networks resulting after disturbance were less resistant to taxa removal. However, the only factor that had a large effect on bacterial diversity was anti-tick immunity. In accordance, the functional microbiome was affected only by anti-tick immunity, while the A. phagocytophilum infection and microbial peptide did not altered the pathway diversity of tick microbiomes. This is consistent with the idea that the tick microbiome evolved through the recruitment of microbial populations of different taxa, but sharing a set of functions with possible impact on tick fitness (51,52).

Immunity to PIXR induced a significant increase in microbial diversity and the abundance of biofilm formers which was associated with dense biofilms α (Figure 1a, (4)). However, the resistance test revealed that the microbial communities of ticks fed on animals immunized with PIXR were less resistant to taxa removal, suggesting that these microbial communities were less organized and connected. This finding match with the assumption that the proportion of interactions within a biological system modulates its stability, but it may become suddenly unstable as the system becomes larger (53,54). This also suggests that despite the functional core may provide essential metabolic pathways, all microbes are not interchangeable in the tick microbial community. This is in agreement with a recent study showing phylosymbiosis between Ixodes ticks and their microbial communities (52). Anti-tick immunity also increased the functional diversity and the importance of pathways involved in biofilm formation including colanic acid biosynthesis, a major component of bacterial biofilms (55). Biofilm creates a favorable environment that increases antibiotic resistance, impairs host immunity, and increases tolerance to nutritional deprivation (56–58). The formation of biofilms α (Figure 1a) due to increase in biofilm formers may be part of a protective response of tick microbiota to anti-tick immunity. This raises an interesting question; can anti-tick vaccination trigger the formation of biofilms α that increase the resistance of ticks to vaccines? It is noteworthy to mention that PIXR is a negative regulator of biofilm’s formation (4). Therefore, future studies should test whether immunity against tick proteins not related with biofilm formation would also induce the formation of biofilms α.

Anaplasma phagocytophilum and antimicrobial peptide were associated with diffuse biofilms (Figure 1b,c, (3)) and produced taxa replacement without impact on the microbial diversity metrics. However, A. phagocytophilum and antimicrobial peptide communities resulted in networks with low resistance to taxa removal. The infection with A. phagocytophilum and the treatment with the antimicrobial peptide did not cause significant modifications on the diversity and abundance of the functional profiles of tick microbiome whereas the relative importance (measured as changes in WD) of several pathways changed in response these disturbing factors. A common property of the microbiome disturbed by A. phagocytophilum and the antimicrobial peptide was a decrease in the importance of polysaccharide biosynthesis pathways (Figure 6b,c) which are important in the formation of bacterial biofilms (55,59,60). Our results are consistent with previous results (3) reporting a decrease in the exopolysaccharide PNAG in tick gut infected with A. phagocytophilum. The decrease of biofilms upon A. phagocytophilum colonization was suggested to be caused by the binding of IAFGP to D-alanine which blocks the formation of gram positive biofilms (3). Our study suggests that the underrepresentation of biofilms forming pathways may be an additional mechanism by which A. phagocytophilum decreases tick guts biofilms.

We conclude that the functional core is an important component of the resistance of tick microbiome to disturbance. The functional core is redundant and is not defined or restricted by taxonomic composition. This is confirmed by the existence of a rich and redundant functional core in opposition to a highly reduced taxonomic core. We strongly suggest to turn the classic taxonomic approach of comparing the tick microbiome between species into a functional framework. Despite functional redundancy, there were changes in abundance and importance of pathways that influenced the functional properties of the core microbiome in response to disturbance. Biofilm formation pathways are part of the response of tick microbiota to disturbance. Tick microbiota is highly susceptible to anti-tick immunity and the increase of biofilm formation pathways may be a response of tick microbiota to this disturbing factor. In contrast, tick microbiota is less sensitive to A. phagocytophilum colonization and tick antimicrobial peptide which may reflect the coevolution and adaptation between tick-borne pathogens, the microbiota and the vector.

This study demonstrated that anti-tick immunity has a large impact on the tick microbiome, enhancing both taxonomic and functional diversity. Microbial community response to anti-tick immunity produced the over-representation of pathways involved in biofilm formation. In contrast, tick microbiota is less sensitive to A. phagocytophilum and antimicrobial peptide, since it alters the taxonomic composition but not the pathway profile of the microbiome. Results also demonstrate that communities of bacteria in ticks are functionally redundant, suggesting a mechanism by ticks selecting the adequate microbiome fulfilling a core set of functions.

Original data sets

We used published 16S rRNA gene (16S) sequencing data about the taxonomic composition of the microbiomes of larvae and nymphs of I. scapularis ticks subjected to disturbance by anti-tick immunity, pathogen infection, or the treatment with an antimicrobial peptide (detailed below). These datasets were generated by the same research group using comparable methodologies and were obtained by amplicon sequencing of the V4 variable region of the bacterial 16S rRNA gene, using barcoded universal primers (515F/806R) and sequenced on an Illumina MiSeq system that generated 251-base paired-end reads. The raw data were obtained in the frame of three different experiments:

i) The study by Narasimhan et al. (4) reported that the infection of larvae of I. scapularis with B. burgdorferi increased the expression of the tick protein PIXR which facilitates B. burgdorferi infection and molting of larvae. The composition of the microbiota of tick larvae fed on mice immunized with recombinant PIXR and infected with B. burgdorferi was compared with that of ticks fed on mice immunized with Ovalbumin (OVA, an antigen not found in ticks) and infected with B. burgdorferi. The authors also showed that anti-tick immunity specific to PIXR produced dense bacterial biofilms in the tick gut. These two 16S datasets are hereinafter referred to as ‘Bb-anti-PIXR’ (n=16) and ‘Bb-anti-Ova’ (n=24). Comparisons were used to represent the microbiota composition in response to anti-tick immunity as disturbing factor.

ii) The report by Abraham et al. (3) studied the changes in gut microbiota composition and biofilms of nymphs of I. scapularis infected or not with A. phagocytophilum. These two 16S datasets are hereinafter referred to as ‘Ap-infected’ (n=20) and ‘Ap-uninfected’ (n=10). Comparisons between them were used to assess the microbiota composition in response to the disturbing factor A. phagocytophilum infection.

iii) In addition, Abraham et al. (3) showed that A. phagocytophilum induces ticks to express an anti-freeze glycoprotein (IAFGP), an antimicrobial protein with the ability to alter microbiota composition (3) inducing the formation of scattered and diffused bacterial biofilms (6). A derivative peptide of IAFGP, P1 (PARKARAATAATAATAATAATAAT) affect the tick gut microbial community compared with a scrambled P1 (sP1) peptide (AATAATATAAARRAAAAPTTAKTT) (3). These two 16S datasets are hereinafter referred to as ‘P1’ (n=16) and represents an example of microbiota composition in response to antimicrobial peptides compared with its control, dataset ‘sP1’ (n=12).

Processing of original raw sequences

We performed de novo taxonomic annotation of all the 16S datasets. To this end, the sequences were downloaded from SRA repository (18), extracted, and de-interlaced in two fastq datasets containing only the first or second mate read (11) using the data analysis platform Galaxy (http://usegalaxy.org). Demultiplexed fastq files were pre-processed and analyzed using QIIME2 pipeline (v. 2019.1) (12). Briefly, the DADA2 software package (13) was used (via q2‐dada2) for denoising the fastq files and merging mate reads. All amplicon sequence variants (ASVs) were aligned with MAFFT (14) (via q2‐alignment) and used to construct a phylogeny with FastTree 2 (15) (via q2‐phylogeny). Alpha diversity metrics (i.e. Faith’s phylogenetic diversity index (16) and Pielou's evenness index (17)) and beta diversity metrics (i.e. weighted- and unweighted UniFrac (18)) were estimated using q2‐diversity plugin after samples were rarefied (subsampled without replacement). Taxonomy was assigned to ASVs using the q2‐feature‐classifier (19) classify‐sklearn naïve Bayes taxonomy classifier against the 16S SILVA database (release 132) (20), trained 99%, hence only the target sequences fragment was used (21). QIIME2 taxa barplot functions were used for viewing the taxonomic profiling by samples. The analyses of alpha and beta diversity were performed by Kruskal‐Wallis and PERMANOVA statistical tests, respectively. For taxa comparisons, abundances based on all obtained reads were used, compared between datasets in each experiment by linear discriminant analysis effect size (LEfSe) method (22).

Prediction of functional traits in tick microbiome

The16S amplicon sequences from each experiment were used to predict the metabolic profiling of each sample. PICRUSt2 (7), a robust bioinformatic tool, was used to predict the genomes of 16S amplicon sequences. The AVSs were placed into a reference tree (NSTI cut-off value of 2) contained 20,000 full 16S sequences from prokaryotic genomes, which is then used to predict individual gene family copy numbers for each AVS. The predictions are based on several gene family catalogs, Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (KO) (23), Enzyme Classification numbers (EC) and Clusters of Orthologous Genes (COGs) (24). Pathway abundances was inferred based on structured pathway mapping based on MetaCyc database (25).

The resulting functional profiles were explored by alpha and beta diversity metrics (via q2‐diversity) and compared between datasets by Kruskal‐Wallis and PERMANOVA tests, respectively. To identify pathways with differential abundance between groups (control vs disturbed), functional profiles were compared by Welch's t-test with multiple test correction by Benjamini – Hochberg FDR (false discovery) method, the analysis were performed whit the Statistical Analysis of Metagenomic Profiles (STAMP) software package (26), and visualized as volcano plot. Samples from all groups were compared by Gneiss differential abundance test (27).

Taxonomic and functional co-occurrence networks and network resistance analysis

The taxonomic profiling resulting from each dataset was used to build co-occurrence networks of interacting organisms (28). We used the SparCC method (28) implemented in the R programming environment to produce the correlations of the abundance among bacterial genera. In taxonomic networks, the strength of correlation between pairs of bacterial genera is the weight of the link between any two co-occurring bacterial genera. Similar methods were used to build functional networks of the pathways profile of the microbiomes for each dataset, as assigned by PICRUSt2. The functional networks were built using pathways co-ocurrence, with weight of the links derived from the correlations between any two pathways according to feature abundance profile. Taxonomic and functional networks were built for each experimental condition ‘Bb-anti-Ova’, ‘Bb-anti-PIXR’, ‘Ap-uninfected’, ‘Ap-infected’, ‘sP1’ and ‘P1’.

Several indexes give an indication of topology and strength of the networks. These included the number of nodes and edges, weighted degree (WD), diameter of the network (a measure of its cohesiveness), modularity (an index that computes how interacting nodes arrange in modules), and the clustering coefficient (a measure of the density of each network). Changes in these indexes can explain the robustness of the network, the density of its links, or the nodes that tend to co-occur more frequently. Modularity maximization was used to extract modules from the networks (29). Most of the calculations were done with the software Gephi 0.9.2 (30).

Attack tolerance was tested only on the taxonomic networks. The purpose was to measure their resistance to systematic removal of nodes, either by a random attack with 100 iterations, or by a directed attack, removing the nodes according to its value of centrality (the highest, the first). The number of secondary removals (i.e. nodes that became disconnected from the network after each node removal) was also calculated. The analyses of the networks resistance were done with the package NetSwan (31) for R.

To test functional redundancy of tick microbial communities, we selected the top three bacterial modules within each taxonomic network. These modules include not less than 75% of the total bacterial genera detected in the network. The ASVs corresponding to each module were filtered in the raw dataset using q2-taxa plugin (13). The functional predictions for each module were produced and the functional profiles compared among modules using non-metric multidimensional scaling (NMDS) based on Bray-Curtis distance.

Identification of the core microbiome

We conducted a step-forward analysis to detect and characterize the taxonomic or functional cores of I. scapularis microbiome. They were defined as the set of bacterial genera or metabolic pathways that were shared by most of the individual ticks in all datasets. First, we used a Venn diagram to identify the features present in at least one sample in each dataset. The second method was based on identifying the core features which persisted across serial fractions of the samples, performed with the Qiime2 plugin feature-table (core-features) in (12). Finally, we integrated presence with feature abundance and connectedness to better characterize the core microbiome. We used a set of fuzzy logic rules (32) to calculate the probability of each bacterial genera or pathway to belong to the taxonomic or functional core. Fuzzy logic rules are not binary, but a range of probabilities. The maximum probability 1 indicates that a genus or a pathway ‘completely belongs’ to the core, while decreasing probabilities indicate a smaller probability to belong to the core. Fuzzy logic rules were built using the software Manifold v.8 (http://www.manifold.net). The degree of membership of a member to the core (either taxonomic or functional) was defined by the rules: (i) appears in the highest percent of ticks of the dataset, (ii) the number of reads is over the median of the distribution of reads detected in the dataset, (iii) the WD and centrality of the member is over the median of all the co-occurring members of the dataset, and (iv) the conditions before are not affected under disturbance. A member will have a maximum probability to be part of the core when all the conditions above are met simultaneously. As one or several of the rules above are not met, the probability of membership decreases.

16S: 16S rRNA gene; ASVs: Amplicon Sequence Variants; OTUs: Operational Taxonomic Units; Bb-anti-OVA: group of ticks infected with Borrelia burgdorferi feeding on mice immunized against Ovalbumin; Bb-anti-PIXR: group of ticks infected with Borrelia burgdorferi feeding on mice immunized against the protein PIXR; Ap-infected: group of ticks infected with Anaplasma phagocytophilum; Ap-uninfected: group of ticks not infected with Anaplasma phagocytophilum; IAFGP: Ixodes scapularis antifreeze glycoprotein; P1: peptide derivative of IAFGP; sP1: scrambled P1 peptide; KEGG: Kyoto Encyclopedia of Genes and genomes; EC number: Enzyme commission number; COG: Clusters of orthologous groups; MetaCyc: Database of metabolic pathways and enzymes; PCoA: Principal Coordinates Analysis; WD: Weighted Degree.

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and material

The datasets analysed during the current study are publicly available at the Sequence Read Archive (SRA) (www.ncbi.nlm.nih.gov/sra), under Bioproject IDs PRJNA232504 and PRJNA353730, respectively.

Competing interests

The authors declare that the research was conducted in the absence of any personal or financial relationship that could be construed as a potential conflict of interest.

Funding

No specific funding was obtained for this study

Author’s contribution

AEP, DO and ACC conceived the study. DO and ACC analysed the core data and obtained the results. AEP worked out the networks. All the authors discussed and integrated the flow of results as presented. All the authors drafted, wrote and approved the final version of the manuscript.

Acknowledgements

We thank members of our laboratories for fruitful discussions on this manuscript.

de la Fuente J, Antunes S, Bonnet S, Cabezas-cruz A, Domingos AG, Estrada-Peña A, et al. Tick-pathogen interactions and vector competence: identification of molecular drivers for tick-borne diseases. Front Cell Infect Microbiol. 2017;7:1–13.
Cabezas-Cruz A, Pollet T, Estrada-Peña A, Allain E, Bonnet SI, Moutailler S. Handling the microbial complexity associated to ticks. In: Abubakar DM, editor. Ticks and Tick-Borne Pathogens. 2018th ed. TechOpen; 2018. p. 1–36.
Abraham NM, Liu L, Lyon B, Yadav AK, Narasimhan S. Pathogen-mediated manipulation of arthropod microbiota to promote infection. PNAS. 2017;781–90.
Narasimhan S, Schuijt TJ, Abraham NM, Rajeevan N, Coumou J, Graham M, et al. Modulation of the tick gut milieu by a secreted tick protein favors Borrelia burgdorferi colonization. Nat Commun. 2017;8(184):1–16.
Narasimhan S, Rajeevan N, Liu L, Zhao YO, Heisig J, Pan J, et al. Gut microbiota of the tick vector Ixodes scapularis modulate colonization of the Lyme disease spirochete. Cell Host Microbe. 2014;15(1):58–71.
Heisig M, Abraham NM, Liu L, Neelakanta G, Mattessich S, Sultana H, et al. Antivirulence properties of an antifreeze protein. Cell Rep. 2014;9(2):417–24.
Douglas GM, Maffei VJ, Zaneveld J, Yurgel SN, Brown JR, Langille CMTCHMGI. PICRUSt2: An improved and extensible approach for metagenome inference Gavin. bioRxiv. 2019;1–42.
Shade A, Peter H, Allison SD, Baho DL, Berga M, Bürgmann H, et al. Fundamentals of microbial community resistance and resilience. Front Microbiol. 2012;3(DEC):1–19.
Allison SD, Martiny JBH. Resistance, resilience, and redundancy in microbial communities. Light Evol. 2009;2:149–66.
Leinonen R, Sugawara H, Shumway M. The sequence read archive. Nucleic Acids Res. 2011;39(SUPPL. 1):2010–2.
Blankenberg D, Gordon A, Von Kuster G, Coraor N, Taylor J, Nekrutenko A, et al. Manipulation of FASTQ data with galaxy. Bioinformatics. 2010;26(14):1783–5.
Bolyen E, Dillon M, Bokulich N, Abnet C, Al-Ghalith G, Alexander H, et al. QIIME 2: Reproducible, interactive, scalable, and extensible microbiome data science. Nat Biotechnol. 2019;37(8):852-7.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13(7):581–3.
Katoh K. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30(14):3059–66.
Price MN, Dehal PS, Arkin AP. FastTree 2 - Approximately maximum-likelihood trees for large alignments. PLoS One. 2010;5(3).
Faith DP. Conservation evaluation and phylogenetic diversity. Biol Conserv. 1992;61(1):1–10.
Pielou EC. The measurement of diversity in different types of biological collections. J Theor Biol. 1966;13(C):131–44.
Lozupone CA, Hamady M, Kelley ST, Knight R. Quantitative and qualitative β diversity measures lead to different insights into factors that structure microbial communities. Appl Environ Microbiol. 2007;73(5):1576–85.
Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome. 2018;6(1):1–17.
Yarza P, Yilmaz P, Pruesse E, Glöckner FO, Ludwig W, Schleifer KH, et al. Uniting the classification of cultured and uncultured bacteria and archaea using 16S rRNA gene sequences. Nat Rev Microbiol. 2014;12(9):635–45.
Werner JJ, Koren O, Hugenholtz P, Desantis TZ, Walters WA, Caporaso JG, et al. Impact of training sets on classification of high-throughput bacterial 16s rRNA gene surveys. ISME J. 2012;6(1):94–103.
Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. 2011;12(6):R60.
Kanehisa M, Goto S. KEGG : Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Tatusov RL, Galperin MY, Natale DA, Koonin E V. The COG database : a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res. 2000;28(1):33–6.
Caspi R, Billington R, Fulcher CA, Keseler IM, Kothari A, Krummenacker M, et al. The MetaCyc database of metabolic pathways and enzymes. Nucleic Acids Res. 2018;46(D1):D633–9.
Parks DH, Tyson GW, Hugenholtz P, Beiko RG. STAMP : statistical analysis of taxonomic and functional profiles. 2014;30(21):3123–4.
Morton JT, Sanders J, Quinn RA, Mcdonald D, Gonzalez A, Vázquez-baeza Y, et al. Balance trees reveal microbial niche differentiation. mSystems. 2017;2(1):1–11.
Friedman J, Alm EJ. Inferring correlation networks from genomic survey data. PLoS Comput Biol. 2012;8(9):1–11.
Newman MEJ. Modularity and community structure in networks. PNAS. 2006;19(23):8577–8582.
Bastian M, Jacomy M. Gephi: An open source software for exploring and manipulating networks. In: Proceedings of the Third International Conference on Weblogs and Social Media, ICWSM 2009, San Jose, California, USA, May 17-20, 2009. 2009. p. 4–6.
Lhomme S. NetSwan: Network strengths and weaknesses analysis. 2015. p. 1–8.
Yazdanbakhsh O, Dick S. A systematic review of complex fuzzy sets and logic. Fuzzy Sets Syst. 2018;338:1–22.
Esteban J, García-Coca M. Mycobacterium biofilms. Front Microbiol. 2018;8:1–8.
Zumsteg A, Urwyler SK, Glaubitz J. Characterizing bacterial communities in paper production—troublemakers revealed. MicrobiologyOpen. 2017;6(4):1–6.
Robert J. Palmer J, Shah N, Valm A, Paster B, Dewhirst F, Inui T, et al. Interbacterial adhesion networks within early oral biofilms of single human hosts. Appl Environ Microbiol. 2017;83(11):1–16.
Leathers TD, Bischoff KM. Biofilm formation by strains of Leuconostoc citreum and L. mesenteroides. Biotechnol Lett. 2011;33(3):517–23.
Espinoza JL, Harkins DM, Torralba M, Gomez A, Highlander SK, Jones MB, et al. Supragingival plaque microbiome ecology and functional potential in the context of health and disease. MBio. 2018;9(6):1–21.
Girvan M, Newman MEJ. Community structure in social and biological networks. Proc Natl Acad Sci U S A. 2002;99(12):7821–6.
Limoli DH, Jones CJ, Wozniak DJ. Bacterial extracellularpolysaccharides in biofilm formation and function. Microbiol Spectr. 2015;3(3):1–19.
Brackman G, Breyne K, De Rycke R, Vermote A, Van Nieuwerburgh F, Meyer E, et al. The quorum sensing inhibitor hamamelitannin increases antibiotic susceptibility of Staphylococcus aureus biofilms by affecting Peptidoglycan biosynthesis and eDNA Release. Sci Rep. 2016;6:1–14.
S. Hathroubi, M. A. Hancock, J. T. Bossé, P. R. Langford, Y. D. N. Tremblay, J. Labrie MJ. Surface polysaccharide mutants reveal that absence of o antigen reduces biofilm formation of actinobacillus pleuropneumoniae. Infect Immun. 2016;84(1):127–37.
May T, Okabe S. Enterobactin is required for biofilm development in reduced-genome Escherichia coli. Environ Microbiol. 2011;13(12):3149–62.
Harrison F, Buckling A. Siderophore production and biofilm formation as linked social traits. ISME J. 2009;3(5):632–4.
Layeghifard M, Hwang DM, Guttman DS. Disentangling interactions in the microbiome: a network perspective. Trends Microbiol. 2017;25(3):217–28.
Berry D, Widder S. Deciphering microbial interactions and detecting keystone species with co-occurrence networks. Front Microbiol. 2014;5:1–14.
Narasimhan S, Fikrig E. Tick microbiome : the force within. Trends Parasitol. 2015;1–9.
Chicana B, Couper LI, Kwan JY, Tahiraj E, Swei A. Comparative microbiome profiles of sympatric tick species from the Far-Western United States. 2019;1–12.
Gall CA, Scoles GA, Magori K, Mason KL, Brayton KA. Laboratory colonization stabilizes the naturally dynamic microbiome composition of field collected Dermacentor andersoni ticks. Microbiome. 2017;5(1):133.
Bonnet SI, Binetruy F, Hernández-jarguín AM, Duron O. The tick microbiome: why non-pathogenic microorganisms matter in tick biology and pathogen transmission. Front Cell Infect Microbiol. 2017;7:1–14.
50. Agler MT, Ruhe J, Kroll S, Morhenn C, Kim ST, Weigel D, et al. Microbial hub taxa link host and abiotic factors to plant microbiome variation. PLoS Biol. 2016;14(1):1–31.
Obregón D, Bard E, Abrial D, Estrada-Peña A, Cabezas-Cruz A. Sex-specific linkages between taxonomic and functional profiles of tick gut microbiomes. Front Cell Infect Microbiol. 2019;9:1–16.
Díaz-Sánchez S, Estrada-Peña A, Cabezas-Cruz A, de la Fuente J. Evolutionary insights into the tick hologenome. Trends Parasitol. 2019;35(9):725–37.
Shade A, Handelsman J. Beyond the Venn diagram: The hunt for a core microbiome. Environ Microbiol. 2012;14(1):4–12.
Gardner MR, Ashby WR. Connectance of large dynamic (Cybernetic) systems: Critical values for stability. Vol. 228, Nature. 1970. p. 784.
Zimmerli W, Low PD, Waldvogel FA. Bacterial extracellular polysaccharides in biofilm formation and function. J Clin Invest. 1984;73(3):1191–200.
Ricciardi BF, Muthukrishnan G, Masters E, Ninomiya M, Lee CC, Schwarz EM. Staphylococcus aureus evasion of host immunity in the setting of prosthetic joint infection: biofilm and beyond. Curr Rev Musculoskelet Med. 2018;11(3):389–400.
González JF, Hahn MM, Gunn JS. Chronic biofilm-based infections: Skewing of the immune response. Pathog Dis. 2018;76(3):1–7.
Scherr TD, Heim CE, Morrison JM, Kielian T. Hiding in plain sight: Interplay between staphylococcal biofilms and host immunity. Front Immunol. 2014;5:1–7.
Flemming HC, Wingender J. The biofilm matrix. Nat Rev Microbiol. 2010;8(9):623–33.
Toole GO, Kaplan HB, Kolter R. Biofilm formation as microbial development. Annu Rev Microbiol. 2000;54:49–79.

Table 1. Topological features of the taxonomic networks under different disturbance factors

Topological features	B. burgdorferi infection and anti-tick immunity		A. phagocytophilum infection		Antimicrobial peptide
Topological features	Bb-anti-OVA	Bb-anti-PIXR	Ap-uninfected	Ap-infected	sP1	P1
Nodes^a	204	329	416	506	378	364
Edges^b	2522	6974	11170	43389	7976	12787
Modules^c	8	8	6	10	8	12
Network diameter^d	3	5	4	9	4	5
Average degree^e	24.73	26.82	53.71	17.15	42.20	26.29
Weighted degree^f	4.38	10.25	23.52	6.51	18.04	10.58
Average clustering coefficient^g	0.393 (8313)	0.581 (34613)	0.528 (113571)	0.376 (12102)	0.521 (61799)	0.490 (22053)

^aMetabolic pathways with at least a significant (P < 0.01) and strong (SparCC > 0.7 or < -0.7) correlation;

^bNumber of connections/correlations obtained by SparCC analysis;

^cModules are formed by a group of nodes densely connected. A higher number of modules means for a higher number of functions that do not co-occur frequently and slightly interact with other functions;

^dThe longest distance between nodes in the network. The longest the distance the less robust is the network;

^eThe average number of connections per node in the network, that is, the node connectivity;

^fWeighted degree is the sum of the weights of all links attached to node.

^gThe average clustering coefficient indicates how nodes are embedded in the network. A higher value, together with a higher number of all possible combinations (in parenthesis) among every three nodes, is indicative of a tigther network, with nodes deeply interconnected among them. High values indicate a “small-world” effect or the presence of modules with well connected nodes but isolated form other modules.

Table 2. Topological features of the functional networks under different disturbance factors

Topological features	B. burgdorferi infection and anti-tick immunity		A. phagocytophilum infection		Antimicrobial peptide
Topological features	Bb-anti-OVA	Bb-anti-PIXR	Ap-uninfected	Ap-infected	sP1	P1
Nodes^a	392	404	427	429	425	427
Edges^b	19284	44362	32654	28067	93162	66514
Modules^c	13	3	5	9	3	3
Network diameter^d	6	4	4	4	4	4
Average degree (D)^e	98.14	124.15	152.95	130.86	125.41	103.75
Weighted degree (WD)^f	74.04	86.51	103.47	85.022	82.62	67.91
Average clustering coefficient^g	0.783 (777010)	0.742 (1087589)	0.768 (1603192)	0.757 (1302711)	0.730 (986125)	0.663 (689790)

^aMetabolic pathways with at least a significant (P < 0.01) and strong (SparCC > 0.7 or < -0.7) correlation;

^bNumber of connections/correlations obtained by SparCC analysis;

^dThe longest distance between nodes in the network. The longest the distance the less robust is the network;

^eThe average number of connections per node in the network, that is, the node connectivity;

^fWeighted degree is the sum of the weights of all links attached to node.

^gThe average clustering coefficient indicates how nodes are embedded in the network. A higher value, together with a higher number of all possible combinations (in parenthesis) among every three nodes, is indicative of a tighter network, with nodes deeply interconnected among them. High values indicate a “small-world” effect or the presence of modules with well connected nodes but isolated form other modules.

Download PDF

Version 1

posted

You are reading this latest preprint version

Resistance of tick microbiome to biological disturbance

Status:

Version 1

Abstract

Figures

Background

Results

Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Tables

Supplementary Files

Status:

Version 1