The Microbiome and Coxiella Diversity Found in Amblyomma hebraeum and Dermacentor rhinocerinus Ticks Sampled from White Rhinoceros

TL;DR

This study explores the bacterial communities in two tick species found on white rhinoceroses and finds that both ticks host diverse bacteria, including the pathogen Coxiella burnetii.

Contribution

The study identifies distinct microbiome profiles and Coxiella diversity in two tick species associated with white rhinoceroses.

Findings

Proteobacteria dominated both tick microbiomes, with Coxiella most abundant in Amblyomma hebraeum and Rickettsia in Dermacentor rhinocerinus.

Dermacentor rhinocerinus harbored significantly greater Coxiella diversity than Amblyomma hebraeum.

Coxiella burnetii prevalence was detected at 66.1% in Dermacentor rhinocerinus and 55.8% in Amblyomma hebraeum.

Abstract

The microbiome and the prevalence of the pathogenic bacterium Coxiella burnetii in ticks associated with white rhinoceros, Ceratotherium simum, is unknown. Targeted Illumina 16S rRNA amplicon sequencing was used to characterize the bacterial microbiome diversity found within 40 Amblyomma hebraeum and 40 Dermacentor rhinocerinus ticks collected from 40 white rhinoceros individuals in the Kruger National Park, South Africa. Specific emphasis was also given to further investigate the prevalence of the pathogenic C. burnetti in these tick species. At the phylum level, Proteobacteria dominated both tick microbiomes, followed by Actinobacteria and Firmicutes; Coxiella was the most abundant genus within A. hebraeum and Rickettsia within D. rhinocerinus. While alpha diversity did not differ significantly between the two tick species, beta diversity revealed significant species-specific differences in bacterial community composition. Additionally, there was no correlation between sampling region and microbiome diversity or composition for either tick species. Twenty-five Coxiella amplicon sequence variants (ASVs) were identified, forming three distinct monophyletic Coxiella clades and a fourth single ASV lineage. The Coxiella clades showed a correlation to tick species identity with D. rhinocerinus harboring significantly greater Coxiella diversity than A. hebraeum—potentially indicative of different coevolutionary pathways between the bacteria and their respective hosts. PCR of the IS1111 transposase gene for 238 ticks detected a 66.1% (56.7–74.4%) prevalence for C. burnetii in D. rhinocerinus compared to 55.8% in A. hebraeum (46.5–64.8%). These findings support a notion that each tick species is characterized by its own microbiome community composition and that both A. hebraeum and D. rhinocerinus may act as reservoirs and potential vectors of C. burnetii to white rhinoceros.

The online version contains supplementary material available at 10.1007/s00248-025-02549-6.

Figures

Relative abundance of the 10 most abundant phyla found across the 80 tick samples sequenced in the study. *Amblyomma hebraeum* samples are highlighted in purple and *Dermacentor rhinocerinus* samples in orangeFig. 2A Relative abundance of the 10 most prevalent genera found within *Amblyomma hebraeum* and *Dermacentor rhinocerinus*, respectively. B The differential abundance of genera significantly more (red) or less (blue) abundant in *Dermacentor rhinocerinus* in comparison to *Amblyomma hebraeum*

Bayesian inference tree of the phylogenetic relationships of partial 16S rRNA *Coxiella* sequences (433 aligned nucleotide sites). Sequences from the current study are in bold, and the accession numbers of reference sequences are indicated in brackets. Branch numbers indicate Bayesian posterior probabilities (bottom) and the maximum-likelihood (ML) percentage bootstrap support values (top) for each lineage (A–D) are indicated above the branches. ASVs specific to *Amblyomma hebraeum* are highlighted in purple, those specific to *Dermacentor rhinocerinus* are highlighted in orange, and those sha

Full text

Introduction

Ticks are important vectors of disease capable of transmitting a variety of microbes to humans, livestock, and wildlife [1, 2]. These microbes play an integral role in tick fitness, survival, and vectoral capacity and include a diverse array of endosymbiotic, commensal, and pathogenic bacteria [1, 3, 4]. Genera such as Anaplasma, Borrelia, Coxiella, Ehrlichia, and Rickettsia are known to have strains that are pathogenic [2, 5]. The gram-negative intracellular bacterium, Coxiella burnetii, which forms the focus of the present study, is well known to cause coxiellosis or Q fever and can affect humans, wildlife, and livestock [6–8]. Many infected individuals remain asymptomatic, but since the main sites of chronic infection are the mammary glands and uterus [9], the disease is often associated with reproductive challenges including mastitis, metritis, weak offspring, premature delivery, infertility, stillbirth, and abortion [6].

In addition to the pathogenic C. burnetti, there are also highly prevalent Coxiella endosymbionts which have a broad tick host range [10, 11]. This lack of knowledge is particularly true for ticks occurring on wildlife [10, 12]. Interestingly, phylogenetic analysis of Coxiella sequences found within ticks revealed that the pathogenic C. burnetii likely evolved from a maternally inherited Coxiella endosymbiont [13]. As obligate blood-feeders, ticks rely on these endosymbiotic microbes to supplement nutrients and vitamins lacking in their hematophagous diet, and apart from Coxiella, other endosymbionts include Rickettsia, Francisella, and Candidatus midichloria [14].

White rhinoceros (Ceratotherium simum) are an iconic flagship species currently listed by the IUCN Red List of Threatened Species (IUCN, 2022) as near threatened [15]. Rhinoceros are susceptible to infection by C. burnetii, but the exact route of infection is poorly understood [16]. Coxiellosis was recently identified and associated with abortions in a managed group of white rhinoceros in the USA [17]. In South Africa, Donnelly et al. [16] reported an overall seropositivity of 71.0% for C. burnetii in white rhinoceros sampled in Kruger National Park (KNP; historically home to the largest white rhinoceros population globally) while individuals sampled in surrounding private reserves had a lower seropositivity of 41.1%.

Since the KNP has served as a source population for regional rhinoceros re-introductions into Africa, it is important to investigate the microbial diversity and abundance and specifically the incidence of C. burnetii in tick species associated with rhinoceros [16]. In South Africa, C. burnetii has been detected in several ixodid ticks associated with domestic animals, and these include Rhipicephalus evertsi evertsi, Rhipicephalus sanguineus, Haemaphysalis elliptica, and* Amblyomma hebraeum* [18–20]. To further our knowledge on the microbiome of ticks associated with wildlife and in particular to gain more insights into the prevalence and diversity of Coxiella, we sampled two ixodid tick species (Dermacentor rhinocerinus and A. hebraeum) commonly associated with white rhinoceros [21]. These ticks were the most abundant during sampling and the life history of the two species also differ. Amblyomma hebraeum is a generalist species (adult stages can occur on multiple different mammalian host species) while adult stages of D. rhinocerinus occur only on rhinoceros [21]. Amblyomma hebraeum is a well-known vector of disease and can harbor an array of pathogenic bacteria such as C. burnetii, Rickettsia africae, and Ehrlichia ruminantium [21–24]. While there is nothing known about the microbiome diversity or vectoral capacity of D. rhinocerinus, a study on the congeneric D. reticulatus in Europe detected Anaplasma marginale, C. burnetii as well as several species of Rickettsia, Borrelia, Francisella, and Bartonella [25].

To obtain novel data describing the microbiome composition found within A. hebraeum and D. rhinocerinus, we sampled ticks from white rhinoceros in KNP and embarked on a 16S rRNA metagenetic study. To further investigate the prevalence of the pathogenic C. burnetti in ticks associated with white rhinoceros in KNP, we used species-specific amplification of the IS1111 transposase gene. Since different tick species often differ in microbiome composition, it was hypothesized that the generalist A. hebraeum and the specialist D. rhinocerinus would differ significantly in microbiome species compositions [26–28]. In addition, since the environment off the host can play a role in shaping microbiome diversity and abundance [28–32], it was also predicted that ticks sampled from different sampling regions (landscape types) within KNP would have different microbiome diversities and abundances. Finally, it was hypothesized that Coxiella bacteria (symbionts and pathogenic strains) would be highly prevalent in both tick species [3, 11, 33] and since members of the Dermacentor genus have been found to be one of the most common carriers of the bacterium, C. burnetii would be more prevalent in D. rhinocerinus [34].

Materials and Methods

Tick Sampling

Ticks were opportunistically sampled from white rhinoceros during dehorning procedures and/or health assessments in the KNP during 2019–2023. Sampled ticks were dry-frozen at –80 °C. Tick identification, sorting, and subsequent DNA extraction took place in a sterile environment at the Veterinary Wildlife Services Laboratory, Skukuza, KNP. Ticks were identified to species level based on morphology [21] using a Nikon stereo microscope (Nikon Microscopy, Tokyo, Japan). For microbiome analysis, 40 A. hebraeum adult males and 40 D. rhinocerinus adult males were respectively sampled from 40 different white rhinoceros individuals (one tick from each species were sampled from the same rhinoceros individual). Ten of the rhinoceros individuals originated from the central region and 30 individuals from the southern region (Supplementary Fig. 1). To determine C. burnetii prevalence, an additional 158 tick samples were included (80 A. hebraeum and 78 D. rhinocerinus sampled randomly from different rhinoceros individuals).

DNA Extraction and Microbiome Sequencing

To remove external surface contamination, ticks were washed five times prior to DNA extraction. First, a 10% sodium hypochlorite (NaClO) rinse was followed by a sterile phosphate buffered saline (PBS) wash and then an absolute ethanol (100% EtOH) wash. The ticks were then double-rinsed with 1 mL of sterile PBS to remove remaining ethanol. Each wash was conducted for 1 min using 1 mL of solution and a vortex. Ticks were cut into pieces in a sterile petri dish and DNA extractions were performed in separate sterile 1.5-mL Eppendorf tubes. To prevent contamination, tick cutting and DNA extractions were performed in a laminar flow hood using separate new sterile scalpels and filter tips.

Genomic DNA was extracted using the NucleoSpin Tissue kit (Macherey–Nagel, Düren, Germany). First, 180 μL of the NucleoSpin Tissue lysis buffer (Macherey–Nagel) and 25 μL of Proteinase K (Macherey–Nagel) were added to the samples and incubated overnight at 56 °C. A negative DNA extraction (blank) without any tick tissue was also included. DNA was eluted in 100 μL elution buffer (Macherey–Nagel) and stored at − 20 °C.

Prior to sequencing and PCR screening, DNA concentration was determined using the Qubit dsDNA Broad Range Assay kit (Thermofisher Scientific, Waltham, MA, USA) and the Qubit 4 fluorometer (Thermofisher Scientific) as per the manufacturer’s instructions. DNA samples to be used for microbiome analysis were equilibrated to represent similar DNA concentration in all samples and these were analyzed by BGI Genomics, Hong Kong. The 16S rRNA V3–V4 hypervariable region was amplified using the 341 F (5′-ACT CCT ACG GGA GGC AG CAG-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′) fusion primers [35] and sequenced on an Illumina HiSeq platform (Illumina, San Diego, CA, USA).

Species-Specific Conventional PCR Detection of Coxiella burnetii

Ticks were washed twice prior to DNA extraction, first using absolute ethanol (100% EtOH), followed by sterile PBS. The primers Trans1 (5′-TAT GTA TCC ACC GTA GCC AGT C-3′) and Trans2 (5′-CCC AAC AAC ACC TCC TTA TTC-3′), which target a 687-bp fragment of the IS1111 transposase element of the C. burnetii genome [36], were used to screen the 238 tick samples (80 that were included in the microbiome study and 158 additional samples) for C. burnetii. We followed the protocol of Kamau et al. [36] with some modifications. Each 12.5-μL PCR reaction consisted of 1 μL equilibrated template DNA (the same concentration), 0.5 μL of each primer (final concentration of 0.5 μM), 6.25 μL of Taq DNA Polymerase Master Mix RED (Amplicon, Denmark), and 4.25 μL of nuclease free water. Amplification was performed on a GeneAmp PCR system 2700 thermal cycler (Applied Biosystems, USA). A total of 5 μL of each PCR product (including the positive and negative controls) was visualized on a 1% agarose gel stained with ethidium bromide. Eight randomly selected positive samples were Sanger sequenced in both directions at the Central Analytical Facility of Stellenbosch University to confirm the authenticity of the pathogenic C. burnetii. Sequences were manually edited using the BioEdit Sequence Alignment Editor (version 7.2.5) and then queried against known C. burnetii reference sequences in GenBank using the Basic Local Alignment Search Tool (BLASTn) online server. The PCR products of samples which produced any visible band at the expected amplicon size of 687 bp (also matching the fragment size of a sequenced positive control) were noted as being positive for C. burnetii. To test for PCR consistency and amplification accuracy, 30 samples (tested as either positive or negative) were randomly selected for reamplification, following the same procedure.

Sequence Data Analyses

Sequenced reads were demultiplexed and reads with a Phred quality score of less than 20 over a 25-bp sliding window were truncated, together with reads whose lengths were less than 75% of their original lengths after truncation. Sequence quality was assessed using FastQC [37] and MultiQC [38]. Thereafter, the Quantitative Insights Into Microbial Ecology (QIIME2 2023.7) bioinformatics pipeline [39] was used to perform sequence analysis on the Stellenbosch University high performance computing cluster (HPC) 2 (http://www.sun.ac.za/hpc). The DADA2 plugin with default parameters was used for denoising, paired-end merging, removal of chimeric reads, sequence dereplication, and generation of amplicon sequence variants (ASVs). The Naive Bayes classifier trained on the Silva 138 database (99% OTUs full-length sequences) was used to assign genus level taxonomy to representative ASV sequences at 99% nucleotide sequence similarity. The ASVs were then aligned using the align-to-tree-mafft-fasttree pipeline and a phylogenetic tree was created using q2-phylogeny. Finally, the filter-features method from the q2-feature-table plugin was used to remove any ASVs with a total count of less than 10 over the whole dataset.

All further downstream processing was conducted in R (version 4.3.3 (2024–02–29)). The ASV table was filtered to exclude “Archaea,” “Eukaryota,” or “Chloroplast” using phyloseq [40]. Decontam v1.22.0 (Davis et al*.*, 2018) was used to identify and remove all contaminant sequences from the dataset under a prevalence method using a threshold of 0.1. To ensure an even sampling depth, the samples were rarefied to 48,000 sequences using the phyloseq::rarefy_even_depth function after inspecting rarefaction curves (Supplementary Fig. 2). All further statistical analyses were conducted on the rarefied ASV table.

Statistical Analyses

Statistical analyses were performed using various R (version 4.3.3 (2024–02–29)) packages. Diversity analyses were performed at the ASV level using phyloseq v1.46.0 [41]. Alpha diversity was assessed for each tick species using (1) observed richness (not accounting for abundance or evenness), (2) Shannon diversity index (accounts for richness and evenness), and (3) Simpson diversity index (accounts for richness and evenness). Statistical significance for these were performed using the Wilcoxon rank-sum test. The relative abundance of the 10 most prevalent bacterial phyla as well as the 10 most prevalent bacterial genera were calculated using ASV counts. Alpha diversity measures and the relative abundances were displayed graphically using ggplot2 v3.5.1 [41]. Beta diversity, the Bray–Curtis dissimilarity index, and a principal coordinate analysis (PCoA) plot was used to visualize the bacterial community differences between the tick species. Significance for the latter was determined by using a permutational multivariate analysis of variance (PERMANOVA) conducted on the Bray–Curtis dissimilarity distances using the adonis2 function in vegan v2.6–6.1 [42]. Differentially abundant ASVs were calculated at the genus level using the Microbiome differential abundance and correlation analysis with bias correction (ANCOMBC v2.4.0) package [43]. Only ASVs present in at least 10% of the samples were included with the Holm’s method selected for p-adjustment. The tick samples were respectively grouped according to sampling region (central and southern regions of the KNP; Supplementary Fig. 1) and similar diversity analyses and statistical tests were performed as outlined above.

Coxiella Phylogenic Analyses

Phylogenetic analysis of the Coxiella ASVs was performed using maximum likelihood (ML) and Bayesian inference (BI) in MEGA version 11 [44] and MrBayes 3.2.7a [45], respectively. Additional 16S rRNA reference sequences for C. burnetii and Coxiella-like organisms (CLO) were downloaded from GenBank (Supplementary Table 1). 16S rRNA sequences of Rickettsiella (two sequences) and Legionella were included as outgroups (Supplementary Table 1). Sequences were aligned using MUSCLE in MEGA and the same software was also used to determine the best-fitting evolutionary model based on the Akaike Information Criterion. For the ML, heuristic search nodal support was estimated using 100 bootstrap replicates, and posterior probabilities for the Bayesian analyses were obtained after running 5 million generations in parallel, sampling every 50 generations with a burn-in of 5000. Consensus trees were visualized using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Coxiella burnetii Summary and Descriptive Statistics

Phyloseq was used to perform diversity analyses of Coxiella ASVs. Alpha and beta diversity calculations followed the same methodology outlined above. The number of ticks which tested positive for C. burnetii was summarized for each species, and the prevalence was expressed as a proportion by dividing the total number of positive samples by the total number of samples screened for each tick species. The prop.test function in R was used to calculate 95% confidence intervals for each species, and a chi-square test was performed to assess statistical significance. Prevalence was also summarized by region (central and southern) for each tick species, using the methods described above, with a Fisher’s exact test performed to test for significant differences.

Results

Microbiome Diversity and Relative Abundance within and among Tick Species

After filtering and correcting for contaminants, a total of 3108 ASVs with a mean sequence read length of 412.33 (± 25.9) bp were characterized, and these belonged to 306 families and 555 genera. All data are available from the NCBI Sequence Read Archive under Bioproject accession number PRJNA1220073. Proteobacteria, Actinobacteria, and Firmicutes were the dominant phyla (Fig. 1). The relative abundance of the 10 most prevalent genera differed between tick species (Fig. 2A). Rickettsia (W = 4.47, Holm’s adjusted p-value = 0.002) and Anaerolinea (W = 2.75, Holm’s adjusted p-value < 0.001) were significantly more abundant in D. rhinocerinus than in A. hebraeum (Fig. 2B). Furthermore, Helococcus, not featuring prominently in A. hebraeum, formed part of the 10 most prevalent genera within D. rhinocerinus. Amblyomma hebraeum had higher relative abundances of Corynebacterium, Coxiella, Proteus, and Streptococcus (Supplementary Fig. 3) with Proteus, nearly absent from Dermacentor. Differential abundance analysis revealed that Rickettsia and Anaerolinea were significantly more abundant in D. rhinocerinus when compared to A. hebraeum (Fig. 2B).Fig. 1. Relative abundance of the 10 most abundant phyla found across the 80 tick samples sequenced in the study. Amblyomma hebraeum samples are highlighted in purple and Dermacentor rhinocerinus samples in orangeFig. 2A Relative abundance of the 10 most prevalent genera found within Amblyomma hebraeum and Dermacentor rhinocerinus, respectively. B The differential abundance of genera significantly more (red) or less (blue) abundant in Dermacentor rhinocerinus in comparison to Amblyomma hebraeum

Although alpha diversity metrics seemed to be generally higher in A. hebraeum (Fig. 3A), there were no significant differences between the two tick species for observed richness (W = 932, p = 0.21), and the Shannon (W = 895, p = 0.37) and Simpson (W = 947, p = 0.16) diversity indexes. Beta diversity analysis using Bray–Curtis dissimilarity distances indicated a grouping or clustering of data points according to tick species, suggesting differences in community composition between the two tick species (Fig. 3B) that was significantly supported by PERMANOVA (species: F = 6.063, R^2^ = 0.073, p = 0.001)*.*Fig. 3A Alpha diversity indices (observed richness, Shannon and Simpson diversity indexes) compared between tick species. For each index, violin plots show the distribution of diversity values, with overlaid boxplots indicating the summary statistics. Individual data points within each violin plot represent sample values. Wilcoxon rank-sum test p-values were > 0.05 for all indices. B Beta diversity of the microbiome in Amblyomma hebraeum (purple) and Dermacentor rhinocerinus (orange) displayed by PCoA for distances across samples calculated using Bray–Curtis dissimilarity. Plot ellipses represent 95% confidence regions for the clusters

Regional Diversity

Alpha diversity indices showed no significant differences between A. hebraeum from the central and from the southern region in the KNP (observed richness: W = 0.163.5, p = 0.68; Shannon: W = 147, p = 0.94 and Simpson: W = 157, p = 0.84). The same result was obtained for D. rhinocerinus (observed richness: W = 158.5, p = 0.80, Shannon: W = 157, p = 0.84 and Simpson: W = 163, p = 0.70). The PCoA plot revealed a complete overlap of regional clusters, with the southern cluster falling within the central cluster for both tick species (Supplementary Fig. 4), and PERMANOVA results (region: F = 0.3738, R^2^ = 0.004, p = 0.98) confirmed that regional bacterial community compositions were not significantly different from each other.

Phylogenetic and Diversity Analyses of Coxiella ASVs

Twenty-five distinct 16S rRNA Coxiella ASV’s were identified. Overall observed richness (W = 531.5, p = 0.007) and Shannon (W = 98, p < 0.001) and Simpson (W = 94, p < 0.001) diversity of these Coxiella ASVs were significantly greater in D. rhinocerinus compared to A. hebraeum (Fig. 4A). The PCoA plot showed clear clustering and separation according to tick species (Fig. 4B) and the PERMANOVA analysis confirmed that the two tick species had significantly different Coxiella community compositions (species: F = 25.131, R^2^ = 0.243, p = 0.001).Fig. 4A Observed richness, Shannon and Simpson diversity indices for Coxiella across tick species. For each index, violin plots show the distribution of diversity values, with overlaid boxplots indicating the summary statistics. Individual data points within each violin plot represent sample values. Significant differences, as indicated by p-values from Wilcoxon rank-sum tests, are displayed above the plots. B Beta diversity of Coxiella found within Amblyomma hebraeum (purple) and Dermacentor rhinocerinus (orange) displayed by PCoA for distances across samples calculated using Bray–Curtis dissimilarity. Plot ellipses represent 95% confidence regions for the clusters

The GTR model of sequence evolution was used in the phylogenetic analyses that supported the clustering of all 25 Coxiella ASVs in a strongly supported monophyletic group (Fig. 5). Both ML and BI grouped the Coxiella ASVs into four distinct lineages (three monophyletic clades and one ASV clustered on its own) with moderate to high support for the nodes defining the four lineages (Fig. 5 and Supplementary Fig. 5). Unfortunately, the phylogenetic analyses could not conclusively resolve which of these four lineages contains the pathogenic C. burnetii (due to low bootstrap support values). The individual Coxiella ASV, labeled lineage A, was most similar to C. cheraxi. The remaining three clades (lineages B–D; Fig. 5) showed some level of association with tick species, with 13 ASVs specific to A. hebraeum, eight specific to D. rhinocerinus, and four shared by both tick species (Supplementary Fig. 6). The 13 ASVs unique to clade B were nearly all derived from A. hebraeum, while clades C and D were mostly derived from D. rhinocerinus (Fig. 5). The shared ASVs had higher prevalences relative to the species-specific ASVs in both tick species (Supplementary Fig. 7). Of the shared ASVs, “Coxiella_51843” had the highest prevalence in A. hebraeum, found in 100% of the samples, while “Coxiella_7c396” was found in 100% of the samples in *D. rhinocerinus.*Fig. 5. Bayesian inference tree of the phylogenetic relationships of partial 16S rRNA Coxiella sequences (433 aligned nucleotide sites). Sequences from the current study are in bold, and the accession numbers of reference sequences are indicated in brackets. Branch numbers indicate Bayesian posterior probabilities (bottom) and the maximum-likelihood (ML) percentage bootstrap support values (top) for each lineage (A–D) are indicated above the branches. ASVs specific to Amblyomma hebraeum are highlighted in purple, those specific to Dermacentor rhinocerinus are highlighted in orange, and those shared between both tick species are highlighted in green. Legionella and Rickettsiella sequences were used as outgroups

Coxiella burnetii Prevalence

Coxiella burnetii was detected in 145/238 individual ticks, 60.9% [95% CI, 54.4–67.1%] (Supplementary Fig. 8). BLASTn of eight sequenced amplicons revealed 99–100% sequence similarity to C. burnetii. Dermacentor rhinocerinus showed a higher mean prevalence of positive samples at 66.1% (95% CI, 56.7–74.4%) when compared to A. hebraeum with 55.8% (95% CI, 46.5–64.8%), albeit not significantly different (χ^2^ = 2.221, p = 0.14) (Supplementary Fig. 9). The overlap in confidence intervals combined with the results from the Fisher’s exact test revealed that there was no significant association between C. burnetii prevalence and sampling region (Supplementary Fig. 10).

Discussion

Bacterial Community Compositions

This study provided first insights into the rich and diverse bacterial microbiome found within A. hebraeum and D. rhinocerinus ticks that are associated with white rhinoceros in the KNP, South Africa. Although it has been predicted in the literature that a higher exposure to different hosts during different life stages (such as for A. hebraeum) will increase microbiome diversities in ticks [46, 47], this was not supported by our data. The generalist A. hebraeum with a larger host range (scrub hares, wild carnivores, small antelope, and game birds) showed a similar species diversity to the host-specific D. rhinocerinus (where adults are confined to rhinoceros and nymphs and larvae attach to rodents; [46, 47]). These findings rather support the notion that tick microbiomes are better correlated to species-specific and environmental factors than to host blood source [28]. Although differences in species diversity were not statistically significant, A. hebraeum exhibited higher species richness and evenness than D. rhinocerinus (Fig. 3A), suggesting a more stable and resilient microbiome. Greater microbial diversity and balance may enhance resistance to disturbance and support functional stability.

In contrast to the lack of significant differences in alpha diversity among the two tick species sampled herein, bacterial community compositions differed significantly between A. hebraeum and D. rhinocerinus. While Rickettsia was the most common bacterium detected in A. hebraeum sampled from cattle in a previous study in South Africa [23], Coxiella was the dominant genus in A. hebraeum collected from rhinoceros hosts in a different geographic region of the country. Rickettsia was, however, the most common and significantly more abundant genus within D. rhinocerinus. The high incidence of Coxiella and Rickettsia in A. hebraeum and D. rhinocerinus, respectively, is not surprising [1, 5, 48] and highlights the importance of these bacteria as endosymbionts for tick survival [3, 33]. While the drivers of the differences in microbial community composition among the two tick species cannot be deduced from our study design, it is tempting to speculate that exposure to different host groups occurring in different environments [47], inter-bacterial competition of the microbes making up the microbiome [49], and the time that various life stages spent feeding on the various hosts [46] could all be invoked as potential mechanisms that play a role in the detected significant differences in beta diversity among D. rhinocerinus and A. hebraeum.

The present study spans two geographic regions (southern and central regions) which are characterized by differences in precipitation and soil types that also culminate in differences in vegetation patterns and wildlife compositions [50]. A previous study on microbes (Feline Immunodeficiency Virus in lion, Panthera leo) showed regional differences among eco-regions [51] and microbiome composition has been reported to differ between geographic regions in other studies [52, 53]. In the present study, there were no significant differences between ticks sampled in the southern and central regions within KNP. At small regional scales, it has been suggested that environmental and host variables play a lesser role in microbiome composition and a few ectoparasite microbiome studies have highlighted the greater influence of species identity on microbiome composition rather than environmental conditions [26, 28, 49]. We propose that the small geographic scale of sampling in our study resulted in the absence of geographical differences in the microbiomes studied herein (also see [54]).

Coxiella Endosymbionts and Pathogens

Dermacentor rhinocerinus ticks appeared to have a higher prevalence of pathogenic C. burnetii (66.1%) when compared to A. hebraeum (55.8%), and interestingly also had a greater diversity of Coxiella ASVs than A. hebraeum (as evident in branch lengths among ASVs in Fig. 5). The in-depth analyses of the 25 distinct Coxiella 16S rRNA ASVs emphasize the high level of sequence diversity among the members of the genus, and the phylogenetic clustering into four distinct lineages/clades that are also associated with specific tick species provides good support for the notion that symbiotic Coxiella most likely shared a long convergent evolutionary history with Amblyomma and Dermacentor, respectively*.* Beta diversity analysis further supported the high level of variation by also displaying significant differences in the community compositions of Coxiella ASVs among the tick species. The high prevalence of Coxiella within Amblyomma and Dermacentor may be explained by the Coxiella endosymbionts being maternally inherited [13]. Indeed, vertical transmission of Coxiella symbionts through the egg has been shown in Amblyomma cajennense [55] and Duron et al. [13] also found closely related Coxiella endosymbionts present among closely related tick species. However, since the host rhinoceroses were not tested in this study for their microbial community, the transfer of microbial communities between ticks and hosts remains unclear and requires further study.

The greater percentage of species-specific Coxiella ASVs detected in A. hebraeum compared to D. rhinocerinus may be linked to the high functional diversity associated with tick life history*. Coxiella* symbionts would be beneficial for a generalist species—enabling it to feed on multiple host blood sources while obtaining the necessary nourishment. Important to note, however, these unique ASVs had a low prevalence and indicates that they are not widespread among all the individuals sampled herein. While D. rhinocerinus had fewer species-specific Coxiella ASVs, alpha diversity measures were significantly greater than those of A. hebraeum. These results provide circumstantial evidence that Coxiella endosymbionts are not passed on from the rhinoceros to the ticks and rather reflect a pattern of host specific co-evolution enforced by vertical transmission within tick species [55].

While it remains untested whether ticks became infected with the pathogenic C. burnetii after feeding on an infected rhinoceros host, and/or whether they were carrying the bacterium before attaching to the rhinoceros, variation in host preference during the juvenile life stages of the ticks, and the associated co-evolution between symbiont and host, could provide a potential reason for the higher overall prevalence of Coxiella observed in D. rhinocerinus [56]. In fact rodents, the preferred hosts for immature D. rhinocerinus, are known to play a significant role in C. burnetii transmission [21, 57] while there is less known about the role of hosts preferred by immature A. hebraeum (e.g., scrub hares, wild carnivores, small antelope, and game birds).

The overall prevalence of C. burnetii positive tick samples at 60.9% is slightly lower than the 71% seropositivity previously detected in white rhinoceros sampled in KNP between 2017 and 2019 [16]. This prevalence is higher than that reported for dogs (41.0%; [19]) and domestic ruminants (20.0–68.0%; [19]) in South Africa but to be expected since wild animals often have higher tick burdens than domestic livestock [58]. A higher tick burden could provide more opportunity for C. burnetii transmission and persistence within the tick vector populations [16]. Although it is not possible to directly link individual rhinoceros seropositivity to the prevalence detected in the ticks in the current study, results indicate that the potential for C. burnetii infection in the KNP white rhinoceros population is high. Circumstantial evidence based on overall Coxiella prevalence throughout the park suggests that the pathogenic bacterium may have a similar environmental load throughout the park; however, more research targeting more tick species as well as wildlife host species (and including the northern region of KNP) is required to confirm this observation.

Supplementary Information

Below is the link to the electronic supplementary material.Supplementary file1 (DOCX 2300 KB)