# Comparison of salivary gland and midgut microbiome in the soft ticks Ornithodoros erraticus and Ornithodoros moubata

Elianne Piloto-Sardiñas, Ana Laura Cano-Argüelles, Apolline Maitre, Alejandra Wu-Chuang, Lourdes Mateos-Hernández, Alexandra Corduneanu, Dasiel Obregón, Ana Oleaga, Ricardo Pérez-Sánchez, Alejandro Cabezas-Cruz

► To cite this version:

Elianne Piloto-Sardiñas, Ana Laura Cano-Argüelles, Apolline Maitre, Alejandra Wu-Chuang, Lourdes Mateos-Hernández, et al. Comparison of salivary gland and midgut microbiome in the soft ticks Ornithodoros erraticus and Ornithodoros moubata. Frontiers in Microbiology, 2023, 14, pp.1173609. 10.3389/fmicb.2023.1173609. [hal-04146807](https://hal.inrae.fr/hal-04146807)

HAL Id: hal-04146807

https://hal.inrae.fr/hal-04146807

Submitted on 30 Jun 2023

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

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

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

Check for updates

## OPEN ACCESS

### EDITED BY

Marcos Rogério André,
São Paulo State University, Brazil

### REVIEWED BY

Isabel Kinney Ferreira De Miranda Santos,
University of São Paulo, Brazil
Sandra Antunes,
NOVA University of Lisbon, Portugal

### *CORRESPONDENCE

Ricardo Pérez-Sánchez

ricardo.perez@irnasa.csic.es

Alejandro Cabezas-Cruz

alejandro.cabezas@vet-alfort.fr

¹These authors have contributed equally to this work

RECEIVED 24 February 2023

ACCEPTED 17 April 2023

PUBLISHED 09 May 2023

### CITATION

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

### COPYRIGHT

© 2023 Piloto-Sardiñas, Cano-Argüelles, Maitre, Wu-Chuang, Mateos-Hernández, Corduneanu, Obregón, Oleaga, Pérez-Sánchez and Cabezas-Cruz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Comparison of salivary gland and midgut microbiome in the soft ticks *Ornithodoros erraticus* and *Ornithodoros moubata*

Elianne Piloto-Sardiñas¹,²†, Ana Laura Cano-Argüelles³†, Apolline Maitre²,⁴,⁵, Alejandra Wu-Chuang², Lourdes Mateos-Hernández², Alexandra Corduneanu⁶,⁷, Dasiel Obregón⁸, Ana Oleaga³, Ricardo Pérez-Sánchez³* and Alejandro Cabezas-Cruz²*

¹Direction of Animal Health, National Center for Animal and Plant Health, Carretera de Tapaste y Autopista Nacional, San José de las Lajas, Mayabeque, Cuba, ²ANSES, INRAE, Ecole Nationale Vétérinaire d'Alfort, UMR BIPAR, Laboratoire de Santé Animale, Maisons-Alfort, France, ³Parasitology Laboratory, Institute of Natural Resources and Agrobiology (IRNASA, CSIC), Salamanca, Spain, ⁴UR 0045 Laboratoire de Recherches Sur Le Développement de L'Elevage (SELMET-LRDE), INRAE, Corte, France, ⁵EA 7310, Laboratoire de Virologie, Université de Corse, Corte, France, ⁶Department of Animal Breeding and Animal Production, University of Agricultural Sciences and Veterinary Medicine, Cluj-Napoca, Romania, ⁷Department of Parasitology and Parasitic Diseases, University of Agricultural Sciences and Veterinary Medicine, Cluj-Napoca, Romania, ⁸School of Environmental Sciences University of Guelph, Guelph, ON, Canada

**Introduction:** *Ornithodoros erraticus* and *Ornithodoros moubata* are the main vectors of African swine fever virus (ASFV) and the human relapsing fever spirochetes *Borrelia hispanica* and *Borrelia crocidurae* in the Mediterranean region and *Borrelia duttoni* in continental Africa. Manipulation of the tick microbiome has been shown to reduce vector fitness and competence in tick vectors, suggesting that the identification of key microbial players associated with tick tissues can inform interventions such as anti-microbiota vaccines to block pathogen development in the midgut and/or salivary glands.

**Methods:** In this study, we analyzed and compared the microbiome of the salivary glands and midgut of *O. erraticus* and *O. moubata*. For the taxonomic and functional characterization of the tissue-specific microbiome, we used 16S rRNA amplicon sequencing and prediction of metabolic profiles using PICRUSt2. Co-occurrence networks were built to characterize the community assembly and identify keystone taxa in each tick species.

**Results:** Our results revealed differences in the composition, diversity, and assembly of the bacterial microbiome of salivary glands and midgut within each tick species, but differences were more noticeable in *O. moubata*. Differences were also found in the microbiome of each tissue, salivary gland and midgut, between species. However, the 'Core Association Networks (CAN)' analysis revealed conserved patterns of interacting taxa in tissues within and between tick species. Different keystone taxa were identified in *O. erraticus* and *O. moubata* tissues, but Muribaculaceae and *Alistipes* were found as keystone taxa in the salivary glands of both tick species which justifies their use as anti-microbiota vaccine candidates to alter the microbiome and reduce tick fitness and/or block pathogen transmission. The high similarity of predicted metabolic pathways profiles between tissues of the two tick species suggests that taxonomic variability of the microbiome is not associated with significant changes in microbial functional profiles.

**Conclusion:** We conclude that the taxonomic structure of the microbiome in *O. erraticus* and *O. moubata* is tissue-specific, suggesting niche partitioning of
bacterial communities associated to these soft ticks. However, shared keystone
taxa and conserved patterns of interacting taxa between tissues and tick species
suggest the presence of key microbial players that could be used as anti-microbiota
vaccine candidates to affect tick physiology and/or pathogen colonization.

KEYWORDS

microbiome, salivary gland, midgut, *Ornithodoros erraticus*, *Ornithodoros moubata*,
networks

# 1. Introduction

Ticks are a growing medical and veterinary concern because tick feeding cause direct harm to their hosts, but more importantly because ticks are efficient vectors of numerous pathogens that cause severe diseases in wildlife, livestock, pets and humans, resulting in significant economic losses worldwide (Rashid et al., 2019). Ticks are haematophagous arthropods belonging to two major families: Ixodidae (hard ticks) and Argasidae (soft ticks). Ixodids are usually exophilic ticks that persist on the soil and in vegetation, actively searching for hosts when the season is suitable. After attaching to a suitable host, they feed for several days, ingest large quantities of blood, and, after becoming engorged, drop from the host and return to the soil, where they moult or, in the case of females, oviposit and die. In contrast, argasid ticks are endophilic/nidicolous parasites. In nature, they live in nests and burrows of their hosts, while in synanthropic environments they colonize animal facilities and human dwellings, hiding in the floor, walls, and ceiling. There they are protected from adverse climatic conditions and have regular access to host blood. Most argasids are fast feeders, taking less than an hour to complete engorgement. They then drop off the host and hide in their shelters to moult or reproduce. Adult specimens can feed and reproduce repeatedly and are very resistant to starvation and can survive for years without food (Sonenshine and Roe, 2014).

Among the Argasid ticks, several *Ornithodoros* species are vectors of microbial pathogens that cause severe diseases. *Ornithodoros erraticus* and *Ornithodoros moubata* are the main vectors of the African swine fever virus (ASFV) and the human relapsing fever spirochetes *Borrelia hispanica* and *Borrelia crocidurae* in the Mediterranean and *Borrelia duttoni* in mainland Africa (Arias et al., 2018; Talagrand-Reboul et al., 2018). The presence of *Ornithodoros* vectors in anthropogenic environments is a challenge for the eradication of ASFV and tick-borne relapsing fever (TBRF) in endemic areas, especially in Africa. These soft ticks also contribute to the spread of ASFV to new territories from Eastern Europe to China (Sánchez-Vizcaíno et al., 2021), as well as increase the risk of reintroduction of ASFV to countries such as Spain and Portugal where this virus has already been eradicated (Arias et al., 2018). Therefore, an effective strategy for the prevention and control of ASFV and TBRF should include the elimination of *Ornithodoros* vectors from at least the anthropogenic environment.

Most tick control strategies rely on the application of chemical acaricides, but this selects resistant tick strains and causes chemical residues to accumulate in animal products and the environment (Abbas et al., 2014). Alternative methods of tick control are urgently needed, and tick vaccines have proven to be an effective and sustainable method of controlling tick infestations and tick-borne diseases (Valle and Guerrero, 2018; Ndawula and Tabor, 2020).

Many tick vaccine development studies have focused on identifying and testing molecular components of tick physiological pathways as targets for vaccine candidates. The tick salivary glands (SG) and midgut (MG) are important components of the interface between host-tick-pathogen and play an important role in tick survival, reproduction and pathogen transmission. Accordingly, numerous tick antigens derived from these tissues have been tested as vaccine candidates (Van Oosterwijk and Wikel, 2021). Regarding *Ornithodoros* sp. vectors, previous studies has identified several salivary and intestinal antigens that could be incorporated into the development of multi-antigen anti-*Ornithodoros* vaccines, furthering the search for new and more effective antigen targets (Díaz-Martín et al., 2015; Obolo-Mvoulouga et al., 2018; Pérez-Sánchez et al., 2019a,b).

Currently, advances in defining the tick microbiome are providing new insights into tick physiology, reproduction, supply of essential molecules not synthesized by the tick itself, tick-host-pathogen interaction and vector competence (Aguilar-Díaz et al., 2021; Narasimhan et al., 2021). This increasing knowledge of microbial ecology and vector-host interactions is leading to new concepts and methods for vector and pathogen control, including anti-tick microbiota vaccines (Mateos-Hernández et al., 2021; Maitre et al., 2022; Wu-Chuang et al., 2022a). The rationale for anti-tick microbiota vaccination is that tick transmission and tick-borne pathogens can be controlled by disrupting bacterial taxa that are central to supporting other microbiome species and essential physiological functions of the tick (Maitre et al., 2022). The innovative study by Mateos-Hernández et al. (2020) showed that host antibodies induced by immunization with a tick microbiome Enterobacteriaceae caused significant mortality of engorging ticks, opening the way for the development of novel anti-tick and transmission-blocking vaccines.

Accordingly, the aim of the present study was to characterize and compare the bacterial microbiome of the SG and MG of *O. erraticus* (abbreviated here as OeSG and OeMG, respectively) and *O. moubata* (abbreviated here as OmSG and OmMG, respectively). Our results revealed that despite differences in the diversity and assembly of the tissue-associated bacterial microbiome of soft ticks there were conserved keystone taxa that could be used as anti-microbiota vaccine candidates to affect tick physiology and/or pathogen colonization.

2. Materials and methods

## 2.1. Ticks and tick material

The specimens of *O. moubata* (Om) and *O. erraticus* (Oe) used
in this study were obtained from two laboratory colonies
maintained at the IRNASA insectarium (CSIC). The *O. moubata*
colony was initiated from specimens kindly donated by Dr. Philip Wilkinson (Institute for Animal Health, Pirbright, United Kingdom), originally from Malawi (13° 59′ 00″ S 33° 47′ 00″ E). The *O. erraticus* colony was formed from specimens captured in the wild in the province of Salamanca (Spain; 40° 58′ 00″ N 5° 39′ 50″ O). Both colonies were established in the IRNASA insectarium in the 1990's. Since then, they have been routinely maintained at 28°C and 85% relative humidity and are regularly fed on rabbits.

The procedures for tick feeding and the experiments with rabbits were approved by the Ethical and Animal Welfare Committee of the Institute of Natural Resources and Agrobiology (IRNASA) and the Ethical Committee of the Spanish National Research Council (CSIC, Spain; Permit Number 742/2017) and complied with the relevant EU legislation (Directive 2010/63/EU).

## 2.2. DNA extraction from midgut and salivary glands

Prior to dissection and deoxyribonucleic acid (DNA) extraction, the ticks were washed once with ethanol 70% and once with hydrogen peroxide 3%, to disinfect their surface. Unfed female ticks of *O. moubata* (*n* = 100) and *O. erraticus* (*n* = 100) were dissected in ice-cold PBS, pH 7.4, and their salivary glands (SG) and midgut (MG) were carefully removed for genomic DNA extraction using the NucleoSpin® Tissue Kit (Macherey-Nagel, Germany). Ten samples of 10 pairs of SG or 10 MG *per* sample from each species were resuspended separately in T1 buffer and Proteinase K, then mechanically homogenized and incubated overnight at 56°C. Genomic DNA was then extracted according to the kit instructions. DNA concentration and quality were determined by spectrophotometry using the NanoDrop 2000c (Thermo Scientific, USA) and agarose gel electrophoresis.

## 2.3. 16S RRNA amplicon sequencing and processing of raw sequences

Genomic DNA was extracted and purified as previously described. A single lane of the Illumina MiSeq system was used to generate 251-base paired-end reads from variable region V4 of the 16S rRNA gene using barcoded universal primers (515F/806R) in ticks. The paired 16S rRNA raw sequences obtained from the OeMG, OeSG, OmMG, and OmSG samples were deposited in the SRA repository (Bioproject No. PRJNA931807).

Analysis of 16S rRNA sequences was performed using the QIIME 2 pipeline (v. 2021.4; Bolyen et al., 2019). The raw sequences (demultiplexed in fatsq fles), were denoised, quality trimmed and merged using the DADA2 software (Callahan et al., 2016) implemented in QIIME2 (Bolyen et al., 2019). The obtain amplicon sequence variants (ASVs) were aligned with q2-alignment of MAFFT (Katoh et al., 2002) and used to generate a phylogeny with q2-phylogeny of FastTree 2 (Price et al., 2010). Taxonomy was assigned to ASVs using a classify-sklearn naïve Bayes taxonomic classifier based on SILVA database (release 138; Bokulich et al., 2018). Only the target sequence fragments were used for the classifier (i.e., the classifier was trained with primers 515F/806R; Werner et al., 2012; Ren and Wu, 2016).

## 2.4. Diversity and differential taxonomic composition

The composition and diversity of the tissue-associated microbiome were compared within and between tick species. Microbial diversity analyses were performed using rarefied ASV tables calculated with the q2-diversity plugin in QIIME 2 (Bolyen et al., 2019). Alpha diversity richness was explored using observed features (DeSantis et al., 2006) and Faith's phylogenetic diversity index (Faith, 1992), while evenness was explored with the Pielou evenness index (Pielou, 1966). Differences in alpha diversity metrics between groups were assessed with the Kruskal-Wallis test ($p \le 0.05$) using QIIME 2 (Bolyen et al., 2019). Bacterial beta diversity was assessed with the Bray-Curtis dissimilarity index (Bray and Curtis, 1957) and compared between groups with the PERMANOVA test ($p \le 0.05$) on QIIME 2. Beta dispersion was calculated using the betadisper function and the Vegan script implemented in RStudio (Oksanen et al., 2021). Statistical analyses of beta dispersion were performed using an ANOVA test ($p \le 0.05$). Cluster analysis was performed with the Jaccard coefficient of similarity using Vegan (Oksanen et al., 2021) implemented in RStudio (RStudio Team, 2020). Jaccard distance is represented between 0 and 2, and lines are proportional to this distance.

Differences in taxa abundance between the four conditions were tested using a Kruskal Wallis test and implemented using the ANOVA-like differential expression (ALDEx2) package (Fernandes et al., 2013) on RStudio (RStudio Team, 2020). The ALDEx2 method involves constructing Monte Carlo samples of Dirichlet distributions to account for the uncertainty in the number of reads in each sample. Proportional data were then transformed using the centered log ratio (clr) transformation, so that standard statistical methods could be used. The identified differentially abundant taxa were used to create a heatmap in the RStudio (RStudio Team, 2020).

## 2.5. Inference of bacterial co-occurrence networks

Co-occurrence networks were created for each dataset using the taxonomic profiles at genera level. The networks provide a graphical representation of the microbial communities, with nodes representing taxa and edges indicating significant correlations between them. Analyses of significant positive (weight > 0.75) or negative (weight < -0.75) correlations were performed using the Sparse Correlations for Compositional data (SparCC) method (Friedman and Alm, 2012), implemented in RStudio (RStudio Team, 2020). Visualization and measurement of topological features (i.e., number of nodes and edges, network diameter, modularity, average degree, weighted degree and clustering coefficient) of the networks were performed using Gephi v0.9.2 (Bastian et al., 2009). Four Core Association Networks (CAN) were computed to search for conserved interactions in networks of different species and same tissues (i.e., OeMG vs. OmMG, and OeSG vs. OmSG) and networks of the same species and different tissues (i.e., OeMG vs. OeSG, and OmMG vs. OmSG) using the software toolbox, anuran (a toolbox with null models for identification of nonrandom patterns in association networks; Röttjers et al., 2021), and in a version available for Python 3.6.
### 2.5.1. Differential network analysis

To compare OmMG, OmSG, OeSG, and OeMG networks, a statistical network estimation analysis was performed using the network construction and comparison for microbiome (NetCoMi) method (Peschel et al., 2021) implemented in RStudio (RStudio Team, 2020). To test for dissimilarities between the two networks (i.e., OmMG vs. OmSG), the Jaccard index was calculated to test for dissimilarities between nodes in the two networks for degree, betweenness centrality, closeness centrality and eigenvector centrality. The Jaccard index tests for the similarity between sets of “most central nodes” of networks, which are defined as those nodes with a centrality value above the empirical 75% quartile. This index expresses the similarity of the sets of most central nodes as well as the sets of hub taxa between the two networks. The Jaccard index ranges from 0 (completely different sets) to 1 (sets equal). The two *value of ps* P ($J \le j$) and P ($J \ge j$) for each Jaccard index are the probability that the observed value of Jaccard's index is ‘less than or equal’ or ‘higher than or equal’, respectively, to the Jaccard value expected at random which is calculated taking into account the present total number of taxa in both sets (Real and Vargas, 1996). The ARI was calculated to test the dissimilarity of clustering in the networks. The ARI values range from -1 to 1. Negative and positive ARI values mean lower and higher than random clustering, respectively. An ARI value of 1 corresponds to identical clustering, and 0 to dissimilar clustering. The value of p tests whether the calculated value is significantly different from zero (Peschel et al., 2021).

### 2.5.2. Keystone taxa identification

Keystone taxa were defined based on three criteria, as previously reported (Mateos-Hernández et al., 2021): (i) ubiquitousness (microbial taxa present in all samples in an experimental group), (ii) eigenvector centrality higher than 0.75, and (iii) high mean relative abundance (i.e., higher than that of the mean relative abundance of all taxa in an experimental group). The eigenvector centrality measures the influence of a node in a network, a high eigenvector score means that a node is connected to many nodes which themselves have high scores (Ruhnau, 2000). Eigenvector centrality values were exported from Gephi (v 0.9.2) software (Bastian et al., 2009). For each sample in both datasets, the mean crl value was calculated and plotted together with the eigenvector centrality values using GraphPad Prism (v 8.0.0; GraphPad Software, San Diego, California USA).

### 2.5.3. Network robustness in nodes removal

In this analysis, the proportion of removed nodes required to achieve a connectivity loss of 0.80 was recorded for the OeMG, OeSG, OmMG, and OmSG networks after directed, cascading, degree or random node removal. Directed removal of nodes consists in removing first the nodes with higher betweenness centrality. The cascading effect is that the nodes with high betweenness centrality are removed first, but recalculated each time a node is removed. The robustness of the network was also measured after removing nodes with very high degree. The last type is a random removal of nodes. The robustness of the networks is calculated using the Network Strengths and Weaknesses Analysis (NetSwan) script (Lhomme, 2015) in RStudio (RStudio Team, 2020).

## 2.6. Prediction and comparison of functional traits in the tick microbiome

The 16S rRNA amplicon sequences from each dataset were used to predict the metabolic profile of each sample. PICRUST2 (Douglas et al., 2020) implemented in QIIME2, was used to predict the metagenomes from the 16S rRNA amplicon sequences. The ASVs were inserted into a reference tree (NSTI cut-of value of 3) containing more than 20,000 complete 16S rRNA sequences from prokaryotic genomes and then used to predict the copy numbers of each gene family for each ASV. The predictions are based on Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (KO; Kanehisa, 2000). Pathways were generated based on the Meta-Cyc database (Caspi et al., 2018).

The differentially frequent paths were identified based on log₂ fold change (LFC), using the Wald test implemented in the DESeq2 (v1.38.2) R package (Love et al., 2014). The DESeq2 pipeline was developed for differential gene analysis on smaller RNA-seq datasets with high overdispersion. The method assumes that the count data has a negative binomial distribution and models it using the maximum likelihood estimation (MLE) method. It includes a dispersion shrinkage estimation and fold-change estimation for each feature, resulting in accurate false discovery rates and FDR-corrected value of ps (Love et al., 2014).

# 3. Results

## 3.1. Differential bacterial composition and abundance in salivary glands and midguts of *Ornithodoros erraticus* and *Ornithodoros moubata*

The alpha diversity of the salivary gland (SG) and midgut (MG) microbiome of *O. erraticus* (Oe) and *O. moubata* (Om) was compared using three diversity metrics (i.e., Pielou evenness index, Faith's phylogenetic diversity index and observed features). Significant differences were found in the bacterial evenness of OeMG and OeSG (*p* < 0.05, Figure 1A), and in the phylogenetic diversity of OeMG and OmMG (*p* < 0.05, Figure 1B; Supplementary Table S1). No significant difference was found in the observed features within and between tissues of the two tick species (Figure 1C). Bray-Curtis index comparison showed that the composition of the tissue-associated microbiome was significantly different between tissues (SG and MG) in *O. moubata* and between species for the same tissue (SG or MG; PERMANOVA, *p* < 0.001, Figure 1D; Supplementary Table S1). However, Bray-Curtis index did not differ between tissues in *O. erraticus*. No significant differences were found in the beta dispersion of pairwise and global comparisons (ANOVA test, *p* > 0.05, Figure 1D), suggesting homogeneity of data dispersion among groups.

Analysis of the bacterial composition in tissues of *O. erraticus* and *O. moubata* revealed a total of 1,146 bacterial taxa, of which 507 genera and families (44.2%) were shared by all samples (Figure 2A). In *O. moubata*, analysis revealed 82 unique bacterial taxa (7.16%), 56 in the MG (4.89%), and 26 in the SG (2.27%). Only two genera (0.17%) were identified as unique to the OeSG (Figure 2A; Supplementary Table S2). Hierarchical clustering of samples based on Jaccard distance showed that the microbiome of OmMG and OmSG
<figure><img src="image_3.png" /><figcaption>FIGURE 1</figcaption></figure>

Comparison of microbial diversity for tissue-associated microbiome within and between *Ornithodoros erraticus* and *Ornithodoros moubata*. Comparison of alpha diversity for OeSG, OeMG, OmMG and OmSG samples, (A) Pielou evenness index, (B) Faith's phylogenetic diversity (PD), (C) observed features. (D) Comparison of beta diversity with Bray-Curtis dissimilarity index for OeSG, OeMG, OmMG, and OmSG. Beta dispersion of four sets of samples (global comparison). Small circles, crosses and triangles represent samples, and ellipses represent dispersion of samples. This test use Principal Coordinate Analysis (PCoA), it is used to explore and to visualize variability in a microbial community. ANOVA test was performed and showed that beta dispersion of the two sets of samples (four conditions) is not significantly different (*p*=0.849).

clustered separately, while OeSG and OeMG clustered closely together (Figure 2B). Analysis of differential relative abundance (expressed as clr) showed that 30 taxa changed significantly between the four datasets (Figure 2C). Overall, these results suggest that soft tick microbiome is species-specific, and that differences in the bacterial composition and diversity between tissues is higher in *O. moubata* than in *O. erraticus*.

## 3.2. Assembly of microbial communities in the microbiome of salivary glands and midguts in *Ornithodoros erraticus* and *Ornithodoros moubata

Co-occurrence networks were constructed to characterize the assembly of the bacterial community in the four groups (Figure 3; Table 1). The networks showed that each tissue of the two tick species had different patterns of bacterial co-occurrence. The highest number of total and connected nodes was observed for OmMG (Figure 3A), while the highest number of edges, and the highest proportion of positive interactions were observed in OmSG (Figure 3B) and OeSG (Figure 3C) networks, respectively. The OeMG network (Figure 3D) has the lowest number of total nodes and connected nodes, positive interactions and edges. The average values for degree and clustering coefficient were higher for the OmMG network. The modularity value of the OmSG network was twice that of the OeSG network, while the diameter values were similar in all networks (Table 1). The highest

number of unique nodes was found in OmMG, followed by OmSG, OeSG and OeMG, while 54 nodes were shared by all networks (Figure 3E; Supplementary Table S3).

The CAN analysis revealed 60 core associated taxa between OeMG and OeSG networks (Supplementary Figure S1A), and 14 core associated taxa between OmMG and OmSG networks (Supplementary Fig. S1B). The core associations between the SG of the two species consisted of 36 nodes (Supplementary Figure S1C), while core associations of 16 nodes were found between the MG networks of the two species (Supplementary Figure S1D). Seven and 11 taxa were classified as keystone in the MG and SG microbiome of *O. erraticus*, respectively. On the other hand, 16 and seven taxa met the keystone criteria in the microbiome of OmMG and OmSG, respectively (Figures 4A–D; Supplementary Table S4). Two taxa (Muribaculaceae and *Alistipes*) were found as keystone in the microbiome of OeMG, OeSG and OmSG. Fourteen keystone taxa (i.e., *Romboustia*, *Alistipes*, *Ralstonia*, *Lachnospiraceae*, *Lachnospiraceae_NK4A136_group*, *Akkermansia*, *Clostridia_UCG-014*, *Muribaculaceae*, *Lactobacillus*, *Cupriavidus*, *Eubacterium coprostanoligenes* group, *Parabacteroides*, *Dubosiella* and *Bifidobacterium*) were found as nodes in the core association networks.

Local centrality measures were compared between networks using the Jaccard index comparison test. The Jaccard index values of the comparisons OmMG vs. OmSG and OeMG vs. OeSG were lower [P (≤Jacc) < 0.05] and higher [P (≥Jacc) < 0.05] than expected by random, respectively, for all of centrality measures except for betweenness centrality in the OeMG-OeSG comparison
<figure><img src="image_4.png" /><figcaption>**FIGURE 2**
Comparison of microbial composition and abundance for the tissue-associated microbiome within and between *Ornithodoros erraticus* and *Ornithodoros moubata*. (A) Venn diagram displaying the comparison of taxa composition for all taxa. (B) Jaccard clusterization of OeMG, OeSG, OmMG and OmSG samples. The Jaccard distance is represented between 0 and 2, the lines are proportional to this distance. Each sample is represented by a color (legend). (C) Comparison of relative abundance for the tissue-associated microbiome within and between *O. erraticus* and *O. moubata*. Dendrogram heatmap resulting from the heatmap.2 function implemented on R studio. The taxa were clustered based on relative abundance (calculated as clr transformed values). Each column represents the clr values for bacterial taxa per sample and per group. Each line represents bacterial taxa with significant changes between the two datasets. Color represent the clr value.</figcaption></figure>

which has a random distribution (Table 2). The centrality measures in the MG of the two tick species were lower than expected by random [P (≤Jacc) < 0.05, Table 3], while in the SG they were higher than expected by random [P (≥Jacc) < 0.05, Table 3], except for betweenness centrality, which followed a random distribution (Table 3). The comparison of node clustering between the networks showed the higher Adjusted Rand Index (ARI) value for OeMG compared to OeSG (ARI = 0.60), while the other comparisons showed low similarities in clustering (ARI ≤ 0.40; Table 4).

The robustness of the networks in the tissues of *O. erraticus* and *O. moubata* was analyzed. In all networks, there is a rapid loss of connectivity after a cascading attack, whereas after a random attack, a larger proportion of nodes had to be removed to achieve a complete loss of connectivity (Supplementary Figure S2). After directed attacks by betweenness (Figure 5A), cascading (Figure 5B), and degree (Figure 5C) the OmSG network was less robust than the other networks, as a smaller proportion of nodes removal resulted in a

greater loss of connectivity. The OmMG network was more robust after directed attacks by betweenness (Figure 5A) and cascading (Figure 5B), while the OeSG network was more robust to degree attack (Figure 5C). All networks showed a loss of connectivity after a random attack in which a similar proportion of nodes were removed (Figure 5D).

## 3.3. Functional profiles of microbiome in salivary glands and midguts of *Ornithodoros erraticus* and *Ornithodoros moubata*

Changes in composition and structure of the microbial community may affect the inferred functional profile of the tick microbiome. We analyzed and compared the composition, diversity and differential abundance of predicted metabolic pathways in SG and MG of *O. erraticus* and *O. moubata*. Significant differences were
<figure><img src="image_5.png" /><figcaption>TABLE 1 Topological features of the networks.</figcaption></figure>

<table><thead><tr><th>Topological features</th><th>OeMG</th><th>OeSG</th><th>OmMG</th><th>OmSG</th></tr></thead><tbody><tr><td>Total nodes</td><td>821</td><td>874</td><td>832</td><td>1,001</td></tr><tr><td>Connected nodes</td><td>181</td><td>196</td><td>266</td><td>294</td></tr><tr><td>Edges</td><td>584</td><td>940</td><td>1783</td><td>793</td></tr><tr><td>Positives</td><td>449 (76.9%)</td><td>889 (94.6%)</td><td>1,312 (73.6%)</td><td>505 (63.7%)</td></tr><tr><td>Negatives</td><td>135 (23.1%)</td><td>51 (5.4%)</td><td>471 (26.4%)</td><td>288 (36.3%)</td></tr><tr><td>Modularity</td><td>0.93</td><td>0.41</td><td>0.87</td><td>2.02</td></tr><tr><td>Network diameter</td><td>12</td><td>11</td><td>11</td><td>15</td></tr><tr><td>Average degree</td><td>6.45</td><td>9.6</td><td>13.41</td><td>5.4</td></tr><tr><td>Weighted degree</td><td>2.86</td><td>7.1</td><td>5.18</td><td>1.2</td></tr><tr><td>Clustering coefficient</td><td>0.49</td><td>0.53</td><td>0.57</td><td>0.49</td></tr></tbody></table>

found in the observed features of metabolic pathways in the microbiome of OeMG and OeSG (Figure 6A; Supplementary Table S5A). The Bray-Curtis index and a PERMANOVA test showed that the beta diversity of metabolic pathways was significantly different between all datasets (PERMANOVA, *p* < 0.05), except in the SG of the two tick species (Figure 6B; Supplementary Table S5B).

The four datasets shared 421 (93.1%, 452 total) metabolic pathways. Six unique metabolic pathways (1.32%, 452 total) were

found in each of the OmSG and OmMG microbiomes (Figure 6C; Supplementary Table S6), while no unique pathway was found for OeSG and OmMG. Three others pathways (0.66%, 452 total) were shared between MG and SG; crotonate fermentation, tylosin biosynthesis and dTDP-6-deoxy-deoxy-α-D-allose biosynthesis. The biphenyl degradation pathway was found in the microbiome of *O. erraticus* (SG-MG; Figure 6C; Supplementary Table S6).

Differential analysis revealed significant changes in the abundance of 146 and 105 metabolic pathways between tissues of *O. erraticus*
<figure><img src="image_6.png" /><figcaption>FIGURE 4</figcaption></figure>

Keystone taxa in microbial communities. (A) OeMG, (B) OeSG, (C) OmMG and (D) OmSG datasets, showing the mean relative abundance (calculated as clr transformed values), and the eigenvector centrality of ubiquitous (present in all samples) bacterial taxa. The vertical dotted line represents the eigenvector centrality value cut-off set at 0.75 and the horizontal dotted line represents the average clr value of all taxa. Taxa were considered keystone when they were ubiquitous and with clr and eigenvector centrality values above these thresholds. Yellow dots represent predicted keystone taxa. (E) Venn diagram displaying the comparison of keystone taxa identified in the different networks.

<table><caption>TABLE 2 Comparison of centrality measures within tick species.</caption><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">OmMG vs. OmSG</th><th colspan="3">OeMG vs. OeSG</th></tr><tr><th>Jacc<sup>†</sup></th><th>P (≤Jacc)</th><th>P (≥Jacc)</th><th>Jacc</th><th>P (≤Jacc)</th><th>P (≥Jacc)</th></tr></thead><tbody><tr><td>Degree</td><td>0.222</td><td>0***</td><td>1.00</td><td>0.529</td><td>1.00</td><td>0***</td></tr><tr><td>Betweenness centrality</td><td>0.186</td><td>0***</td><td>1.00</td><td>0.380</td><td>0.89</td><td>0.14</td></tr><tr><td>Closeness centrality</td><td>0.222</td><td>0***</td><td>1.00</td><td>0.529</td><td>1.00</td><td>0***</td></tr><tr><td>Eigenvector centrality</td><td>0.222</td><td>0***</td><td>1.00</td><td>0.529</td><td>1.00</td><td>0***</td></tr><tr><td>Hub taxa</td><td>0.222</td><td>0***</td><td>1.00</td><td>0.529</td><td>1.00</td><td>0***</td></tr></tbody></table>

<table><caption>TABLE 3 Comparison of centrality measures between same tissues in different tick species.</caption><thead><tr><th rowspan="2">Local centrality measures</th><th colspan="3">OeMG vs. OmMG</th><th colspan="3">OeSG vs. OmSG</th></tr><tr><th>Jacc<sup>†</sup></th><th>P (≤Jacc)</th><th>P (≥Jacc)</th><th>Jacc</th><th>P (≤Jacc)</th><th>P (≥Jacc)</th></tr></thead><tbody><tr><td>Degree</td><td>0.207</td><td>0.0e+00***</td><td>1.00</td><td>0.382</td><td>0.97</td><td>0.03*</td></tr><tr><td>Betweenness centrality</td><td>0.199</td><td>1.4e-05***</td><td>0.99</td><td>0.287</td><td>0.09</td><td>0.93</td></tr><tr><td>Closeness centrality</td><td>0.207</td><td>0.0e+00***</td><td>1.00</td><td>0.382</td><td>0.97</td><td>0.03*</td></tr><tr><td>Eigenvector centrality</td><td>0.207</td><td>0.0e+00***</td><td>1.00</td><td>0.382</td><td>0.97</td><td>0.03*</td></tr><tr><td>Hub taxa</td><td>0.207</td><td>0.0e+00***</td><td>1.00</td><td>0.382</td><td>0.97</td><td>0.03*</td></tr></tbody></table>

(Figure 7A) and O. moubata (Figure 7B), respectively (Wald test, *p* < 0.05, Supplementary Table S7). In the MG (Figure 7C) and SG (Figure 7D) of the two species, 122 and 189 pathways, respectively, were found with significant changes in abundance (Wald test, *p* < 0.05, Supplementary Table S7).

## 4. Discussion

The composition and variability of the tick microbiome have been associated with several factors including tick species, tissue type, sex, stage and environmental factors, among others (Wu-Chuang et al.,

<footnote>†Jaccard index. ***p < 0.001.</footnote>

<footnote>†Jaccard index. ***p < 0.001, *p < 0.05.</footnote>
2021). In our study, the use of ticks of the same sex (i.e., female) reared in colonies and fed on the same host (i.e., rabbits) is expected to reduce confounding factors. The reduction of covariables to test our hypothesis (i.e., the microbial communities in *O. erraticus* and *O. moubata* microbiome are tissue-specific) is a strength of our work. However, a limitation of this study is the absence of negative controls that would have allowed us to exclude potential bacterial contaminants from the analysis and obtain a better resolution of the microbial communities in the soft ticks studied. Contamination throughout the different stages of the analytical process, from DNA extraction to amplification, can potentially affect analysis of the DNA sequences as well as the interpretation of the results (Lejal et al., 2020). However, it

is noteworthy mentioning that the disinfection process used in our study may have reduced potential tick surface contaminants.

Our analysis on diversity, relative abundance and composition revealed that the bacterial communities in *O. erraticus* and *O. moubata* are species and tissue-specific. These results are consistent with previous microbiome studies on hard and soft ticks. Chicana et al. (2019) found significant differences in bacterial diversity and composition between different hard tick species including *Ixodes pacificus*, *Ixodes angustus*, *Dermacentor variabilis*, *Dermacentor occidentalis*, *Dermacentor albipictus*, and *Haemaphysalis leporispalustris*. Rojas-Jaimes et al. (2022) found differences in bacterial diversity and richness between species of the genus *Amblyomma* collected from different hosts. Particularly, the bacterial microbiome of *Amblyomma scaldaturatum* and *Amblyomma ovale* collected from *Tapirus terrestris* and *Amblyomma sabanerae* collected from *Chelonoidis denticulate* presented significant differences in diversity (Rojas-Jaimes et al., 2022). Differences in microbial community diversity, composition and relative abundance were also found in the soft tick species including *Ornithodoros cf. hasei* (Carvajal-Agudelo et al., 2022), *Ornithodoros turicata* (Barraza-Guerrero et al., 2020), *Ornithodoros maritimus* (Gomard et al., 2021) and *Argas japonicus* (Yan et al., 2019). Differences in bacterial diversity, composition and relative abundance were also found in MG and SG of *Dermacentor andersoni* (Gall et al., 2017), *Dermacentor silvarum* (Duan et al., 2020) and *Ixodes ricinus*

<table><caption>TABLE 4 Network clustering comparisons.</caption><thead><tr><td>Network comparisons</td><td>Adjusted rand index (ARI)</td><td>p-Value</td></tr></thead><tbody><tr><td>OeMG vs. OeSG</td><td>0.60</td><td>0</td></tr><tr><td>OmMG vs. OmSG</td><td>0.14</td><td>0</td></tr><tr><td>OeMG vs. OmMG</td><td>0.18</td><td>0</td></tr><tr><td>OeSG vs. OmSG</td><td>0.40</td><td>0</td></tr></tbody></table>

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 6
Pathway diversity and composition in the different groups for OeSG, OeMG, OmMG and OmSG samples. (A) Observed features (B) Comparison of beta diversity with Bray–Curtis dissimilarity index. Beta dispersion of the four datasets. ANOVA test was performed and showed that beta dispersion of the four sets of samples is not significantly different (ANOVA, *p*=0.83). (C) Venn diagram displaying the comparison of pathway.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 7
Differential abundance of predicted metabolic pathways. Volcano plot showing the differential abundance of pathways in OeSG-OeMG (A), OmMG-OmSG (B), OeMG-OmMG (C) and OmSG-OeSG (D). The pathways with significant differences (Wald test, *p*<0.05) between the groups are represented by blue dots.</figcaption></figure>
(Wiesinger et al., 2023). Differences in bacterial diversity and composition between tissues of the same tick species could be associated to differential microbe-microbe interactions (e.g., competition caused by differences in nutrient availability) and/or host-microbe interactions (e.g., tissue-specific immunity or immune evasion by microbes; Wu-Chuang et al., 2021). In our study, several genera of Proteobacteria (i.e., Rickettsiella, Cellvibrio, Methylophaga, Xanthomonadaceae, Cupriavidus, and Ralstonia) presented significant changes in relative abundance between the four datasets. Kumar et al. (2022) showed a significant difference in the relative abundance of Proteobacteria in SG and MG of unfed and partially engorged *Ixodes scapularis* ticks. In addition, Proteobacteria was also reported as the most abundant taxon in soft and hard ticks such as *A. japonicus* (Yan et al., 2019), *O. cf. hasei* (Carvajal-Agudelo et al., 2022), *Dermacentor marginatus*, *I. ricinus* and *Rhiphicephalus sanguineus* (Portillo et al., 2019). The abundance of Proteobacteria in different hard and soft tick species may be associated to the fact that the bacterial species within the ecosystem use a wide range of metabolic strategies (Zhou et al., 2020), which may decrease competition with other microorganisms for nutrient availability allowing them to survive in diverse environments. Altogether, our study supports high taxonomic variability in hard and soft tick microbiome (Wu-Chuang et al., 2021).

Bacterial assembly within arthropod vectors influence pathogen ecology and transmission (Aželytė et al., 2022). The detection of bacterial taxa in the same physical space (i.e., tick tissues) could be an indication of direct interactions between microbes (Gomard et al., 2021). Here, bacterial co-occurrence networks, a methodological framework to construct ideally defined associations, were used to characterize the bacterial associations in SG and MG of *O. erraticus* and *O. moubata*. Network analysis showed that positive correlations were more frequent in all microbiomes except for OeMG where negative correlations were more frequent. The unequal number of positive and negative interactions between *O. erraticus* and *O. moubata* suggests differences in microbial-microbial interactions in these two species. Positive associations could indicate cooperation and/or functional complementarity between different bacterial taxa in a shared niche within the host organism, while negative associations would indicate co-exclusion or competition between bacterial taxa (Kumar et al., 2022). Our results suggest that microbial competition in the MG of soft ticks is stronger in *O. erraticus* compared to *O. moubata*. Differences in community assembly in the MG of *O. moubata* and *O. erraticus* concur with dissimilarities in the MG transcriptome (Oleaga et al., 2017a, 2018) and the anatomo-physiology of the MG in these two species. Particularly, the digestive tube of *O. moubata* lacks anatomical connection between the MG and the hindgut, which prevents the passage of haematin and other toxic digestion waste to the hindgut to be excreted (Sonenshine and Roe, 2014). Different concentrations of haematin and other digestion waste may explain differences in the bacterial community assembly in the MG of *O. erraticus* and *O. moubata*. To further characterize conserved patterns of microbial interactions, a CAN analysis was performed. Results showed the existence of core associations in networks of different species and tissues. Despite differences in the networks of OeMG and OmMG, we found the existence of a conserved CAN in the MG of *O. erraticus* and *O. moubata*, suggesting conserved patterns of microbe-microbe interactions in the microbiome of the genus Ornithodoros. It is important to mention that since microbial networks alone do not provide information on the biological nature of the microbial interactions; experimental studies are needed to specifically test the biological basis of the positive and negative interactions between the taxa encountered.

Networks can also be compared in terms of their resistance to taxa removal. In our study, we used robustness, an important property of networks (Liu et al., 2017), to test the ability of the networks to withstand directed and random attacks. Loss in connectivity after directed attacks was comparable between all networks; however OmSG network displayed the lowest robustness. Significant reduction of network robustness after directed attacks was found in the microbiota of mice exposed to antibiotics (Mahana et al., 2016). This analysis suggests that specific factors decreasing network robustness might be readout of ecosystem collapse (Mahana et al., 2016). The lower robustness suggests higher susceptibility of OmSG microbiome to disturbing factors. A previous report showed that network robustness decreased only marginally in *I. scapularis* ticks infected with the tick-borne pathogen *Anaplasma phagocytophilum* (Estrada-Peña et al., 2020a), compared with other disturbance factors such as anti-tick immunity (Estrada-Peña et al., 2020a), or anti-microbiota vaccines (Mateos-Hernández et al., 2021), which perturbed the network significantly. Altogether, the results suggest that OmSG microbiome might be highly susceptible to perturbation by anti-microbiota vaccines targeting keystone taxa (see below). In addition to co-occurrence patterns, microbial networks can be used to identify keystone taxa (Banerjee et al., 2018). Keystone species play an important role in the stability of microbial ecosystems. In our study, keystone taxa were identified based on their ubiquitousness, high relative abundance, and high eigenvector centrality in the networks, as previously reported for ticks (Mateos-Hernández et al., 2020; Maitre et al., 2022). Most keystone taxa identified in the microbiome of *O. erraticus* and *O. moubata* were tissue-specific, while others were found in both the MG and SG of each species (Supplementary Table 4). Identification of the same keystone taxa in both tissues suggests the possibility of migration of key microbes between these different niches. Notably, Muribaculaceae and *Alistipes* were found as keystone taxa in the microbiome of OeMG, OeSG and OmSG. Four keystone taxa were identified in *I. scapularis* ticks exposed to different temperatures (Wu-Chuang et al., 2022b), and the functional diversity of predicted metabolic pathways of the tick microbiome was contained in these keystone bacteria and their associated taxa (Wu-Chuang et al., 2022b). The involvement of the keystone Muribaculaceae and *Alistipes* and their associated taxa in the functional stability of the microbial communities in OeMG, OeSG and OmSG guarantees further research. In addition, these two taxa could be used as candidates in anti-microbiota vaccine formulations. Previous studies showed that targeting *I. ricinus* microbiome caused significant mortality in ticks fed on α-1,3-galactosyltransferase-deficient-C57BL/6 mice immunized against the keystone taxon *Escherichia-Shigella* (Mateos-Hernández et al., 2020). Furthermore, immunization against *Escherichia-Shigella* elicited the production of host antibodies the modulated the tick microbiome (Mateos-Hernández et al., 2021).

The metabolic activity of bacterial populations in specific tissues could influence tick physiology. Differences in the beta diversity and abundance of predicted metabolic pathways suggest functional specialization between SG and MG within the same species. Differential gene expression (Oleaga et al., 2017a, 2018) and protein representation (Oleaga et al., 2017b, 2021; Pérez-Sánchez et al., 2022) between MG (Oleaga et al., 2017a,b, 2018) and SG (Oleaga et al., 2021; Pérez-Sánchez et al., 2022) of *O. erraticus* (Oleaga et al., 2018; Pérez-Sánchez et al., 2022) and *O. moubata* (Oleaga et al., 2017a,b, 2021) demonstrates functional differentiation and specialization in these tissues, which may be influenced by the composition of the microbial community and its functional contribution as a result of its metabolic activity. Despite evidence of metabolic specialization in the SG and MG of *O. erraticus* and *O. moubata*,
93.1% of the predicted metabolic pathways were shared by all the datasets. Likewise, a previous study found similar biochemical pathway composition across tick ontogeny (Zolnik et al., 2016), suggesting high functional redundancy in the tick microbiome (Estrada-Peña et al., 2020a,b).

## 5. Conclusion

This study shows that the structure and composition of the tick microbiome in the soft ticks *O. erraticus* and *O. moubata* is species and tissue specific. Ticks may be exposed to different bacterial species; but, this does not guarantee colonization. Biological interactions exist between the bacteria that are part of the tick microbiome. Aspects such as the availability of nutrients or the ability to evade the immune response in the different tissues of the tick, may explain the differences in diversity and abundance of the bacterial community. Using bacterial assemblies, we were able to characterize the bacterial associations in SG and MG of *O. erraticus* and *O. moubata*, and these were found to depend on tick species and tissue. Bacterial taxa were detected in the same tissue, which could indicate a direct association from a biological point of view. The robustness of the OmMG and OeSG networks suggests that there is no trophic dependence or competition between coexisting genera, while the similarity of biochemical pathways suggests that the loss of bacterial taxa does not result in significant changes in microbial function. Pathogens and other microorganisms coexist within ticks, and bacterial communities residing inside ticks can modulate their vectorial capacity by affecting tissue colonization by pathogens (Wu-Chuang et al., 2021). Studies on the composition of internal microbial communities provide the basis for developing new strategies to interrupt pathogen transmission by modulating the tick microbiome. For example, keystone taxa identified in the microbiome of *I. ricinus* were used as vaccine candidates to reduce tick infestation (Mateos-Hernández et al., 2020) and modulate the tick microbiome (Mateos-Hernández et al., 2021). The finding of Muribaculaceae and *Alistipes* as keystone taxa in the SG of the two tick species and the MG of *O. erraticus* justify their use as anti-microbiota vaccine candidates to target all tissues and/or for the tissue-specific manipulation to block pathogen transmission. Analysis of the differences in bacterial community composition in different tick tissues could be useful to determine the major key players in pathogen colonization and to make the manipulation of the tick microbiome more specific.

## Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/, PRJNA931807.

## References

Abbas, R. Z., Zaman, M. A., Colwell, D. D., Gilleard, J., and Iqbal, Z. (2014). Acaricide resistance in cattle ticks and approaches to its management: the state of play. *Vet. Parasitol.* 203, 6–20. doi: 10.1016/j.vetpar.2014.03.006

Aguilar-Díaz, H., Quiroz-Castañeda, R. E., Cobaxin-Cárdenas, M., Salinas-Estrella, E., and Amaro-Estrada, I. (2021). Advances in the study of the tick cattle microbiota and the influence on Vectorial capacity. *Front. Vet. Sci.* 8:710352. doi: 10.3389/fvets.2021.710352

## Author contributions

EP-S: formal analysis, visualization and writing-original draft. AC-A: investigation, review and editing. AM: formal analysis and review and editing. AW-C: formal analysis and writing-review and editing. LM-H: data curation, writing-review and editing. AC: writing-review and editing. DO: software, writing-review and editing. AO: investigation and writing-review and editing. RP-S: conceptualization, resources, supervision and writing-original draft. AC-C: conceptualization, resources, visualization, supervision and writing-original draft. All authors contributed to the article and approved the submitted version.

## Funding

UMR BIPAR was supported by the French Government's Investissement d'Avenir program, Laboratoire d'Excellence "Integrative Biology of Emerging Infectious Diseases" (grant no. ANR-10-LABX- 62-IBEID). AW-C was supported by Programa Nacional de Becas de Postgrado en el Exterior "Don Carlos Antonio López" (grant no. 205/2018). AM was supported by the 'Collectivité de Corse', grant: 'Formations superieures' (grant code. SGCE-RAPPORT No. 0300). AC-A was supported by project "CLU-2019-05-IRNASA/CSIC Unit of Excellence," granted by the Junta de Castilla y León and co-financed by the European Union (ERDF "Europe drives our growth").

## Conflict of interest

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

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2023.1173609/full#supplementary-material

Arias, M., Jurado, C., Gallardo, C., Fernández-Pinero, J., and Sánchez-Vizcaíno, J. M. C. (2018). Gaps in African swine fever: analysis and priorities. *Transbound. Emerg. Dis.* 65, 235–247. doi: 10.1111/tbed.12695

Azelytė, J., Wu-Chuang, A., Žiegytė, R., Platonova, E., Mateos-Hernandez, L., Maye, J., et al. (2022). Anti-microbiota vaccine reduces avian malaria infection within mosquito vectors. *Front. Immunol.* 13:841835. doi: 10.3389/fimmu. 2022.841835
Banerjee, S., Schlaeppi, K., and van der Heijden, M. G. A. (2018). Keystone taxa as drivers of microbiome structure and functioning. *Nat. Rev. Microbiol.* 16, 567–576. doi: 10.1038/s41579-018-0024-1

Barraza-Guerrero, S. I., Meza-Herrera, C. A., García-De la Peña, C., González-Alvarez, V. H., Vaca-Paniagua, F., Díaz-Velásquez, C. E., et al. (2020). General microbiota of the soft tick *Ornithodoros turicata* parasitizing the bolson tortoise (*Gopherus flavomarginatus*) in the Mapimi biosphere reserve, Mexico. *Biology* 9:275. doi: 10.3390/biology9090275

Bastian, M., Heymann, S., and Jacomy, M., (2009). Gephi: An open source software for exploring and manipulating networks. in: *Proceedings of the Third International Conference on Weblogs and Social Media*. (Menlo Park, California: AAAI Press), 4–6.

Bokulich, N. A., Kaehler, B. D., Rideout, J. R., Dillon, M., Bolyen, E., Knight, R., et al. (2018). Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2's q2-feature-classifier plugin. *Microbiome* 6:90. doi: 10.1186/s40168-018-0470-z

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

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

Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J., and Holmes, S. P. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. *Nat. Methods* 13, 581–583. doi: 10.1038/nmeth.3869

Carvajal-Agudelo, J. D., Ramírez-Chaves, H. E., Ossa-López, P. A., and Rivera-Páez, F. A. (2022). Bacteria related to tick-borne pathogen assemblages in Ornithodoros cf. hasei (Acari: Argasidae) and blood of the wild mammal hosts in the Orinoquia region, Colombia. *Exp. Appl. Acarol.* 87, 253–271. doi: 10.1007/s10493-022-00724-9

Caspi, R., Billington, R., Fulcher, C. A., Keseler, I. M., Kothari, A., Krummenacker, M., et al. (2018). The MetaCyc database of metabolic pathways and enzymes. *Nucleic Acids Res.* 46, 633–639. doi: 10.1093/nar/gkx935

Chicana, B., Couper, L. I., Kwan, J. Y., Tahiraj, E., and Swei, A. (2019). Comparative microbiome profiles of sympatric tick species from the far-Western United States. *Insects* 10:353. doi: 10.3390/insects10100353

Clow, K. M., Weese, J. S., Rousseaude, J., and Jardine, C. M. (2018). Microbiota of field-collected *Ixodes scapularis* and *Dermacentor variabilis* from eastern and southern Ontario, Canada. *Ticks Tick Borne Dis.* 9, 235–244. doi: 10.1016/j.ttbdis.2017.09.009

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

Díaz-Martín, V., Manzano-Román, R., Obolo-Mvoulouga, P., Oleaga, A., and Pérez-Sánchez, R. (2015). Development of vaccines against Ornithodoros soft ticks: an update. *Tick Tick Borne Dis.* 6, 211–220. doi: 10.1016/j.ttbdis.2015.03.006

Douglas, G. M., Maffei, V. J., Zaneveld, J. R., Yurgel, S. N., Brown, J. R., Taylor, C. M., et al. (2020). PICRUSt2 for prediction of metagenome functions. *Nat. Biotechnol.* 38, 685–688. doi: 10.1038/s41587-020-0548-6

Duan, D. Y., Liu, G. H., and Cheng, T. Y. (2020). Microbiome analysis of the saliva and midgut from partially or fully engorged female adult *Dermacentor silvarum* ticks in China. *Exp. Appl. Acarol.* 80, 543–558. doi: 10.1007/s10493-020-00478-2

Elbir, H., Almathen, F., and Alhumam, N. A. (2019). A glimpse of the bacteriome of *Hyalomma dromedarii* ticks infesting camels reveals human *Helicobacter pylori* pathogen. *J. Inf. Dev. Countries* 13, 1001–1012. doi: 10.3855/jidc.11604

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

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

Estrada-Peña, A., Cabezas-Cruz, A., Pollet, T., Vayssier-Taussat, M., and Cosson, J. F. (2018). High throughput sequencing and network analysis disentangle the microbial communities of ticks and hosts within and between ecosystems. *Front. Cell. Infect. Microbiol.* 8:236. doi: 10.3389/fcimb.2018.00236

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

Fernandes, A. D., Macklaim, J. M., Linn, T. G., Reid, G., and Gloor, G. B. (2013). ANOVA-like differential expression (ALDEX) analysis for mixed population RNA-Seq. *PLoS One* 8:e67019. doi: 10.1371/journal.pone.0067019

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

Gall, C. A., Scoles, G. A., Magori, K., Mason, K. L., and Brayton, K. A. (2017). Laboratory colonization stabilizes the naturally dynamic microbiome composition of field collected *Dermacentor andersoni* ticks. *Microbiome* 5:133. doi: 10.1186/s40168-017-0352-9

Gomard, Y., Flores, O., Vittecoq, M., Blanchon, T., Toty, C., Duron, O., et al. (2021). Changes in bacterial diversity, composition and interactions during the development of

the seabird tick *Ornithodoros maritimus* (Argasidae). *Microb. Ecol.* 81, 770–783. doi: 10.1007/s00248-020-01611-9

Kanehisa, M. (2000). KEGG: Kyoto encyclopedia of genes and genomes. *Nucleic Acids Res.* 28, 27–30. doi: 10.1093/nar/28.1.27

Katoh, K., Misawa, K., Kuma, K., and Miyata, T. (2002). MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. *Nucleic Acids Res.* 30, 3059–3066. doi: 10.1093/nar/gkf436

Kumar, D., Downs, L. P., Adegoke, A., Machtinger, E., Oggenfuss, K., Ostfeld, R. S., et al. (2022). An exploratory study on the microbiome of northern and southern populations of *Ixodes scapularis* ticks predicts changes and unique bacterial interactions. *Pathogens* 11:130. doi: 10.3390/pathogens11020130

Lejal, E., Estrada-Peña, A., Marsot, M., Cosson, J. F., Rué, O., Mariadassou, M., et al. (2020). Taxon appearance from extraction and amplification steps demonstrates the value of multiple controls in tick microbiota analysis. *Front. Microbiol.* 11:1093. doi: 10.3389/fmicb.2020.01093

Lejal, E., Moutailler, S., Šimo, L., Vayssier-Taussat, M., and Pollet, T. (2019). Tick-borne pathogen detection in midgut and salivary glands of adult *Ixodes ricinus*. *Parasit. Vectors* 12:152. doi: 10.1101/489328

Lhomme, S. (2015). Analyse spatiale de la structure des réseaux techniques dans un contexte de risques. *Cybergeo: Eur. J. Geogr.* 71. doi: 10.4000/cybergeo.26763

Liu, J., Zhou, M., Wang, S., and Liu, P. (2017). A comparative study of network robustness measures. *Front. Comput. Sci.* 11, 568–584. doi: 10.1007/s11704-016-6108-z

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol.* 15:550. doi: 10.1186/s13059-014-0550-8

Mahana, D., Trent, C. M., Kurtz, Z. D., Bokulich, N. A., Battaglia, T., Chung, J., et al. (2016). Antibiotic perturbation of the murine gut microbiome enhances the adiposity, insulin resistance, and liver disease associated with high-fat diet. *Genome Med.* 8:48. doi: 10.1186/s13073-016-0297-9

Maitre, A., Wu-Chuang, A., Aželytė, J., Palinauskas, V., Mateos-Hernández, L., Obregón, D., et al. (2022). Vector microbiota manipulation by host antibodies: the forgotten strategy to develop transmission-blocking vaccines. *Parasit. Vectors* 15. doi: 10.1186/s13071-021-05122-5

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

Maldonado-Ruiz, L. P., Neupane, S., Park, Y., and Zurek, L. (2021). The bacterial community of the lone star tick (*Amblyomma americanum*). *Parasit. Vectors* 14:49. doi: 10.1186/s13071-020-04550-z

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

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

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

Narasimhan, S., Swei, A., Abounameah, S., Pal, U., Pedra, J. H. F., and Fikrig, E. (2021). Grappling with the tick microbiome. *Trends Parasitol.* 37, 722–733. doi: 10.1016/j. pt.2021.04.004

Ndawula, C. Jr., and Tabor, A. E. (2020). Cocktail anti-tick vaccines: the unforeseen constraints and approaches toward enhanced efficacies. *Vaccine* 8:457. doi: 10.3390/ vaccines8030457

Obolo-Mvoulouga, P., Oleaga, A., Manzano-Román, R., and Pérez-Sánchez, R. (2018). Evaluation of the protective efficacy of *Ornithodoros moubata* midgut membrane antigens selected using omics and in silico prediction algorithms. *Ticks Tick Borne Dis.* 9, 1158–1172. doi: 10.1016/j.ttbdis.2018.04.015

Oksanen, J., Blanchet, G., and Friendly, M., (2021). Vegan: community ecology package. R package version 2.6-0. Available at: https://github.com/vegandevs/vegan.

Oleaga, A., Carnero-Morán, A., Valero, M. L., and Pérez-Sánchez, R. (2021). Proteomics informed by transcriptomics for a qualitative and quantitative analysis of the sialoproteome of adult *Ornithodoros moubata* ticks. *Parasit. Vectors* 14:396. doi: 10.1186/s13071-021-04892-2

Oleaga, A., Obolo-Mvoulouga, P., Manzano-Román, R., and Pérez-Sánchez, R. (2017a). Functional annotation and analysis of the *Ornithodoros moubata* midgut genes differentially expressed after blood feeding. *Ticks Tick Borne Dis.* 8, 693–708. doi: 10.1016/j.ttbdis.2017.05.002

Oleaga, A., Obolo-Mvoulouga, P., Manzano-Román, R., and Pérez-Sánchez, R. (2017b). A proteomic insight into the midgut proteome of *Ornithodoros moubata* females reveals novel information on blood digestion in argasid ticks. *Parasit. Vectors* 10:366. doi: 10.1186/s13071-017-2300-8

Oleaga, A., Obolo-Mvoulouga, P., Manzano-Román, R., and Pérez-Sánchez, R. (2018). De novo assembly and analysis of midgut transcriptome of the argasid tick Ornithodoros
erraticus and identification of genes differentially expressed after blood feeding. *Ticks*
*Tick Borne Dis.* 9, 1537–1554. doi: 10.1016/j.ttbdis.2018.06.018

Pérez-Sánchez, R., Carnero-Morán, A., Luz Valero, M., and Oleaga, A. (2022). A proteomics informed by transcriptomics insight into the proteome of *Ornithodoros erraticus* adult tick saliva. *Parasit. Vectors* 15:1. doi: 10.1186/s13071-021-05118-1

Pérez-Sánchez, R., Manzano-Román, R., Obolo-Mvoulouga, P., and Oleaga, A. (2019a). Function-guided selection of midgut antigens from *Ornithodoros erraticus* ticks and an evaluation of their protective efficacy in rabbits. *Vet. Parasitol.* 272, 1–12. doi: 10.1016/j.vetpar.2019.06.016

Pérez-Sánchez, R., Manzano-Román, R., Obolo-Mvoulouga, P., and Oleaga, A. (2019b). In silico selection of functionally important proteins from the mialome of *Ornithodoros erraticus* ticks and assessment of their protective efficacy as vaccine targets. *Parasit. Vectors* 12:508. doi: 10.1186/s13071-019-3768-1

Perveen, N., Muzaffar, S. B., Vijayan, R., and Al-Deeb, M. A. (2020). Microbial communities associated with the camel tick, *Hyalomma dromedarii*: 16S rRNA gene-based analysis. *Sci. Rep.* 10, 1–11. doi: 10.1038/s41598-020-74116-7

Peschel, S., Müller, C. L., von Mutius, E., Boulesteix, A. L., and Depner, M. (2021). NetCoMi: network construction and comparison for microbiome data in R. *Brief. Bioinform.* 22:bbaa290. doi: 10.1093/bib/bbaa290

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

Portillo, A., Palomar, A. M., de Toro, M., Santibáñez, S., Santibáñez, P., and Otejo, J. A. (2019). Exploring the bacteriome in anthropophilic ticks: to investigate the vectors for diagnosis. *PLoS One* 14:e0213384. doi: 10.1371/journal.pone.0213384

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

Rashid, M., Rashid, M. I., Akbar, H., Ahmad, L., Hassan, M. A., and Ashraf, K. (2019). A systematic review on modelling approaches for economic losses studies caused by parasites and their associated diseases in cattle. *Parasitology* 146, 129–141. doi: 10.1017/S0031182018001282

Real, R., and Vargas, J. M. (1996). The probabilistic basis of Jaccard's index of similarity. Olmstead R, ed. *Syst. Biol.* 45, 380–385. doi: 10.1093/sysbio/45.3.380

Ren, T., and Wu, M. (2016). PhyloCore: a phylogenetic approach to identifying core taxa in microbial communities. *Gene* 593, 330–333. doi: 10.1016/j.gene.2016.08.032

Rojas-Jaimes, J., Lindo-Seminario, D., Correa-Núñez, G., and Diringer, B. (2022). Characterization of the bacterial microbiome of Amblyomma sculpturatum and *Amblyomma ovale* collected from *Tapirus terrestris* and *Amblyomma sabanerae* collected from *Chelonoidis denticulata*, Madre de Dios, Peru. *BMC Microbiol.* 22:305. doi: 10.1186/s12866-022-02717-5

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

RStudio Team (2020). RStudio: Integrated Development for R. RStudio, PBC, Boston, MA

Ruhnau, B. (2000). Eigenvector-centrality—a node-centrality? *Soc. Networks* 22, 357–365. doi: 10.1016/S0378-8733(00)00031-9

Sánchez-Vizcaíno, J. M., Laddomada, A., and Martínez Avilés, M. (2021). Editorial: African swine fever. *Front. Vet. Sci.* 7:632292. doi: 10.3389/fvets.2020.632292

Sonenshine, D. E., and Roe, R. M. (2014). “Overview: ticks, people and animals” in *Biology of Ticks*. eds. D. E. Sonenshine and R. M. Roe, vol. 1 (Oxford, UK: Oxford University Press), 3–16.

Talagrand-Reboul, E., Boyer, P. H., Bergström, S., Vial, L., and Boulanger, N. (2018). Relapsing fevers: neglected tick-borne diseases. *Front. Cell. Infect. Microbiol.* 8:98. doi: 10.3389/fcimb.2018.00098

Valle, M. R., and Guerrero, F. D. (2018). Anti-tick vaccines in the omics era. *Front. Biosci.* 10, 122–136. doi: 10.2741/e812

Van Oosterwijk, J. G., and Wikel, S. K. (2021). Resistance to ticks and the path to anti-tick and transmission blocking vaccines. *Vaccine* 9:725. doi: 10.3390/vaccines9070725

Werner, J. J., Koren, O., Hugenholtz, P., Hugenholtz, P., DeSantis, T., Walters, Z., et al. (2012). Impact of training sets on classification of high-throughput bacterial 16s rRNA gene surveys. *ISME J.* 6, 94–103. doi: 10.1038/ismej.2011.82

Wiesinger, A., Wenderlein, J., Ulrich, S., Hiereth, S., Chitimia-Dobler, L., and Straubinger, R. K. (2023). Revealing the tick microbiome: insights into Midgut and salivary gland microbiota of female *Ixodes ricinus* ticks. *Int. J. Mol. Sci.* 24:1100. doi: 10.3390/ijms24021100

Wu-Chuang, A., Hodžić, A., Mateos-Hernández, L., Estrada-Peña, A., Obregon, D., and Cabezas-Cruz, A. (2021). Current debates and advances in tick microbiome research. *Curr. Res. Parasitol. Vector Borne Dis.* 1:100036. doi: 10.1016/j.crpvbd.2021.100036

Wu-Chuang, A., Obregon, D., Estrada-Peña, A., and Cabezas-Cruz, A. (2022b). Thermostable keystone Bacteria maintain the functional diversity of the *Ixodes scapularis* microbiome under heat stress. *Microb. Ecol.* 84, 1224–1235. doi: 10.1007/s00248-021-01929-y

Wu-Chuang, A., Obregon, D., Mateos-Hernández, L., and Cabezas-Cruz, A. (2022a). Anti-tick microbiota vaccines: how can this actually work? *Biologia* 77, 1555–1562. doi: 10.1007/s11756-021-00818-6

Yan, P., Qiu, Z., Zhang, T., Li, Y., Wang, W., Li, M., et al. (2019). Microbial diversity in the tick *Argas japonicus* (Acari: Argasidae) with a focus on Rickettsia pathogens. *Med. Vet. Entomol.* 33, 327–335. doi: 10.1111/mve.12373

Zhou, Z., Tran, P. Q., Kieft, K., and Anantharaman, K. (2020). Genome diversification in globally distributed novel marine Proteobacteria is linked to environmental adaptation. *ISME J.* 14, 2060–2077. doi: 10.1038/s41396-020-0669-4

Zolnik, C. P., Prill, R. J., Falco, R. C., Daniels, T. J., and Kolokotronis, S. O. (2016). Microbiome changes through ontogeny of a tick pathogen vector. *Mol. Ecol.* 25, 4963–4977. doi: 10.1111/mec.13832