Exploring the bacteriome in anthropophilic ticks: To investigate the vectors for diagnosis

Objective The aim of this study was to characterize the bacterial microbiome of hard ticks with affinity to bite humans in La Rioja (North of Spain). Methods A total of 88 adult ticks (22 Rhipicephalus sanguineus sensu lato, 27 Haemaphysalis punctata, 30 Dermacentor marginatus and 9 Ixodes ricinus) and 120 I. ricinus nymphs (CRETAV collection, La Rioja, Spain), representing the main anthropophilic species in our environment, were subjected to a metagenomic analysis of the V3-V4 region of the 16S rRNA gene using an Illumina MiSeq platform. Data obtained with Greengenes database were refined with BLAST. Four groups of samples were defined, according to the four tick species. Results Proteobacteria was the predominant phylum observed in all groups. Gammaproteobacteria was the most abundant class, followed by Alphaproteobacteria for R. sanguineus, H. punctata and D. marginatus but the relative abundance of reads for these classes was reversed for I. ricinus. This tick species showed more than 46% reads corresponding to ‘not assigned’ OTUs (Greengenes), and >97% of them corresponded to ‘Candidatus Midichloriaceae’ using BLAST. Within Rickettsiales, ‘Candidatus Midichloria’, Rickettsia, Ehrlichia, ‘Candidatus Neoehrlichia’ and Wolbachia were detected. I. ricinus was the most alpha-diverse species. Regarding beta-diversity, I. ricinus and H. punctata samples grouped according to their tick species but microbial communities of some R. sanguineus and D. marginatus specimens clustered together. Conclusions The metagenomics approach seems useful to discover the spectrum of tick-related bacteria. More studies are needed to identify and differentiate bacterial species, and to improve the knowledge of tick-borne diseases in Spain.


Introduction / Objective
The identification of microorganisms from biological samples has been dominated by the use of traditional culture-dependent methods and conventional molecular biology techniques (mostly polymerase chain reaction, PCR). The isolation of most tick-borne bacteria in synthetic media or in cell culture is difficult to obtain, and a high number of microbes remain uncultured. For the last two decades, the identification of Rickettsia spp. and other tick-associated pathogens has been mainly based on the use of specific PCR assays and sequence analysis [1,2]. Until recently, most studies focused on the detection of pathogens in vectors were able to detect a unique or a few microorganisms in a single assay. Metagenomic approaches, based on the development of the Next Generation Sequencing (NGS) techniques, and primary focused on the 16S rRNA study combined with bioinformatics tools, is revolutionizing the research in the fields of epidemiology and diagnosis of infectious diseases, among others, overcoming the limitation of detecting only one or few microorganisms at a time [3]. Metagenomic analysis can reveal the complexity of the microbiota of a given sample [4]. The number of pathogens associated with ticks has increased over the last years. Currently, there is a worldwide rising incidence of patients with a history of a tick-bite [5,6]. The importance of tick-borne diseases (TBDs) as a growing threat for public health has been recently underlined, and 'what is not sought, is not found' [7]. As ticks are able to transmit different microorganisms at one bite, it is necessary to be aware of possible co-infections. To investigate the microbial community composition harbored by ticks can facilitate the knowledge about the interactions among tick-associated microorganisms, the discovery of new uncultured microorganisms and subsequently, their implications as human pathogens.
Up to date, reports about metagenomics to investigate bacterial diversity of tick species are scarce. Our aim was to characterize the bacterial microbiome of hard ticks with affinity to bite humans in La Rioja (North of Spain).

Tick samples
A total of 280 questing ticks (130 adults and 150 nymphs) belonging to the main species with affinity to bite humans in La Rioja (Ixodes ricinus, Rhipicephalus sanguineus sensu lato, Dermacentor marginatus and Haemaphysalis punctata) were selected from the -80˚C freezer of the CRETAV collection (CIBIR, La Rioja, Spain) for the study of their bacterial profile. Ixodes ricinus is the most common arthropod vector of human diseases, and particularly nymphs of this species are the most frequent stage attacking humans in La Rioja [8]. Therefore, I. ricinus nymphs were also included in the study, in addition to adult specimens.
Ticks had been obtained from the field in La Rioja by flagging methods or by direct capture, either in urban habitats or in natural areas where outdoor activities are usually practiced, with the subsequent risk of infestation for humans (S1 Table). Specimens had been classified using taxonomic keys [9,10] and kept frozen at -80˚C while still alive. Before DNA extraction, a half from every adult tick (longitudinally cut) was immediately frozen again at -80˚C.

DNA extraction
For DNA extraction, ticks were manipulated under sterile conditions in a Class II biosafety cabinet using cycles of UV light prior and between uses to prevent contamination. All the tools were also irradiated with UV light for at least 15 min. Sterile single-use instruments were used whenever possible. Non-disponsable material was sterilized between samples (e.g. forceps were rinsed in 70% ethanol and flamed). Ticks were surface-sterilized by immersion and shaking in 70% ethanol for two min. followed by rinsing twice in sterile deionized water (one min. each). All the solutions were sterile. Ticks were dried on autoclaved sterile filter paper, transferred to sterile petri dishes and cut into small fragments that were collected in sterile tubes. The DNA was extracted using DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany), following the manufacturer's instructions except for an overnight digestion and a final elution in 25 μL of warm (at 56˚C) elution buffer. All the kit reagents had been previously tested for the absence of microorganisms using a pan-bacterial PCR [11]. Moreover, negative controls of extraction corresponding to extraction tubes without tick samples were included in parallel. DNA was quantified with a Qubit 3.0 fluorometer (Thermo Scientific) using Qubit dsDNA HS (High Sensitivity) assay kit. The quality of DNA was assessed by capillary electrophoresis with Fragment Analyzer (AATI) using Genomic DNA 50kb kit. DNA of enough quantity and quality for the NGS study was obtained from 88 adult ticks (22 Rhipicephalus sanguineus s.l., 27 Haemaphysalis punctata, 30 Dermacentor marginatus and 9 Ixodes ricinus) and 120 I. ricinus nymphs (in pools of ten individuals each) (S1 Table).
DNA extraction, preparation of PCR master mix, and amplification were performed in separate rooms to prevent contamination.

16S rRNA gene amplification, library preparation and sequencing
A total of 12.5 ng DNA per sample were used for the amplification step. Primers targeting the hypervariable V3 and V4 regions of 16S rRNA gene were used [12]. Amplified regions were purified and indexed with Nextera XT Index kit (Illumina). The library quality was assessed on a Qubit 3.0 Fluorometer (Thermo Scientific) and Fragment Analyzer (AATI) using dsDNA reagent (35-5000bp) kit. Paired-end 300 bp sequences were obtained on an Illumina MiSeq platform.

Sequence processing and analysis
Quality controls of raw reads were carried out with FastQC software [13], and trimmed with the Trimmomatic software [14]. The V3-V4 amplified region (550-580 bp) was reconstructed through paired reads according the Quantitative Insights Into Microbial Ecology (QIIME) protocol (v1.9.1) [15]. Operational Taxonomic Units (OTUs) were defined as sequences with at least 97% similarity versus Greengenes database [16] using UClust clustering algorithm [17] and following the open-reference method described by QIIME [18]. OTUs with <0.01% relative abundance of the total read counts on a per-sample basis were removed (spurious and chimeric reads). Data were refined with BLAST tool against GenBank database using the consensus sequence from each OTU [19].
Four groups of samples were defined, according to the four tick species. Rarefaction curves were calculated prior to all analytical techniques in order to assess species richness from the samples. OTU abundance was normalized by Cumulative Sum Scaling (CSS) method with metagenomeSeq software [20] and barplots were constructed.

Statistical analysis
Alpha diversity and relative evenness of communities' analyses were calculated by Chao1, Fisher, Margalef, Observed OTUs, Phylogenetic diversity (PD) whole tree, Shannon, Simpson, and Singles indexes with QIIME. Similarity distance matrixes between species groups were calculated following Bray-Curtis, Weighted Unifrac and Unweighted Unifrac beta-diversity metrics. Principal Coordinate Analysis (PCoA) and Hierarchical Clustering Dendrograms (UPGMA) for each beta-diversity metric were drawn to visualize sample groupings. The Kruskall-Wallis (KW) test was calculated to study significant differences between species groups. Analysis was also performed with MicrobiomeAnalyst software [21].

Results
A total of 19,977,253 read counts (average counts per sample = 201,790) and 227 OTUs were observed. The rarefaction curves reached a plateau, demonstrating that bacterial diversity had been satisfactorily detected for all samples (S1-S4 Figs).
Sequences assigned to 'Candidatus Midichloriaceae' using BLAST (that corresponded to not assigned OTUs against Greengenes) appeared in all I. ricinus nymph pools, with relative abundance of reads that ranged from 2.01 to 52.02% depending on the sample (S5 Table).
Bacteria belonging to the order Borreliales were minority (0.52% abundance of reads). Specifically, Borrelia spp. belonging to B. burgdorferi sensu lato (B. garinii) and relapsing fever According to alpha-diversity measures, the mean alpha diversity was greater for I. ricinus, followed by D. marginatus, R. sanguineus and H. punctata. Differences in alpha diversity between I. ricinus and D. marginatus, between I. ricinus and H. punctata and between I. ricinus and R. sanguineus were statistically significant (p<0.01) using all but Singles index. The highest standard deviation of the mean appeared in R. sanguineus for all but Shannon, Simpson and Singles indexes, which showed the highest standard deviation of the mean in I. ricinus ( Table 2). Group differences using Chao1 are showed in Fig 3. Regarding beta diversity metrics (distance measure), using PCoA with Bray-Curtis or Weighted Unifrac distance index and Analysis of Group Similarities (ANOSIM) method at genus level, I. ricinus and H. punctata samples gathered according to their tick species. On the contrary, microbial communities of several specimens of R. sanguineus and D. marginatus groups clustered together, suggesting profile similarity (Fig 4).
At OTU level, the best correlation between samples and tick species was showed using Bray-Curtis index. All samples grouped according to their tick species, except four R. sanguineus specimens that clustered within H. punctata (n = 2) or I. ricinus (n = 2) (Fig 5).

Discussion
Many TBDs have been recognized for the first time in the last few years, and emerging tickborne pathogens are being detected [3, [22][23][24]. Not only the clinical observation but also the application of new diagnostic methods (based on culture and molecular biology assays) has contributed to this progress [3]. Nevertheless, TBDs are dangerously expanding and they constitute underestimated causes of human illness worldwide [5]. The implementation of NGS platforms aimed to diagnosis is being developed, although reports about the contribution of this technique to the clinical diagnostic of TBDs are sporadic [25]. Herein, the bacteriome of tick species with affinity to bite humans was analysed using the 16S metagenomic approach to investigate tick-related microorganisms and to improve the diagnosis of TBDs, particularly in cases with unknown etiologic agents. Data generated with NGS studies for tick microbiome characterization allow us to delve into microorganism interactions [26]. As reported by Estrada-Peña and Cabezas-Cruz (2019), recent findings about the tick microbiome are driving to a change of paradigm: 'most bacteria found in tick microbiome are fundamental for tick biological processes' [27]. We agree with that statement, although the 'traditional' point of view should not be forgotten. We have learned throughout history that microorganisms first detected in ticks, and for a long time considered non-pathogenic to humans, have been later implicated in human diseases (e.g. Rickettsia parkeri), even though some of them do not fulfil Koch's Postulates (e.g. 'Ca. N mikurensis', a notyet-cultivated bacterium) [3, [28][29][30][31]. The finding of an infectious agent in a vector could enable its involvement in human pathology, especially if repeatedly detected. Metagenomics can allow the identification of microorganisms carried by arthropod vectors in people with suspected TBDs, thus contributing to the etiological diagnostics. Clinicians should consider that infection with multiple TBDs is possible, especially in tick-endemic areas. In cases of co-infection with more than one pathogen the clinical symptoms may be longer and more severe than expected, and the diagnosis can be even more difficult [32,33]. Data analysis obtained from NGS methods can be promising for the simultaneous detection of tick-borne pathogens in patients suffering TBDs of unknown aetiology.
Bacteria corresponding to genus Wolbachia were also detected in our samples. Wolbachia spp. are obligate intracellular endosymbionts of arthropods and nematodes. There is evidence about the capacity of these bacteria to affect biology, physiology, immunity, ecology and evolution and reproduction of the hosts, and to influence other infectious diseases due to viruses, protozoa and filariae [42]. The co-occurrence of these microorganisms considered endosymbionts can constitute a valuable research field of future studies because the viability of ticks, or even of the pathogens that ticks are able to transmit, may depend on these endosymbionts.
According to our data, other examples of OTUs that could be better identified using BLAST were those that matched with Anaplasma, Borreliaceae and Entomoplasmatales. However, the identification was not possible for other 'not assigned' OTUs that showed 91% identity (the highest) with known sequences of Rickettsia spp. or Coxiella endosymbionts. These findings can be useful for a future targeted search of unknown bacteria associated with ticks, and their potential implications for human health.
According to our results, the composition of the microbiota of ticks was affected by sex and geography, as previously reported [45]. For instance, on the one hand, H. punctata males from our study showed higher relative abundance of reads for Rickettsiales than females of the species or other tick species. This pattern could be explained by different host preferences between males and females and/or influence of host hormones and/or higher adaptive capacity of the microorganism to the tick and/or relationships between tick microorganisms, among other factors. Nevertheless, our data refer to the abundance of reads but not to prevalence, and a bias may have occurred since females generally have much more of the endosymbiont than males. On the other hand, D. marginatus and R. sanguineus showed overlapping PCoA plots, maybe because specimens of both species were collected from the same site (Villalba de Rioja). There is preliminary evidence that ticks that are geographically close share microbes [45].
Of particular interest is our observation of the highest I. ricinus alpha diversity over the other tick species analyzed herein. The generalist behavior in host choice of I. ricinus could have played a major role in the great variability of this tick-associated microbiota. Nearly all the life cycle of this tick species is spent in the surface layers of soil or forest litter where environmental conditions influence its development. I. ricinus is the primary vector of a wide variety of pathogens with considerable impact on human and animal health [46]. Contrary to experiments that have demonstrated higher mortality rates of R. sanguineus infected with Rickettsia conorii than non-infected when exposed at low or high temperature [47], I. ricinus is a tick with potential to adapt to new climates as they change [48].
Unfortunately, the comparison of data between studies that evaluate tick microbiomes is complex since every research team is focused on different research interests. Variations in techniques, target regions of the 16S rRNA gene, reference taxonomic databases or source of tick samples may hinder comparisons. As an example, relevant information about the ecology of tick-associated microorganisms in ticks and in voles from a French area has been recently published [39]. However, our reads from Coxiella endosymbionts could not be accurately compared to those obtained by Estrada-Peña et al. (2018) due to differences in length of reads (V3-V4 vs. V4 region). A review of NGS strategies for the study of the microbiome of ticks shows an updated view of the current scene [49]. As the authors conclude, further studies aimed to assess the influences of the environments, the hosts or the ticks themselves on the diversity of the tick microbiomes are required. According to the authors, bacteriome tick findings must be completed with new ones focused on viruses and eukarya in ticks [49]. Herein, a picture of bacteriome of ticks in a certain environment is showed, although ticks also harbour viruses, protozoa, fungi, helminths, etc. [50] and plenty of questions remain unresolved. The technique has difficulties and possible bias due to: storage of samples, DNA extraction method, reagents contamination, amplified 16S rRNA regions, updating and maintenance of curated sequences by reference databases or multiple repeated partial sequences of GenBank database, among others. However, the metagenomic approach seems useful to discover the spectrum of bacteria carried by ticks. More studies are needed to identify and differentiate bacterial species, and to improve the knowledge of TBDs in Spain.