Four different datasets were utilized in this study. Therefore, for clarity, each dataset has been summarized in Table 1.
Table 1. All datasets used in this study.
Dataset (name)
|
Ticks (n)
|
Life stages
|
Locations in the Netherlandsa
|
Year of sampling
|
Molecular technique
|
Source
|
Microbiome
|
655 in 133 pools
|
Larvae, nymphs, adult females & males
|
AW, DK, ST, BU, HD, ZM (n = 6)
|
2016
|
16S rRNA sequencing
|
This study
|
Symbiont
|
16,555
|
Nymphs
|
AW, DK, ST, BU, HD, ZM, PD, SD, VH, VA, MH, HM, BB, PW, DW, EN, RB, VL, KB (n = 19)
|
2013-2014
|
qPCR targeting: R. helvetica, S. ixodetis, M. mitochondrii, Rickettsiella spp.
|
This study
|
Pathogen
|
13,967
|
Nymphs
|
AW, DK, ST, BU, HD, ZM, PD, SD, VH, VA, MH, HM, BB, PW, DW, EN, RB, VL, KB (n = 19)
|
2013-2014
|
qPCR targeting: B. burgdorferi s.l., A. phagocytophilum, N. mikurensis
|
Takumi et al., 2019
|
Transmission mode
|
1,130
|
Larvae, nymphs, adult females & males
|
AW, ST
(n = 2)
|
2019
|
qPCR targeting: R. helvetica, S. ixodetis, M. mitochondrii, Rickettsiella spp. B. burgdorferi s.l., A. phagocytophilum, N. mikurensis
|
This study
|
aThe abbreviations of locations are explained in Additional File 1: Table S1.
Microbiome dataset
Sample collection and preparation for microbiome analyses of I. ricinus
We utilized pools of I. ricinus larvae, nymphs, and individual adult females and males from six locations across the Netherlands for microbiome profiling (Table 1; Fig. 1, triangles). These locations were selected based on pre-existing knowledge of B. burgdorferi s.l. prevalence, the density of ticks, vegetation profile, and vertebrate community obtained in a cross-sectional study (22). The full names of the sites, their coordinates, and vegetation descriptions are provided in Additional file 1: Table S1.
Pooled and individual ticks from six forest sites (triangles) were used for a 16s rRNA amplicon sequencing analysis. Individual nymphs from these and 13 (points) other forest sites were tested by qPCR for the presence of tick symbionts. A box marks the sampling site by two letters, and a linear colour gradient represents latitude. Full coordinates, habitat, vegetation cover, tick densities, and a number of vertebrate species per locations are provided in Additional file 1: Table S1.
Questing I. ricinus were collected in 2016 by blanket dragging. In total, 655 ticks were combined into pools by life stage and location. All ticks were washed three times in 70% ethanol, and DNA from 32 pools of 10 larvae, 18 pools of 5 nymphs, 18 pools of 10 nymphs, 34 individual female, and 31 individual male ticks was extracted using the Qiagen DNeasy Blood & Tissue Kit according to the manufacturer’s protocol (Qiagen, Venlo, The Netherlands). A sampling scheme and sample metadata are provided in Additional file 1: Table S2 and Table S3.
Microbial profiling and taxonomic clustering
Illumina MiSeq V3-V4 region of 16S rRNA amplicon libraries were generated and sequenced by BaseClear (Leiden, the Netherlands). The description of the method has been published previously (30). Sequenced reads were imported to CLC Genomics Workbench 10.0.1 supplemented with CLC Microbial Genomics Module 3.6.1 (www.clcbio.com). Overlapping pairs of raw reads were merged into single longer reads and trimmed with a quality score limit of 0.05 and 2 ambiguous nucleotides. At this stage, primer sequences were trimmed. Subsequently, reads were fixed-length trimmed (~ 400 bp). To identify operational taxonomic units (OTUs), reads were clustered using the reference databases SILVA 16S version 128 with 97% identity as the clustering criterion; chimeras were removed.
16S rRNA quantification and total bacterial load
Total bacterial load in all samples was quantified, and proportions were multiplied by the load to convert relative into absolute abundances. Quantification of total bacterial DNA load was determined by 16S rRNA qPCR as previously described (31). The total bacterial load is a cost-effective and scalable solution for datasets of this size, since quantification methods through flow cytometry are not compatible with the sampling technique. Samples were normalized to control for arbitrary variation in sequencing depth, and the normalized abundances were scaled by the 16S rRNA qPCR values of each sample. For diversity analyses, samples were rarefied to the lowest sequencing depth.
Microbiome analyses
All analyses were carried out in R 3.6.0 (32). We used the R package vegan (version 2.5-6) for ordinations, diversity indices, and fitting of environmental vectors or factors onto ordinations (33). We also computed silhouette scores from the Bray-Curtis dissimilarities and Sheldon evenness (34). All principal coordinate analyses were carried out using Bray-Curtis dissimilarities, and envfit correlation to the principal components was corrected for multiple-testing with the Benjamini-Hochberg correction. Additionally, we carried out PERMANOVA with the adonis function from R package vegan to assess whether tick life stage significantly affected community variation. We tested for multivariate spread through the betadisper function.
To test how well different factors explained clusters observed on the PCoA plots, we calculated the Bray-Curtis dissimilarities and evaluated cluster quality as the silhouette score with the factors as cluster labels. The silhouette score, bounded by -1 and 1, takes both cluster cohesion and cluster separation into account. We used k-means clustering of the log-transformed Rickettsia abundances to compute silhouette scores for Rickettsia.
To investigate correlations between OTUs and tick life stage, we fitted proportional odds models with OTU abundances scaled by the total bacterial load as the dependent variable and used tick life stage as independent variables. Details on models are provided in Additional file 2: Text S1. We compared these models with the likelihood ratio test, using an implementation of Nagelkerke’s pseudo-R2 from the R package rcompanion (version 2.3.7) (35).
Transmission mode dataset
A total of 1,130 I. ricinus ticks of all developmental stages were collected in 2019 from two locations (ST and AW) and tested individually with qPCR for the presence of S. ixodetis, R. helvetica, B. burgdorferi s.l., A. phagocytophilum, and N. mikurensis for which primers and probes have been developed and published before (36–39). In addition, ticks were tested for Rickettsiella spp., M. mitochondrii for which primers and probes were designed in this study. All targeted genes and qPCR protocol are provided in Additional file 2: Text S2. Subsequently, pathogens and symbionts were assigned the transmission mode based on their presence or absence in the larval stage indicating vertical or horizontal transmission, respectively.
Dataset on vertically-transmitted symbionts (Symbiont dataset)
To determine the geographic distribution and prevalence of tick symbionts, we analysed a total of 16,555 ticks, which were collected in a previous cross-sectional study (22, 40). Briefly, questing nymphs of I. ricinus were collected from 19 locations in forested areas in the Netherlands in 2013 and 2014 (Fig. 1, triangles and dots). Details on data collection were described previously (22, 40). We tested questing individual nymphs of I. ricinus for the presence of S. ixodetis, R. helvetica, Rickettsiella spp., and M. mitochondrii.
Dataset on horizontally-transmitted pathogens (Pathogen dataset)
In addition, we included in our analysis data on the prevalence and distribution of tick-borne pathogens (Table 1). This data was generated from the same tick collection, from 13,967 of 16,555 ticks, in the study of Takumi et al. (29). The pathogens included A. phagocytophilum, N. mikurensis, B. miyamotoi, and three genospecies of B. burgdorferi s.l. as follows: B. afzelii, B. garinii, and B. valaisiana. Data on B. garinii and B. valaisiana were combined for further analysis as they are both considered bird-borne pathogens.
Relative occurrence of vertically-transmitted tick symbionts and horizontally-transmitted pathogens
Based on the presence and absence of microorganisms detected with qPCR (Symbiont and Pathogen datasets), we assigned each tick a haplotype, an integer between 0 and 2n − 1. Each integer corresponds to one of 2n distinct outcomes for a series of n qPCR tests performed on the tick. Haplotype frequencies were arranged in a table where rows are the sampling sites and columns are the observed haplotypes. Each row was divided by the row sum. The column mean was subtracted from each column. Principal component analysis of the data table was performed by applying the singular value decomposition (41).