The composition of tick microbiomes varies both within and among tick species. Whether this variation is intrinsic (related to tick characteristics), or extrinsic (related to vertebrate host and habitat) is poorly understood but important, as microbiota can influence the reproductive success and vector competence of ticks. We aimed to uncover what intrinsic and extrinsic factors best explain the microbial composition and taxon richness of 11 species of Neotropical ticks, collected from eight species of small mammals in 18 forest fragments across central Panama. Microbial richness varied among tick species, life stages, and collection sites, but was not related to host blood source. Microbiome composition was best explained by tick life stage, with bacterial assemblages of larvae being a subset of those of nymphs. Collection site explained most of the bacterial taxa with differential abundance across intrinsic and extrinsic factors. Francisella and Rickettsia were highly prevalent, but their proportional abundance differed greatly among tick species and we found both positive and negative co-occurrence between members of these two genera. Other tick endosymbionts (e.g. Coxiella, Rickettsiella) were associated with specific tick species. In addition, we detected Anaplasma and Bartonella in several tick species. Our results indicate that the microbial composition and richness of Neotropical ticks are principally related to intrinsic factors (tick species, life stage) and collection site. Taken together, our analysis informs how tick microbiomes are structured and can help anchor our understanding of tick microbiomes from tropical environments more broadly.

Tick sampling

Ticks were collected at 18 sites in central Panama between 2012 – 2014 (Figure 1, Table 1). All sites were lowland tropical rainforest, but varied in elevation from 30 to 460 m; (Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM); V2) and annual rainfall from 1825 to 3975 mm yr^-1 (48). The study sites also differed in disturbance histories and patch sizes, ranging from small forest fragments to large national parks, and therefore varied widely in wildlife community composition (49). Ticks were collected directly from small mammals during the dry season from January to May that were trapped using 100 Tomahawk live traps (40x13x13 cm) per site that were placed in a grid with 10 m interspacing and baited with ripe banana. Trapped mammals were identified to species using Reid (50) and subjected to extensive search over their entire body for ticks. We collected and analyzed ticks from 179 individual small mammals of eight species, including four species of opossums and four species of rodents (Table 1). All ticks collected were larvae or nymphs.

Ticks were stored in 95% ethanol and identified to species either morphologically or molecularly. Larvae and nymphs of Amblyomma ovale, Haemaphysalis juxtakochi, Ixodes affinis, and Ornithodoros puertoricensis were morphologically identified using an extensive reference collection at the Gorgas Memorial Institute and taxonomic keys (51). All other immature Amblyomma ticks were molecularly identified using either the 16S rRNA gene or the mtDNA COI barcode fragment (52). We also determined the DNA sequence of the 16S rRNA gene for a subset of the morphologically identified ticks to confirm correct species identification. In total, we identified 815 ticks (456 larvae and 359 nymphs) to 15 species: Amblyomma auricularium (n=82), A. dissimile (n=31), A. geayi (n=28), A. longirostre (n=1), A. mixtum (n=44) , A. naponense (n=5), A. oblongoguttatum (n=1), A. ovale (n=103), A. pacae (n=45), A. sabanerae (n=257), A. tapirellum (n=1), A. varium (n=9), H. juxtakochi (n=165), I. affinis (n=2), and O. puertoricensis (n=41). In addition, one Amblyomma larva yielded a sequence that did not match to any known species in GenBank.

DNA extraction and sequencing

We used molecular methods to assess tick microbiomes. All identified ticks (n=815) were individually processed under open benchtop conditions.  Prior to DNA extractions, we surface sterilized ticks in 10% bleach, then rinsed them with sterile molecular grade water, as described by Clay et al. (26). Genomic DNA was extracted from whole ticks using the Qiagen DNeasy kit (Qiagen, Valencia, CA, USA) following the manufacturer’s protocol. PCR products were obtained targeting the V1 to V3 hypervariable regions of the 16S rRNA gene using the forward primer 5’-ATTACCGCGGCTGCTG-3’ and reverse primer 5’-GTTTGATCCTGGCTCAG-3’, following Hawlena et al. (30). These hypervariable regions have previously been shown to be suitable for distinguishing bacterial species to the genus level (53). DNA extraction and PCRs control samples showed no isolation or significant amplification and were subsequently excluded. The 16S rRNA PCR products of the 815 samples were multiplexed and sequenced on an Illumina MiSeq sequencer. The dataset was deposited in the Short Read Archive (Bioproject XXXX).

Sequence and sample processing

The 16s rRNA sequences were demultiplexed and quality-filtered using Quantitative Insights Into Microbial Ecology (QIIME) 1.9 (54). All sequences were trimmed to 150 bp. Amplicon sequence variants (ASVs) — previously described as sub-operational taxonomic units (sOTU) — were determined using the deblur workflow (minimum q=19), which allows for unique taxa to be differentiated based on single-nucleotide differences (Deblur 1.1.0) (55). Bacterial taxonomic assignments were made using the Ribosomal Database Project Classifier (56).  All chloroplast and mitochondria sequences were removed and a phylogenetic tree was built from representative sequences in QIIME using fasttree2 (57). After filtering, the dataset comprised 10,689,269 high-quality 16S rRNA sequence reads from 733 tick samples of 11 species, with a median of 8,414 sequences per sample. ASVs with less than 10 sequence reads and/or less than 0.05 percent reads across all samples reads were removed (58). Samples were subsequently rarefied (without replacement using single_rarefaction.py in Qiime) at 2,014 sequences per sample to account for variable numbers of raw sequences. Analyses of alpha and beta diversity utilized the rarefied datasets, whereas correlation networks, analysis of composition of microbiomes (ANCOM) (59), and balance trees (60) were done on the unrarefied dataset, after any samples from libraries of sizes less than 1000 were excluded. In analysis of the unrarefied dataset, ASVs with low prevalence (total count less than half the number of samples) and presence (in less than 5% of samples) were also removed as described by McMurdie et al. (61).

Presumptively Rickettsia-positive samples were subjected to conventional PCR and DNA sequencing. Amplification of the rickettsial outer membrane protein A (ompA) gene was performed using primers R190-70 and R190-602 (62), modified to use GoTaq Green Master Mix (Promega, Madison, WI, USA) in 25 microliter reactions containing 1.0 M of each primer and 3 microliter of template DNA. Results of PCR were assessed by electrophoresis and UV-transillumination of GelStar (Lonza, Rockland, ME USA)-stained 1% agarose gels. Bands of the expected size were excised and cleaned with a QIAquick gel extraction kit (Qiagen) per manufacturer instructions. Products were sequenced in the forward direction in an ABI Prism 3730 Genetic Analyzer (UC DNA Sequencing Facility, Davis, California, USA). DNA sequences were manually trimmed and corrected if the nucleotide could be unambiguously determined, and then compared with sequences in a large database (GenBank, NCBI, Bethesda, MD, USA) by BLAST search.

Statistical analyses

We characterized within-community (α) ASV diversity (richness) and between-community (β) diversity of tick microbiomes in R version 3.3.3 (63). We included tick species, tick life stage, host blood source, and collection site as predictors for tick microbial diversity and composition. We first summarized the taxonomic composition of all bacteria at the genus level by each of these predictors (Figure 2). We also summarized the composition of bacterial genera that harbor well-known tick endosymbionts and/or pathogenic members (i.e. Anaplasma, Bartonella, Coxiella, Francisella, Rickettsia, and Rickettsiella) by each of these predictor variables (Figure 3). Additionally, we summarized the bacterial phyla present in each individual tick (Supplemental Figure 1).

To assess which factors best explain ASV richness, we performed a Generalized Linear Mixed Model (GLMM) with a Gaussian distribution and log link function using the glmmADMB package (64). Tick life stage was modeled as a fixed factor, while host species, tick species, and collection site were modeled as random factors. We used likelihood-ratio tests to test for significance of the random factors.

We used balance trees to investigate the groups of bacteria that change with intrinsic and extrinsic factors (60). Specifically, we ran a multivariate regression on balances to test the relative importance of each factor. Balance trees naturally correct for differences in sequencing dept, and use nonoverlapping subcommunities to account for the compositional nature of the data (65). Accounting for compositional data is not possible with previously common methods for analyzing tick microbial composition (e.g. PERMANOVA). Thus, for ease of comparison with previously published tick microbiome studies, we provide the results of PERMANOVA analyses (999 permutations) in the supplement (Supplemental Discussion, Supplemental Table S1, Supplemental Figure 7).

Analysis of composition of microbiomes (ANCOM) was used to investigate ASVs that were significantly over- or under-represented in each tick life stage, species, host order and collection site (59). Because machine learning within ANCOM requires greater that 25 samples for sufficient data training and testing, two tick species with fewer than 25 samples (A. varium and A. naponense) were removed prior to ANCOM analysis. Additionally, host species were summarized at the order level (i.e. Rodentia or Didelphimorphia) because several host species had fewer than 25 samples.

To test for pairs of ASVs that occurred together significantly more or less often than expected by chance, we performed a Sparse Correlations for Compositional Data (SparCC) network analysis using PySurvey (66). Using a Mantel test, we examined spatial autocorrelation of bacterial ASV by collection site.  Bipartite interaction networks were constructed using the R package ‘bipartite’ (67) to visualize notable bacterial genera (i.e. Anaplasma, Bartonella, Coxiella, Francisella, Rickettsia, Rickettsiella) that were detected per tick species and host order (Didelphimorphia and Rodentia).