Integrating population genetic structure, microbiome, and pathogens presence data in Dermacentor variabilis

Introduction

As vector-borne diseases continue to emerge and re-emerge in several regions of the world, there is an urgent need for a better understanding of their ecology, including the tripartite pathogen-vector-host relationship, to ultimately develop effective control strategies. Ticks represent an increasing threat to One Health given the range expansion of some species, and the frequent identification of new tick-borne pathogens (Paddock et al., 2016; Eisen et al., 2017). Due to the complex nature of diseases transmitted by ticks, their control requires interdisciplinary studies and collaboration (Murgia et al., 2019).

In recent years, the microbiome (referring to bacteria only throughout the manuscript) of vectors (mostly mosquitoes and ticks) has received much attention. The main reason is that associations between native microbes and pathogens may provide a new promising path towards the disruption of pathogen transmission (Narasimhan & Fikrig, 2015; Bonnet et al., 2017). In the case of ticks, microbes can be acquired from the environment (e.g., through the spiracles), horizontally (from the host, during blood feeding), or vertically (maternally inherited). Relationships between the tick and a particular microbe or set of microbes can range from mutualistic to parasitic (Casadevall, Fang & Pirofski, 2011; Bonnet et al., 2017). Microbes can affect tick fitness, vector competence, and pathogen transmission (Bonnet et al., 2017; Budachetri et al., 2018). For example, a few maternally inherited microbes have been proposed as primary symbionts necessary for tick success by providing vitamins that are lacking in blood (Smith et al., 2015b; Duron et al., 2017; Guizzo et al., 2017). Without these microbes, the ticks’ fitness is negatively impacted (Zhong, Jasinskas & Barbour, 2007). Other symbionts interact with pathogenic microorganisms, either positively or negatively. High proportions of the endosymbiotic Rickettsia bellii inhibit the transmission of the pathogenic Anaplasma marginale (Gall et al., 2016), and R. parkeri is likely excluded from A. maculatum by Candidatus Rickettsia andeanae (Paddock et al., 2015), presumably through a phenomenon known as interference. On the other hand, cases of facilitation have also been reported; such as a positive relationship between the proportion of the microbiome occupied by Francisella endosymbiont (FLE) and the infection level of F. novicida (Gall et al., 2016).

Even though a considerable amount of information and knowledge has been generated during the last decade or so, microbiome research on ticks is still in its infancy. For example, we still do not fully understand how host genetics and environmental factors interact to shape the microbiome of organisms (Spor, Koren & Ley, 2011; Goodrich et al., 2014; Steury et al., 2019), or how pathogenic microorganisms affect the microbiome and vice versa. Integration of different lines of evidence may be the key to improve our understanding of TBDs ecology. As highlighted by Griffiths et al. (2018), research exploring the links between, for example, vector genetics, microbiome composition and structure, and pathogen susceptibility may enable a better understanding of the factors governing disease in vulnerable populations. Unfortunately, to the best of our knowledge, there are no such studies on ticks. There are a few studies based on vertebrate taxa that have attempted to incorporate genomic and microbial data with data on environmental variation (Fietz et al., 2018; Griffiths et al., 2018; Steury et al., 2019). Overall, these studies have shown that genetically divergent host populations, exhibited more divergent microbiomes (Smith et al., 2015a; Griffiths et al., 2018; Steury et al., 2019). In amphibians, Griffiths et al. (2018) found that the genetic distance among hosts was correlated with microbial community dissimilarity when controlling for geographic distance. Nevertheless, the same study showed that at the site-level the microbiome did not mirror the host population genetic structure. Steury et al. (2019) investigated populations of Threespine Stickleback and concluded that the microbiome composition was better predicted by fish population genetic divergence than by geographic distance and environment. It is worth noting that the global trends appear to be driven by a subset of the microbiome. In other words, the influence of host genetic factors on the microbiome composition depends on the bacterial taxa in question; host genetics may affect the presence of some microbes, but presence of others may be better explained by, for example, the environment (Fietz et al., 2018). In the same way, differences between populations could be the result of differences in the relative abundance of a small subset of microorganisms (Steury et al., 2019).

The vertebrate results suggest that integrating microbiome and population genetic data may lead to better understanding of the ecology of tick-borne diseases. Do ticks from different genetic clusters harbor a different microbiome? Do infected ticks belong to a specific genetic cluster? Is the microbiome of infected ticks different from that of non-infected ticks? To start answering these questions, it is paramount to first generate high quality data, and then to develop an appropriated framework to integrate the different lines of evidence. This exploratory study focuses on Dermacentor variabilis (Say), a North American tick species that commonly bites humans, and that historically has been implicated in the transmission of several pathogens (Hecht et al., 2019). Dermacentor variabilis is also one of the most widely distributed ticks in the United States (Eisen et al., 2017), one for which population genetic data was already available. During this investigation, we first generated high-quality microbiome data and determined the presence/absence of several pathogens. Next, we combined those data sets with an existing population genetics data set for D. variabilis (see Lado et al., 2019). All three pieces of data (microbiome, pathogens, and population genetics) were derived from each of the individual ticks included in this study. The goal of this preliminary study was to integrate the different lines of evidence for D. variabilis ticks using clustering methods and multivariate statistics to identify trends and patterns of variation. Our expectation was that D. variabilis ticks that are genetically more alike (i.e., that belong to the same genetic cluster), will also be more alike at the microbiome level. Additionally, we hypothesized that infected and non-infected ticks would have a subset of microbes in their microbiome that are characteristic for each category.

Methods

Tick samples

The tick specimens employed in this study are the same individual ticks from Lado et al. (2019), with the exception of two specimens for which we did not have enough genetic material. The sample includes 64 adult D. variabilis collected from California (n = 3), Georgia (n = 2), Indiana (n = 3), Maryland (n = 4), Maine (n = 3), Minnesota (n = 5), North Dakota (n = 1), Ohio (n = 31), Oklahoma (n = 5), Tennessee (n = 5), Virginia (n = 1), and Washington (n = 1) (Fig. 1). All specimens were wild caught by dragging a 1m ×0.8m cloth thorugh the vegetation, and preserved in 95% ethanol. We collected ticks from vegetation (unfed) to minimize the amount of genetic material from the host. Collection data for the specimens used in this study can be accessed through the Ohio State Acarology Collection (OSAL) online database (https://acarology.osu.edu). Metro Parks permitted collection of tick species throughout the Metro Parks. The Division of State Parks (DNR), Indiana Department of Natural Resources authorized tick collection at Brown County State Park, Clifty Falls State Park, O’Bannon Woods State Park, Monroe Lake, Splinter Ridge Fish and Wildlife Area, Clark State Forest, and Ferdinand State Forest.

Results

Microbiome

Microbiome data were successfully generated for all ticks. However, some issues were encountered when attempting to merge forward and reverse reads (only a minimal fraction of the reads merged correctly). Because of that and given the large amount of data obtained, only forward reads (R1) were employed for all downstream analysis. It is possible that the difficulties encountered when merging the reads is due to the primers pair selected, as other research groups have faced similar issues with these primers for different taxa (Argonne National Lab personal communication). The number of Illumina forward (R1) reads obtained was 2,131,680. After quality filtering the number of reads retained was 2,129,331 with an average length of 251 bp. The average number of sequences per sample was 32,795 (minimum = 9,331; maximum = 44,533). The number of reads for the extraction blank was 2349, considerably lower than all samples; and no reads were generated from the library blank.

A total of 41 orders of bacteria were identified (13 remaining after filtering), and three of them were dominant: Rickettsiales, Legionellales, and Enterobacteriales (Fig. 2A). At the genus level, 61 taxa were detected (File S1). Five taxa were removed since their relative abundance was greater in the negative control than in tick samples. Of the 56 genera remaining, 36 were eliminated during the filtering steps according to minimum relative abundance and the number of samples (see Methods section for details). Thus, 20 genera were kept for downstream analyses (File S1). Three taxa were highly abundant at the genus level: Rickettsia spp., Francisella spp., and Enterobacteriaceae “other”. (Fig. 2B). Francisella spp. were present in all ticks, with a relative abundance ranging from 0.4 to 100% (mean = 81%, median = 98.5%). Rickettsia was detected in 12 (18.8%) of the ticks. In all cases where Rickettsia spp. had a relative abundance of at least 0.5% by microbiome analysis, tick specimens were also PCR positive, thus microbiome and PCR results are overall, congruent (Table 2). Rickettsiella, an intracellular endosymbiont of arthropods (Leclerque, 2008), was found in a tick from Ohio (1/64; tick ID 119572D) with a relative abundance of 11%.

Statistical analysis: integrating individual data sets

Clustering approach

The dendrogram resulting from hierarchical clustering showed that the two most distinct samples (clustering together) were from the same location and collection event: Fernald Preserve, Hamilton Co., Ohio (110559-1 and 110559-2) (Fig. 3). The second clade grouped two samples from CA, the two samples collected in Lake Co (115102A,B). The next pair of samples corresponded to the two remaining samples from the west coast, a sample from Napa Co, CA (115101C); and one from Whitman Co, WA (115105B). These two samples were closely associated with a clade that included all Northern samples. Lastly, all Eastern samples clustered together, without any internal grouping according to geography (Fig. 3).

No differences were evident upon visualization of community membership or structure. PCA did not show any specific clustering pattern between samples (File S2). In fact, most of the ticks were clustered in a single aggregate. The ticks that clustered somewhat separated were firstly the two that also appeared separated in the clustering analysis (110559-1 and 110559-2); and a tick from Ohio, 119241-3. It is important to note that the first two principal components explained only approximately 40% of the variance present in our sample.

As for k , the optimal number of clusters, silhouette method and elbow method suggested k = 2 and k = 4 respectively (File S3). Both values were close to the number of genetic clusters, which was also the number of clusters we used in the analysis. The C-index between the predicted clusters and genetic clusters was 0.61, indicating some weak but nonnegligible consistency between these two pairs.

Ordination methods

When visualizing the samples in the space through NMDS (stress value 0.064), the samples were spread out without any specific pattern (Fig. 4A). In PCoA with Gower’s distance as metric, most samples appeared in a single tight aggregate, while a few others dispersed randomly along the axes (Fig. 4B).

Test for significance for different variables.

The number of bacterial genera found in ticks belonging to the Eastern genetic cluster was higher than that of the Western, and Northern clusters. The genus Francisella was dominant in almost all the ticks, regardless of geography. Its relative abundance did not differ significantly between males and females (p-val 0.8); nor between ticks belonging to different genetic clusters (Eastern vs Western p-val 0.57, Eastern vs Northern p-val 0.96, Northern vs Western p-val 0.60) (Fig. 5). The relative abundance of this genus was also not significantly different between infected and uninfected ticks (p-val 0.14), although it was generally lower in infected ticks. Infected ticks were defined as those that were PCR positive for either Anaplasma or Rickettsia species.

Rickettsia relative abundance did not significantly differ between males and females (p-val 0.44), nor between ticks belonging to different genetic clusters (Northern vs. Eastern p-val = 0.26; Northern vs. Western p-val = 0.72; and Eastern vs. Western p-val = 0.45) (Fig. 5). Rickettsia rhipicephali and Anaplasma spp. were present only in ticks belonging to the Eastern genetic cluster, and R. montanensis was present in ticks belonging to the Northern and Eastern clusters. Rickettsia bellii was found infecting ticks from the Western and Northern genetic clusters, although of these ticks were collected in the Western region (Fig. 6). One of the ticks positive for this species was 115101C, a tick collected in California but with the genetic make-up corresponding to Northern ticks.

Overall, microbial beta diversity differed between regions (p-val < 0.05), but not between sexes (p = 0.63) when using Bray–Curtis distance. Post hoc tests indicated that ticks collected in the Western region harbored a different microbiome from those collected in the Eastern or Northern regions (p < 0.05). In the case of presence/absence of Bacteria (Jaccard distance), the results were the same. A significant difference between ticks collected from different regions was detected (p < 0.05), and the post hoc tests showed that the differences were significant between the Western and both the Eastern and Northern regions.

The Adonis test for genetic clusters showed that there were differences (p < 0.05) only when considering relative abundance of bacterial taxa (Bray–Curtis distance). When the pairwise Adonis test was performed to determine what genetic clusters differed in microbiome structure (p < 0.05) between Eastern and Western genetic clusters, but its correction was not (plFDR 0.096). Eastern and Western clades are genetically, the two more distant (Lado et al., 2019).

Discussion

This work corresponds, to the best of our knowledge, to the first attempt in tick research to integrate population genetic, microbiome, and pathogens presence data to better understand the ecology of TBDs. To optimize comparisons, it was our goal to generate those three pieces of data for each of the individual ticks included in the data set. We show that this is possible, even when techniques based on next generation sequencing require high quantities of high quality DNA.

We developed an effective pre-processing and processing procedure for researchers interested in microbiomes of small organisms, or parts of organisms with limited amounts of genetic material. This was motivated by a need for consensus regarding the pre-processing of ticks for generating microbiome data. For example, it is likely that many of the incongruencies between microbiome studies arise from differences in approaches to decontamination of the ticks’ surface. Several studies in the literature performed surface sterilization of the ticks before DNA extraction (e.g., Lado et al., 2018; Trout Fryxell & DeBruyn, 2016), whereas others did not (e.g., Clow et al., 2018); and the differences in bacterial communities may be a reflection of extensive environmental “contamination” in the latter. Even between studies that “washed” the ticks, comparisons need to be done cautiously, since not all sterilization methods are equally effective. The incorporation of bleach appears necessary for effective decontamination (Binetruy et al., 2019). Increased consistency in results of microbiome studies of ticks requires some community consensus related to the generation and analysis of microbiome data. This is especially important for non-model organisms such as ticks, given the “noise” (Alpha and Beta diversity inflation) generated by the presence of environmental microorganisms.

Consistent with most other studies of tick microbiomes (Van Treuren et al., 2015; Gall et al., 2016; Chicana et al., 2019), our analysis showed that the microbiome of many species is heavily dominated by a few genera, while the microbiomes of others, such as Ixodes angustus (Chicana et al., 2019) and Haemaphysalis lemuris (Lado et al., 2018), are more diverse. The microbiome of D. variabilis is dominated by Francisella (Chicana et al., 2019; Travanty et al., 2019; current study). Clow et al. (2018) presented, to some extent, opposing results (far higher microbial diversity), but they analyzed only 9 specimens, and the ticks were not surface sterilized. These results casts some doubt on the concept of a “core microbiome”. For example, Chicana et al. (2019) report a D. variabilis core microbiome as composed by Francisella, Sphingomonas, and Methylobacterium. On the other hand, Travanty et al. (2019) reported a D. variabilis core microbiome that includes Francisella spp., Sphingomonas spp., Delftia spp., and Hymenobacter spp. Our study recovered all of these taxa but only Fransicella spp. is nearly universally present. Francisella has been well established as endosymbiont and dominant in D. variabilis, as well as other Dermacentor species. All other taxa reported as “core microbiome” are not consistent across studies (e.g., Chicana et al., 2019; Clow et al., 2018; Rynkiewicz et al., 2015; Travanty et al., 2019). Thus, the true “core” of taxa overlapping across studies, can be reduced to Francisella only. The concept of “core microbiome” is generally used in microbiome research to refer to a suite of microbes, and not to refer to only one taxon. The utility of the concept of “core microbiome” for D. variabilis is therefore unclear.

Recent studies focusing on the pathogen transmission by D. variabis and on its role as vector of human diseases have consistently reported low prevalence of infections with known pathogenic microorganisms. For example, R. rickettsii, is usually present in ∼1% of the ticks analyzed e.g., (Hecht et al., 2019). The results obtained in the present study are consistent with these literature results, and show an overall Rickettsia spp. prevalence of 17.2%, including R. montanensis, R. bellii, and R. amblyommatis infections. Rickettsia rickettsii DNA was not found in any of the ticks analyzed. We detected DNA of Anaplasma spp. in one tick sample collected in Oklahoma. Its DNA sequence matched 100% with GenBank sequences submitted by researchers at the CDC Fort Collins, who isolated the agent from a human blood sample. No additional information is available in the literature about this case. This study comprises the first report of that specific bacterial agent in D. variabilis. However, an Anaplasma spp. bovis-like agent has been previously reported in D. andersoni ticks from Canada (Dergousoff & Chilton, 2011; Chilton, Dergousoff & Lysyk, 2018). Unfortunately, the region of 16S amplified differs between our study and that of D. andersoni (Dergousoff & Chilton, 2011), making it difficult to determine if the sequences correspond to the same specific agent. Regardless, it appears that uncharacterized Anaplasma agents are circulating in nature (ticks and humans), and their characterization, together with the determination of whether or not they are capable of causing disease in humans should be further explored. Furthermore, due to the finding of this agent’s DNA in a D. variabilis sample in this study, the role of D. variabilis as potential vector should be considered and further explored.

With respect to the primary goal of this preliminary study, the integration of the different types of data was done successfully, although the interpretation of the results is somewhat challenging. The latter problem is most likely caused by to the low levels of genetic diversity and moderate levels of population structure (see details in Lado et al., 2019), and a microbiome highly dominated by a few taxa. As a result, ordination methods failed to show patterns of variation across different groups. Nonetheless, the dendrogram resulting from the clustering analysis, in which all lines of evidence and geographic location were considered, was largely congruent with the observations at the population genetic level. Thus, the addition of microbiome, pathogen presence, ticks sex, and geographic location led to conclusions that are, in general, consistent with ticks genetics. All Eastern samples clustered together (with the exception of two samples from Hamilton Co, OH), and separated from both Western and Northern samples. The distinctiveness and separation of those two samples from Hamilton Co, OH reflects their distinctiveness at the microbiome level (Fig. 2). The fact that both samples appeared very similar to all other Eastern samples at the population genetic level (Lado et al., 2019), supports that hypothesis. It is possible that these two specimens were not well surface sterilized, or that the environment at that collection location is different.

The higher number of bacterial genera found in ticks belonging to the Eastern genetic cluster when compared to either Northern or Western ticks, is likely a product of the higher number of ticks analyzed, and the greater geographic area covered. Even though in the broader sense the microbiome is not diverse (dominated by three genera), it differed between the three main geographic regions: Eastern, Western, and Northern. Statistical results looking at the beta diversity between genetic clusters are less clear: initial test demonstrated a difference in the microbial communities between genetic clusters, but pairwise tests with their corresponding corrections failed to reach the same result. This may be a result of one of the study’s limitations, the number of samples. Thus, these findings should be further explored using more samples, and more samples per genetic cluster. From our preliminary analyses, geography explains the differences in the microbial communities better than host genetics. This could be the result of certain microorganisms occurring only in certain geographic areas.

Despite challenges arising during the interpretation of the results, the integration of lines of evidence and metadata, revealed that Eastern ticks can be separated from Northern and Western ticks. It can also be noted that some Rickettsia species were associated with a certain geographic area. For example, R. bellii was found only in samples collected in CA, a finding consistent with that reported in Hecht et al. (2019). It seems possible that R. bellii is more common along the west coast, although more samples should be analyzed to confirm this. And while the two R. bellii positive tick samples were collected in CA, one of them (ID 115101C) did not belong to the Western genetic cluster; it belonged to the Northern cluster. From the host genetics perspective, this tick is more similar to ticks from Northern locations; whereas from the rickettsial agents perspective, it is more similar to Western samples. Looking at the microbiome composition, this tick could belong to any geographic region. This underscores the value of integrating different types of data when thinking about disease ecology.

The approaches taken during this study, both to generate and analyze data, can be applied to a wide variety of taxa, and groups of organisms. It is likely that in organisms with a higher level of variation, the clustering methods will be more informative, and their interpretation more straightforward. Therefore, the approaches used herein have potential, and could greatly improve future studies looking at different aspects of diseases ecology. In the particular case of D. variabilis, the approaches used herein failed to detect clear tendencies or patterns. This was likely due to the low levels of variation (both in the microbiome and population genetics), and the very low prevalence of pathogenic microorganisms.

Conclusions

In conclusion, this is the first study that successfully generated microbiome, population genetics, and pathogens presence data for the same individual ticks. General methodologies and pre-processing steps are replicable, and applicable to different groups of organisms across the tree of life. This work also comprises, at present, one of the few studies aiming at integrating population genetics and microbiome data to better understand ecological processes and disease, and it is the first one to do so for ticks. The integration of different lines of evidence allows a more holistic approach; and clustering and ordination methods are very helpful to summarize and visualize the results. Finally, this study comprises the first report of “Uncultured Anaplasma sp. clone 15-3642” in D. variabilis. This agent was previously isolated from human blood, and it is important to determine if it is capable of causing disease in humans, and if so, the role of Dermacentor ticks as potential vectors.

Supplemental Information