Drivers and patterns of microbial community assembly in a Lyme disease vector

Abstract Vector‐borne diseases constitute a major global health burden and are increasing in geographic range and prevalence. Mounting evidence has demonstrated that the vector microbiome can impact pathogen dynamics, making the microbiome a focal point in vector‐borne disease ecology. However, efforts to generalize preliminary findings across studies and systems and translate these findings into disease control strategies are hindered by a lack of fundamental understanding of the processes shaping the vector microbiome and the interactions therein. Here, we use 16S rRNA sequencing and apply a community ecology framework to analyze microbiome community assembly and interactions in Ixodes pacificus, the Lyme disease vector in the western United States. We find that vertical transmission routes drive population‐level patterns in I. pacificus microbial diversity and composition, but that microbial function and overall abundance do not vary over time or between clutches. Further, we find that the I. pacificus microbiome is not strongly structured based on competition but assembles nonrandomly, potentially due to vector‐specific filtering processes which largely eliminate all but the dominant endosymbiont, Rickettsia. At the scale of the individual I. pacificus, we find support for a highly limited internal microbial community, and hypothesize that the tick endosymbiont may be the most important component of the vector microbiome in influencing pathogen dynamics.


| INTRODUC TI ON
Tick-borne diseases pose a serious and increasing threat to global human and animal health as tick distributions expand and new tick-borne pathogens are identified (Eisen, Kugeler, Eisen, Beard, & Paddock, 2017;Paddock, Lane, Staples, & Labruna, 2016).
It is evident that tick microbes can be maternally inherited (vertical transmission), environmentally acquired through the tick spiracles, mouth, or anal pore (environmental transmission), or obtained from host-blood feeding (horizontal transmission) (Narasimhan & Fikrig, 2015). The microbial community resulting from these processes consists largely of tick symbionts and guest commensals (Clay & Fuqua, 2010;Greay et al., 2018), which can affect vector fitness (reviewed in Bonnet et al., 2017) and have been associated with variation in vector competence (Budachetri et al., 2018;Civitello, Rynkiewicz, & Clay, 2010;Gall et al., 2016;Narasimhan et al., 2014;Telford, 2009). The latter observation motivated an increased number of studies of the role of ecological and environmental factors in shaping the microbiome and the relationship between microbial symbionts and pathogens (Fryxell & DeBruyn, 2016;Gall, Scoles, Magori, Mason, & Brayton, 2017;Kwan, Griggs, Chicana, Miller, & Swei, 2017;Narasimhan et al., 2014;Van Treuren et al., 2015;Zolnik, Prill, Falco, Daniels, & Kolokotronis, 2016). These reports implicate a complex suite of factors affecting the tick microbiome including geography, sex, and life stage, and provide opposing evidence about the effect of native tick microbiota on pathogen acquisition, highlighting the inherent complexity of tick microbial ecology and obscuring generalizable patterns.
In particular, studies investigating interactions between and among tick endosymbionts and pathogens have found competitive, facultative, or nonexistent microbial associations driven by a variety of underlying mechanisms (reviewed in de la Fuente et al., 2017).
For example, numerous studies have found negative correlations between vertically transmitted endosymbionts and the colonization or transmission of tick-borne pathogens (Gall et al., 2016;Macaluso, Sonenshine, Ceraul, & Azad, 2002;Narasimhan et al., 2014;Steiner et al., 2008;Telford, 2009). However, the suspected underlying mechanisms have varied in each case, ranging from differences in tick gene expression to microclimate and habitat-related factors.
Other studies have found coinfection rates among various tick-associated microbes to exceed random expectation (Civitello et al., 2010;Mather, Ribeiro, & Spielman, 1987;Zeidner et al., 2000), suggesting mutual facilitation. Still others have found no difference in the microbiome composition of infected and uninfected ticks, suggesting minimal interaction between the native microbiota and pathogens Kwan et al., 2017).
In this study, we attempt to address this complexity by investigating the fundamental processes shaping the larval microbiome of Ixodes pacificus, the vector of the Lyme disease pathogen in the western United States. We apply a community ecology framework, using each tick microbiome as a community, to explore the assembly and interactions of tick-associated microbes. Using larval I. pacificus, placed in permeable bags in the field, we confine the inputs to the microbiome to environmental and vertical transmission routes and analyze the effects of these routes on microbial diversity, composition, loads, and function. Our results demonstrate a generally limited I. pacificus microbiome that is dominated by the vertically transmitted endosymbiont, Rickettsia, with populationlevel diversity and composition patterns driven by vertically acquired microbes and, to a lesser degree, environmentally acquired microbes.
F I G U R E 1 Experimental setup and main findings. Larvae from 3 maternal adults (clutches) were randomly assigned to environmental exposure treatments which included remaining in laboratory, or spending 2, 4, or 6 weeks buried in soil in oak woodland habitat. These vertical and environmental transmission routes generated variation in microbial diversity and composition across clutch and exposure groups, but not in predicted microbial function or microbial loads 2 | MATERIAL S AND ME THODS

| Sample collection
Adult female I. pacificus were field-collected from China Camp State Park in Marin County, CA, and fed to repletion on New Zealand white rabbits in the laboratory. Three engorged females were obtained from these means, and collectively produced 93 surviving I. pacificus larvae, hatched in the laboratory. These larvae were randomly assigned to spend 0, 2, 4, or 6 weeks in the field environment in an oak woodland forest in northern California (Figure 1). Ticks assigned to the 0-week group were maintained in the laboratory and, upon hatching, were surface-sterilized with hydrogen peroxide and ethanol, then flash-frozen with liquid nitrogen and stored in 70% ethanol. Ticks assigned to all other groups were placed inside sealed but permeable silk, mesh bags following the procedures of Padgett and Lane (2001). Separate mesh bags were used for ticks from each clutch and environmental exposure treatment to ensure we could later identify the treatment group of each tick sample. The silk, mesh bags were hung inside metal cages constructed from wire mesh that was 10 × 10 × 10 cm to prevent predation on ticks from insectivores and birds while still allowing the larvae to move vertically in the soil to avoid desiccation. All cages were buried with the tops of the tick bags exposed to the soil surface. All cages were buried at Pepperwood Preserve in Sonoma County in oak woodland habitat (see Appendix S1 for further details). After 2, 4, or 6 weeks, larvae were gathered, surface-sterilized, and flash-frozen on site, then stored in 70% ethanol.

| Microbiome sample preparation and sequencing
All ticks were thoroughly surface-sterilized with successive washes of hydrogen peroxide, ethanol, and de-ionized H 2 0 to remove environmental contamination (Rudolf et al., 2009). Whole ticks were individually pulverized using a sterilized pestle. Genomic DNA was then extracted using a Qiagen DNeasy Extraction Kit (Qiagen) following the manufacturers' specifications and using an elution volume of 100 μl. Libraries from individual ticks were then prepared for 16S sequencing following the Illumina MiSeq 16S Metagenomic Sequencing Library Preparation Protocol (Klindworth et al., 2013) with amplicon primers targeting the V3-V4 hypervariable region.
Each sample was amplified in triplicate to reduce PCR bias and then pooled for DNA purification using paramagnetic beads (Appendix S1). To enable differentiation postsequencing, purified amplicons were then barcoded using dual-index primers with Illumina adapters supplied in a Nextera XT Index Kit (Illumina). A combined library was then prepared by combining equimolar concentrations of all purified, barcoded samples. This combined library contained all 93 larvae, the 3 maternal adults, and 3 negative controls originating from the DNA extraction step (Appendix S1). The library was sequenced on an Illumina MiSeq using the V3 reagent cartridge (300 base pair, paired-end).

| Sequence analysis
Sequence reads were quality-filtered and processed using QIIME (for more details, see Appendix S1). In total, 84% of reads passed quality filter with an average of 44,995 reads per sample. Sequences were clustered at 97% sequence similarity and rarefied to a depth of 10,182 reads per sample (Appendix Figure S1) to correct for uneven sampling. Rarefying to this depth retained 68 of the original 93 tick samples (see Appendix Table S1 for sample size by treatment).
Reads were then assigned to operational taxonomic units (OTUs) using an open reference picking strategy with the NCBI taxonomic database (Benson, Karsch-Mizrachi, Lipman, Ostell, & Wheeler, 2006 (Glassing, Dowd, Galandiuk, Davis, & Chiodini, 2016). After these quality-filtering steps, a total of 23 genera, including the rare genera category, were retained for downstream analysis on the remaining 65 larvae and 3 adult samples. While the quality-filtering steps may have removed or pooled rare microbes with important functions (Jousset et al., 2017), we cannot distinguish these microbes from suspected contaminants, and find it more defensible to pool these genera than to leave in hundreds or thousands of dubious OTUs.

| Bacterial load quantification
To measure the overall abundance, or load, of microbiota present in the tick microbiome, we used a SYBR-based quantitative PCR on the 16S rRNA gene (Appendix S1) (Bacchetti De Gregoris, Aldred, Clare, & Burgess, 2011). Of the original 93 larval samples, only 71 had sufficient volume of DNA extract remaining for testing. The qPCR was performed on these 71 samples which spanned all treatment groups, and each sample was tested in triplicate with each replicate requiring a template volume of 7 μl (Appendix Table S5). Another qPCR protocol targeting the genes encoding outer membrane protein A was employed to quantify the load of the dominant I. pacificus endosymbiont, Rickettsia phylotype G021 (Appendix S1) (Cheng, Vigil, Schanes, Brown, & Zhong, 2013b). This qPCR was conducted on the 54 samples with remaining DNA extract, and each sample was tested in triplicate with each replicate requiring a template volume of 5 μl. These 54 samples did not include any larvae from clutch 2, as there was no DNA extract remaining from samples in this group.
Thus, analysis of Rickettsia abundance by clutch pertains only to clutch 1 and clutch 3.

| Community ecology analysis
We defined the microbiome of an individual tick as a community and applied standard community ecology analysis techniques to measure microbial presence, abundance, and composition within ticks and compare these metrics across ticks from different clutches (i.e., arising from different adult females) and field exposure times.
To characterize microbial alpha diversity within a tick, species richness and evenness were calculated manually using the fully quality-filtered dataset. Shannon's diversity, which is calculated from both richness and evenness of OTUs, was also calculated for all ticks using the vegan package in R (Dixon, 2003). To incorporate phylogenetic differences between species, we also calculated Faith's phylogenetic diversity using the R package picante (Kembel et al., 2010) (see Appendix Figure S10 for phylogenetic tree). Diversity values were then log-transformed to meet the normality assumption of the ANOVA used to detect differences between groups.
To compare tick microbiome composition across treatments, the distance between microbial communities was measured using two metrics targeting different community features (see Appendix S1 for further details). The Jaccard dissimilarity index was used to measure differences in microbial presence-absence, and the Bray-Curtis dissimilarity index was used to capture differences in microbial abundance. These pairwise dissimilarity indices were calculated for all possible pairs of ticks using the vegdist function from the vegan package in R. We then partitioned these dissimilarity matrices based on clutch and field exposure times through permutational multivariate analysis of variance (PERMANOVA) using the vegan function adonis, and obtained significance values through permutation of the raw data (n = 5,000). To determine where the differences occurred, post hoc tests were conducted using the pairwise.adonis function. All p-values were corrected for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) procedure (Benjamini & Hochberg, 1995).
To visualize these results, the weighted and unweighted dissimilarity indices were plotted using nonmetric multidimensional scaling (NMDS), an ordination technique which transforms highly dimensional data into a two-dimensional representation, implemented using the phyloseq package in R (McMurdie & Holmes, 2013).

| Core microbiome analysis
To determine which OTUs constituted the core I. pacificus microbiome, we used indicator species analysis using the indicspecies package in R (Cáceres & Legendre, 2009). This analysis distinguishes the genera specific to certain treatment groups (i.e., specific to a certain clutch or exposure time) from those shared among all I. pacificus microbiomes.
Indicator OTUs were those genera with a statistically significant association (p < 0.01, FDR corrected) to a particular clutch or exposure period. p-values were then generated through a permutation test using 1,000 permutations of the original data test, and were corrected for multiple testing using the Benjamini-Hochberg procedure. Core microbiota were defined as those OTUs occurring in all treatments, an association which cannot be statistically tested as there is no external group for comparison (Cáceres, 2013).

| Co-occurrence and network analysis
To test general assembly rules governing the tick microbiome and measure the types of interactions occurring between OTUs therein, we examined the co-occurrence patterns of the observed OTUs. We calculated checkerboard scores, or C-scores, which measure the average number of mutual exclusions for each pair of OTUs across a set of communities (Stone & Roberts, 1990). For example, the number of tick samples in which only Bacillus or only Pseudomonas occur (not both Bacillus and Pseudomonas) constitutes the number of mutual exclusions, or "checkerboard units," for this OTU pair. The number of checkerboard units is similarly measured for all possible pairs of OTUs, and averaged to determine the C-score for the set of communities. These calculations were performed using the bipartite package in R (Dormann, Gruber, & Fründ, 2008).
This value was compared to the average C-score calculated from 5,000 permutations of null communities preserving the marginal totals and connectance from the original data. This approach enabled comparisons of the degree of competitive structuring of the tick microbiome against a null hypothesis of random assortment (Stone & Roberts, 1990). Actual and simulated C-scores were calculated at the genus level and phylum level to ensure that observed outcomes were not an artifact of the taxonomic rank selected. To further examine microbial interactions, we calculated correlation coefficients for each pair of genera based on their distributions across tick samples. Positively correlated genera were pairs which tended to co-occur in tick samples, while negatively correlated genera were pairs which co-occurred less frequently than expected by chance.

| Functional role analysis
To estimate the functional role of the microbiota present in I. pacificus, the functional gene content was estimated using PICRUSt (v1.1) (Langille et al., 2013). This computational approach uses evolutionary modeling to predict metagenome function from 16S rRNA sequence data and a reference genome database. Briefly, the normalized OTU table was multiplied by the known or inferred gene function for each OTU. Gene functions were then classified into gene family predictions using the KEGG orthology database (Kanehisa, Goto, Sato, Furumichi, & Tanabe, 2012). Differences in gene family predictions between clutches and exposure times were measured via ANOVA using an alpha level of 0.01 and FDR correction. Using PICRUSt, gene predictions can be collapsed into various hierarchical levels; we compared gene predictions between treatment groups at two different hierarchical levels (Level 1 and Level 2) to validate results. The overall feasibility of this approach was assessed by calculating the weighted Nearest Sequenced Taxon Index (NSTI), a measure of the availability of nearby genome representatives for the given OTUs.

| Bacterial load quantification
The absolute abundance of microbes was measured for ticks from all treatment groups via qPCR. Total microbial loads for all groups ranged from 160,438 to 2,182,181 copies/µl, but no significant differences in loads were observed between ticks from different clutches or from different environmental exposure times (Appendix Figure S6). The abundance of Rickettsia phylotype G021, the dominant endosymbiont in I. pacificus, was also directly measured using a targeted qPCR.
Rickettsial loads ranged from 3,052 to 47,987 copies/µl, but again, no significant differences in loads were observed between clutches or between environmental exposure times (Appendix Figure S7).
F I G U R E 2 Ixodes pacificus microbiome composition after sequence quality filtering. Each of the 5 vertical bars represents an averaged microbiome for samples from that treatment. Each color represents a different OTU present in the microbiome, shown here at the genus level, and the heights of each bar represent the proportion of reads attributed to that OTU. The "Other" category represents all OTUs which did not meet the criteria of presence in at least one sample at ≥1% relative abundance. For similar community composition barplots with genera pooled at a higher level, see Appendix Figure

| Core microbiome analysis
Indicator species analysis was used to identify OTUs associated with ticks from specific clutches or times (indicator species) and those associated with ticks from all treatments (core microbiota).
Indicator species included microbes from the genera Fusobacterium,

| Co-occurrence and network analyses
Community assembly rules were tested by examining the co-occurrence patterns of OTUs across tick microbiome samples. The observed checkerboard score, an index measuring the number of species which never co-occur, was significantly lower than that of the simulated communities, indicating that there were fewer negative interactions among OTUs than expected by chance (Stone & Roberts, 1992

| Functional role analysis
PICRUSt, a functional gene content prediction approach, was utilized to estimate the functional role of the I. pacificus microbiome based on the 16S rRNA sequencing data. The average weighted Nearest Sequenced Taxon Index (NSTI) was 0.0134 (SD = 0.0004), indicating that our samples were highly tractable for metagenome prediction (Appendix S1) (Langille et al., 2013). Comparisons of the predicted gene pathways across clutches and exposure times yielded no significant differences between groups regardless of the hierarchical level at which the gene pathways were collapsed (Appendix Figure   S5). That is, the estimated functional gene content of the I. pacificus microbiome did not vary by treatment, despite differences in species composition and diversity.

| D ISCUSS I ON
We examined the process of community assembly in the initial, and those aforementioned is orders of magnitude greater than that typically generated by species-and life stage-specific differences.
Evidence suggests that the high diversity estimates previously reported are inflated due to contamination (Salter et al., 2014;, sequencing error (Huse, Welch, Morrison, & Sogin, 2010), or the presence of transient microbes (Wang, Gilbreath, Kukutla, Yan, & Xu, 2011;Zolnik, Falco, Daniels, & Kolokotronis, 2018). While the raw sequence data obtained in this study similarly suggested the presence of thousands of OTUs, we used rigorous statistical methods to control for potential contamination and sequencing error.
Further, we used high sample replication and sequencing results from multiple time points to better differentiate rare or transient microbes from common tick microbial residents. Using these methods, we found strong support for a highly limited internal I. pacificus microbiome.
Despite a limited microbiome within individual ticks, drivers and patterns of microbial diversity were discernible at the population level. Specifically, we found that maternal identity and environmental exposure period generated significant variation in microbiome diversity and composition. However, these vertical and environmental transmission routes did not generate differences in overall microbial abundance, endosymbiont abundance, or functional gene content.
These results indicate that microbial diversity and composition patterns do not necessarily reflect differences in microbial function and that abundance and function appear to be more highly conserved features of the I. pacificus microbiome.
We further found that while ticks can acquire internal microbiota from their environment, microbiome diversity decreased with increasing environmental exposure time. Additionally, microbial diversity was lower in environmentally exposed larvae compared to laboratory-maintained larvae, but composition did not differ between these groups. These results suggest that environmental transmission plays a relatively fleeting role in shaping the microbiome, in contrast to prior findings that environmentally acquired microbes contribute substantially to the tick microbiome (Carpi et al., 2011;Greay et al., 2018;Kwan et al., 2017;Zolnik et al., 2016).
To investigate the processes underlying the observed reduction in microbial diversity over time, we examined microbial interactions through co-occurrence analysis. We found fewer negative interactions than expected by chance, suggesting that the I. pacificus microbiome is not strongly shaped by microbial competition. Specifically, we found fewer microbial pairs with nonoverlapping distributions in the true community than in a null community model, and fewer negative pairwise correlations than positive correlations. These results suggest that the observed decline in I. pacificus microbial alpha diversity over time was not driven by competition among microbes but were potentially imposed by environmental conditions within the tick itself (Costello, Stagaman, Dethlefsen, Bohannan, & Relman, 2012;Faust et al., 2012). The corresponding decrease in phylogenetic diversity further suggests that conditions within the tick may shift over time allowing only certain more closely related taxa to persist in the environment (Webb, 2000;Webb, Ackerly, McPeek, & Donoghue, 2002). Environmental filtering, in which habitat conditions allow only taxa with particular traits or phenotypes to persist, is a common phenomenon in bacterial communities (Horner-Devine & Bohannan, 2006;Newton, Jones, Helmus, & McMahon, 2007).
While the process of environmental filtering was not directly explored here, evidence of the host niche dictating microbial assembly, similar to that seen here, has previously been observed across diverse study systems including C. elegans (Berg et al., 2016), fish (Schmidt, Smith, Melvin, & Amaral-Zettler, 2015), cacti (Fonseca-García et al., 2016), and humans (Levy & Borenstein, 2013). In ticks, F I G U R E 5 Microbial co-occurrence analysis patterns. (a) The distribution of C-scores, a measure of species co-occurrence, is shown here with OTUs aggregated at the genus level. As indicated by the arrow, the actual C-score is significantly lower than that of the simulated communities, indicating a nonrandom species assembly process. (b) Network analysis reveals a lack of microbial competition in the Ixodes pacificus microbiome. The interaction network shows the pairwise correlations between all OTUs (genera) found across all samples and treatments. Nodes represent genera found in the microbiome and are colored by phylum. Edges represent significant pairwise correlations between genera. Red edges denote negative correlations, and green edges denote positive correlations the microbial selection process may be driven by internal tick characteristics such as redox conditions and pH (Hyde, Trzeciakowski, & Skare, 2007;Narasimhan et al., 2014), low quantities of essential nutrients (Narasimhan & Fikrig, 2015), innate immune response used to control microbes (Hajdušek et al., 2013;Sonenshine & Roe, 2013), or microbe-microbe competition (Devevey, Dang, Graves, Murray, & Brisson, 2015;Gall et al., 2016). As tick microbes vary greatly in their contribution to vector fitness, we hypothesize that I. pacificus creates favorable conditions for microbes providing strong fitness benefits, perhaps at the expense of other microbes (Burgdorfer, 1988;Macaluso et al., 2002).
Identifying the core microbiome, those microbes present within all or the vast majority of assemblages (Turnbaugh & Gordon, 2009), allowed for further exploration of tick-microbe associations.
Although 22 OTUs were found across tick samples after rigorous sequence quality control measures, only Clostridium and Rickettsia were consistently found in ticks across clutches and time points (Appendix Table S3). The role of Clostridium within the tick remains unknown, but this common environmental bacterium has been reported in various Ixodid tick species around the world (Andreotti et al., 2011;Carpi et al., 2011;Clow, Weese, Rousseau, & Jardine, 2018 efficiency, to increase significantly after engorgement, and to play a nutritional role within ticks (Cheng, Vigil, et al., 2013b). Specifically, this Rickettsia phylotype has the demonstrated capability of synthesizing folate (Hunter et al., 2015), an essential vitamin absent from the ticks' blood-based diet and which ticks lack the genetic capacity to synthesize de novo (Hill & Wikel, 2005 and Enterobacteriaceae may limit pathogen colonization in Ixodes (Ross et al., 2018), these microbes are not reliably vertically transmitted, and their ability to persist in the tick during a molting event is unknown. Endosymbionts such as Rickettsia, however, are efficiently vertically transmitted and have been associated with reductions in pathogen acquisition in some tick-borne pathogen systems (Gall et al., 2016;Macaluso et al., 2002;Telford, 2009). However, the mechanisms by which Rickettsia, an intracellular microbe (Winkler, 1990), interact with pathogens remain unknown, and associations between Rickettsia and pathogen presence or loads are not consistently found . Our results found that, while negative microbial interactions were relatively rare within the I. pacificus microbiome, Rickettsia was involved in the majority of the negative interactions.
It remains to be determined whether Rickettsia is directly outcompeting these other microbes, or opportunistically or independently increasing as others decrease within the tick over time. Elucidating the role of Rickettsia, and tick endosymbionts generally, in tick-borne pathogen dynamics, constitutes a critical next step forward in this field.
Overall, we demonstrate the use of a community ecology framework for investigating features of the vector microbiome.
Characterizing the initial microbial community assembly patterns using thoroughly filtered sequence data to avoid spurious conclusions provided a framework for evaluating vector microbial variation and interactions. In this study, we demonstrate that vertical transmission routes contribute to population-level patterns of microbial diversity and composition within I. pacificus, with environmental transmission playing a minimal role. Despite discernible patterns at the population level, we found that individual ticks harbor lowdiversity microbiomes, increasingly dominated by the vertically transmitted endosymbiont, Rickettsia. Given the dominance of this endosymbiont as well as tick-associated endosymbionts in general (Narasimhan & Fikrig, 2015), future efforts to investigate tick microbiomes for applications in disease prevention may want to focus on the role of vector-associated endosymbionts in pathogen acquisition, maintenance, and transmission. In particular, investigating how vector competence varies with obligate endosymbiont burdens or facultative endosymbiont presence may identify microbes capable of disrupting pathogen transmission that may be harnessed in novel strategies to prevent tick-borne disease.

ACK N OWLED G M ENTS
We wish to thank Enxhi Tahiraj for laboratory assistance. This research was supported by NSF grants 1745411, and 1427772, the Bay Area Lyme Foundation, the Presidio Trust Foundation, and a CSUPERB New Investigator Grant to AS.

CO N FLI C T O F I NTE R E S T
None declared.

DATA ACCE SS I B I LIT Y
The project OTU