DOI: 10.5061/dryad.h9w0vt4h0
title: Tick microbiomes in neotropical forest fragments are best explained by tick-associated and environmental factors rather than host blood source

description: 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.
description: <strong><em>Tick
 sampling</em></strong> 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
 (<i>ASTER</i>) Global Digital Elevation Model
 (<i>GDEM</i>); V2) and annual rainfall from 1825 to 3975 mm
 yr<sup>-1</sup> (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
 <i>Amblyomma ovale</i>,<i> Haemaphysalis
 juxtakochi</i>, <i>Ixodes affinis</i>, and
 <i>Ornithodoros puertoricensis</i> were morphologically
 identified using an extensive reference collection at the Gorgas Memorial
 Institute and taxonomic keys (51). All other immature <i>Amblyomma
 </i>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:
 <i>Amblyomma auricularium</i> (<i>n</i>=82),
 <i>A. dissimile </i>(<i>n</i>=31), <i>A.
 geayi </i>(<i>n</i>=28), <i>A.
 longirostre</i> (<i>n</i>=1), <i>A. mixtum
 </i>(<i>n</i>=44) , <i>A. naponense
 </i>(<i>n</i>=5), <i>A. oblongoguttatum</i>
 (<i>n</i>=1), <i>A. ovale</i>
 (<i>n</i>=103), <i>A. pacae
 </i>(<i>n</i>=45), <i>A. sabanerae
 </i>(<i>n</i>=257), <i>A. tapirellum</i>
 (<i>n</i>=1), <i>A. varium
 </i>(<i>n</i>=9),<i> H. juxtakochi
 </i>(<i>n</i>=165), <i>I. affinis
 </i>(<i>n</i>=2), and <i>O. puertoricensis
 </i>(<i>n</i>=41). In addition, one <i>Amblyomma
 </i>larva yielded a sequence that did not match to any known species
 in GenBank.  
 <strong><em>DNA extraction and
 sequencing</em></strong> We used molecular
 methods to assess tick microbiomes. All identified ticks
 (<i>n</i>=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).  
 <em><strong>Sequence and sample
 processing</strong></em> 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
 <i>Rickettsia</i>-positive samples were subjected to
 conventional PCR and DNA sequencing. Amplification of the rickettsial
 outer membrane protein A (<i>ompA</i>) 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.  
 <em><strong>Statistical
 analyses</strong></em> 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. <i>Anaplasma</i>,
 <i>Bartonella</i>,<i> Coxiella</i>,<i>
 Francisella</i>, <i>Rickettsia</i>, and
 <i>Rickettsiella</i>) 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 <i>glmmADMB</i> 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 (<i>A. varium</i> and <i>A. naponense</i>)
 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.
 <i>Anaplasma</i>, <i>Bartonella,
 Coxiella</i>,<i> Francisella</i>, <i>Rickettsia,
 Rickettsiella</i>) that were detected per tick species and host
 order (Didelphimorphia and Rodentia).