Skip to main content

A cross-sectional screening by next-generation sequencing reveals Rickettsia, Coxiella, Francisella, Borrelia, Babesia, Theileria and Hemolivia species in ticks from Anatolia



Ticks participate as arthropod vectors in the transmission of pathogenic microorganisms to humans. Several tick-borne infections have reemerged, along with newly described agents of unexplored pathogenicity. In an attempt to expand current information on tick-associated bacteria and protozoans, we performed a cross-sectional screening of ticks, using next-generation sequencing. Ticks seeking hosts and infesting domestic animals were collected in four provinces across the Aegean, Mediterranean and Central Anatolia regions of Turkey and analyzed by commonly used procedures and platforms.


Two hundred and eighty ticks comprising 10 species were evaluated in 40 pools. Contigs from tick-associated microorganisms were detected in 22 (55%) questing and 4 feeding (10%) tick pools, with multiple microorganisms identified in 12 pools. Rickettsia 16S ribosomal RNA gene, gltA, sca1 and ompA sequences were present in 7 pools (17.5%), comprising feeding Haemaphysalis parva and questing/hunting Rhipicephalus bursa, Rhipicephalus sanguineus (sensu lato) and Hyalomma marginatum specimens. A near-complete genome and conjugative plasmid of a Rickettsia hoogstraalii strain could be characterized in questing Ha. parva. Coxiella-like endosymbionts were identified in pools of questing (12/40) as well as feeding (4/40) ticks of the genera Rhipicephalus, Haemaphysalis and Hyalomma. Francisella-like endosymbionts were also detected in 22.5% (9/40) of the pools that comprise hunting Hyalomma ticks in 8 pools. Coxiella-like and Francisella-like endosymbionts formed phylogenetically distinct clusters associated with their tick hosts. Borrelia turcica was characterized in 5% (2/40) of the pools, comprising hunting Hyalomma aegyptium ticks. Co-infection of Coxiella-like endosymbiont and Babesia was noted in a questing R. sanguineus (s.l.) specimen. Furthermore, protozoan 18S rRNA gene sequences were detected in 4 pools of questing/hunting ticks (10%) and identified as Babesia ovis, Hemolivia mauritanica, Babesia and Theileria spp.


Our metagenomic approach enabled identification of diverse pathogenic and non-pathogenic microorganisms in questing and feeding ticks in Anatolia.


Ticks (class Arachnida, subclass Acari) are the most significant arthropod vectors, along with mosquitoes, participating in the transmission of pathogens to humans [1]. A diverse group of infectious agents including viruses, bacteria and protozoans can be transmitted by ticks, surpassing most arthropods in terms of vector potential [2]. Tick-borne infections of humans are of zoonotic origin, with pathogens maintained in natural cycles involving tick vectors and animal hosts [3]. Frequently, humans are accidental, dead-end hosts that do not significantly contribute to the pathogen’s life-cycle. Various tick species occupy distinct ecological niches that define their distribution patterns and risk areas for tick-borne infections [4]. The past decades have witnessed the emergence and resurgence of several tick-borne infections with considerable impact on human and animal welfare [1, 5]. A deeper understanding of the epidemiology and potential public health threats of tick-borne infections rely on effective surveillance programmes to identify circulating pathogens in vectors and reliable diagnosis of vertebrate infections.

In addition to the tick-borne pathogens, a diverse group of commensal and symbiotic bacteria are described in ticks, usually co-circulating with the infectious agents [6]. Their biology and effect on tick life-cycle remain largely unexplored, despite evidence suggesting their involvement in fitness, nutritional adaptation, defense and immunity [7]. These microorganisms are also likely to interact with the replication and transmission of tick-borne pathogens, with potential implications for human and animal health [7, 8].

Turkey is located in Asia Minor and maintains a natural transmission zone for vector-borne infections between Asia, Africa and Europe [9]. The geographical regions of Anatolia, with diverse climate conditions, vegetation patterns, domestic animals and wildlife provide suitable habitats for perpetuating several arthropod vectors of disease, including ticks [9]. Several species of the families Ixodidae and Argasidae are present in the tick fauna of Turkey [10]. Human tick-borne infections have also been documented, caused by protozoans, nematodes, bacteria and viruses [9, 11]. We have recently reported the presence of several RNA viruses in ticks collected from various regions of Anatolia [12]. In the present study, we aimed to perform a cross-sectional screening by using next-generation sequencing (NGS) to characterize tick-associated bacteria and protozoans.


Specimen collection and processing

Ticks collected in several locations from Ankara and Cankiri provinces (central Anatolia), Mugla Province (western Anatolia, Aegean region) and Mersin Province (southern Anatolia, Mediterranean region) from April to October from 2014 to 2016 were evaluated. Questing ticks were captured on site by flagging as well as from infested domesticated animals: dogs (Canis familiaris); cattle (Bos taurus); and goats (Capra aegagrus hircus). The ticks were kept alive individually in vials, transferred to the laboratory and identified morphologically to the species level using several taxonomic keys [13,14,15,16,17]. Following identification, the specimens were pooled according to species and collection site up to a maximum of 22 individuals per pool and stored at -80 °C for further analysis.

Individual and pooled ticks with up to five specimens were homogenized using the SpeedMill PLUS (Analytik Jena, Jena, Germany), and total nucleic acid purification was performed by using BlackPREP tick DNA/RNA kit (Analytik Jena) according to the manufacturer’s instructions. Pools with six or more specimens were kept in 500–700 μl of Eagle’s minimal essential medium, supplemented with 1% L-glutamine and 5% fetal bovine serum. These pools were homogenized by vortexing with tungsten carbide beads (Qiagen, Hilden, Germany) and clarified by centrifugation for 4 min at 4000× rpm. Subsequently, the ground pools were aliquoted and subjected to nucleic acid extraction using High Pure Viral Nucleic Acid Kit (Roche Diagnostics, Mannheim, Germany)

Next-generation sequencing (NGS) and phylogenetic analysis

Purified nucleic acids from tick pools were reverse transcribed with random hexamer primers to double-stranded cDNA using SuperScript IV Reverse Transcriptase (Thermo Fisher Scientific, Hennigsdorf, Germany) and NEBNext mRNA Second Strand Synthesis Module (New England Biolabs, Frankfurt am Main, Germany). Agencourt AMPure XP Reagent (Beckman Coulter Biosciences, Krefeld, Germany) and Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany) were employed for cleanup, yield and size distribution determination. Fragmentation, adaptor ligation and amplification were carried out using NexteraXT DNA Library Preparation Kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s protocols. Sequencing runs were performed on an Illumina HiSeq (Illumina Inc.) instrument in paired-end mode.

The raw sequencing data was de-multiplexed and extracted in fastq format. Trimmomatic software was employed for trimming for quality and length with a phred score of 33 and a minimum length of 30 base pairs (bp) and removal of Illumina adaptors [18]. Obtained reads were aligned to the GenBank RefSeq databases of the National Center for Biotechnology Information (NCBI) for bacteria (v.17.03.2017), 16S ribosomal RNA (rRNA) (RefSeq rRNA, v.01.08.2017) and selected protozoa (in-house curated database, sequences available upon request, v.29.09.2017) using MALT (MEGAN alignment tool, v0.3.8) and MEGAN (Metagenome Analyzer, v. 6.12.3) [19, 20]. Aligned reads were extracted and assembled into contigs using Velvet (v.1.2.10) with a k-mer length of 31 [21]. The contigs were checked for heterogeneity by visual inspection and via pairwise identity values using Geneious software v.11.1.5 (Biomatters Ltd, Auckland, New Zealand). The 16S rRNA gene sequences were scanned for chimeras with divergence of > 3% from the closest parent using UCHIME2, implemented at the NCBI database [22]. For the near-complete genome and plasmid sequences, contigs and remaining reads were mapped to closely related strains. BLASTn, BLASTn optimized for highly similar sequences (MEGABLAST) and BLASTp algorithms were used for nucleotide and deduced amino acid similarity searches in the public databases implemented in the NCBI website ( [23]. Nucleotide and putative amino acid alignments and pairwise sequence comparisons were generated by using the CLUSTAL W program implemented within Geneious software [24]. Conserved protein domain and motif searches were performed using the web search tool ( and MOTIF Search ( in the PFAM database [25, 26]. The models for the phylogenetic and molecular evolutionary analyses were selected using the best-fit DNA/protein-substitution model tools of the MEGA v.6.06 software [27]. Phylogenetic trees were constructed using the maximum-likelihood method with the Tamura-Nei substitution model. The reliability of the inferred trees was evaluated by bootstrap analysis of 1000 replicates.


Two-hundred eighty ticks, comprising 179 female (63.9%), 100 male (35.7%) and 1 nymph (0.4%) specimens were evaluated in 40 pools, prepared according to species and collection site (Additional file 1: Table S1). A total of 10 tick species were identified among which Rhipicephalus bursa (n = 76; 27.1%), Hyalomma aegyptium (n = 49; 17.5%) and Haemaphysalis parva (n = 46; 16.4%) represented the most abundant species. A total of 8 pools with Ha. parva (n = 3), Rhipicephalus sanguineus (s.l.) (n = 3), Dermacentor marginatus (n = 1) and Rhipicephalus annulatus (n = 1) specimens were collected from animal hosts, whereas the remaining pools (n = 32, 80%) comprised questing/hunting ticks (Additional file 1: Table S1).

NGS provided trimmed read numbers of 67,753–53,026,910 (mean = 5,061,800; median = 2,430,793) in the tick pools (Additional file 1: Table S1). Tick-associated microbial sequences were detected in 26/40 pools (65%), with multiple microorganisms identified in 12 (46.2%) of these positive pools. Coxiella, Francisella, Rickettsia, Babesia, Borrelia, Theileria and Hemolivia sequences were characterized in reactive tick pools. Coxiella spp. were the most frequently detected microorganism, identified in 16 of the 40 pools (40%), followed by Francisella spp. (22.5%), Rickettsia spp. (17.5%) and other microorganisms (Table 1, Additional file 1: Table S2).

Table 1 Tick pools with detectable microorganism sequences

Rickettsia findings

Rickettsia spp. sequences were identified in a total of 7 tick pools comprising Ha. parva (n = 3), Hyalomma marginatum (n = 2), R. bursa (n = 1) and R. sanguineus (s.l.) (n = 1) specimens (Table 2). Rickettsia spp. were detected in 4 pools (57.1%) of questing/hunting and 3 pools (42.9%) of feeding ticks.

Table 2 Tick pools with Rickettsia spp. contigs. Size and GenBank accession numbers are provided

A near-complete Rickettsia genome was assembled from the pool P4 that comprised 13 feeding Ha. parva ticks. A total of 82,002 reads from this pool were aligned to the genomes of two Rickettsia strains, R. hoogstraalii strain Croatica and Rickettsia felis strain URRWXCal2 (CP000053). These were further assembled into 1516 contigs with an N50 length of 1012 bp and a total length of 1,176,263 bp. Pairwise comparison of this sequence revealed 98.3 and 89.6% identity with R. hoogstraalii and R. felis, respectively. The sequence was disrupted by several gaps of varying length and further sequencing to complete the genome was not feasible. Therefore, we extracted intact contigs for comparison, including complete 16S and 23S rRNA genes, citrate synthase (gltA), surface cell antigen 1 (sca1) and outer membrane protein A (ompA). These sequences were submitted to the GenBank database (Table 2) and the assembled genome sequence is available in FASTA format as Additional file 2.

A 2376 bp section of the Rickettsia putative conjugative plasmid was also detected in pool P4. Pairwise comparison showed 80.8 and 97.8% identity with R. australis and R. hoogstraalii plasmids, respectively. A comparative alignment is provided in Additional file 3. Motifs of TraA_Ti conjugative transfer protein, MobA/MobL family mobilization protein and TraA conjugal transfer relaxase were identified within the sequence.

In addition to the complete 16S rRNA gene sequence in pool 4, 16S rRNA gene contigs of 1232–1433 bp were obtained in feeding Ha. parva pools (P5 and P6), hunting H. marginatum pools (P35 and P37), a questing R. bursa pool (P18) and a questing R. sanguineus (s.l.) pool (P21) (Table 2). Contigs in pools P5 and P6 revealed 97–98% identity with R. hoogstraalii in BLASTn and MEGABLAST searches. In the maximum-likelihood tree, the P4 and P6 contigs grouped with R. hoogstraalii, with separate clustering of P18, P5-P21 and P35–P37 (Fig. 1).

Fig. 1
figure 1

Maximum-likelihood analysis of the Rickettsia partial 16S rRNA gene sequences (1294 nucleotides). The tree was constructed using the Tamura-Nei model, with a bootstrap analysis of 1000 replicates. Sequences characterized in this study are given in bold and indicated with a symbol, GenBank accession number, pool code and host tick species. Rickettsia strains are indicated by GenBank accession number, microorganism and strain/isolate name. Bootstrap values lower than 60 are not shown. Coxiella burnetii was included as the outgroup

The complete gltA-coding region extracted from the Rickettsia genome in pool P4 demonstrated 99 and 100% identities to R. hoogstraalii prototype isolate, in nucleotide and deduced amino acid comparisons, respectively. The gltA-coding sequences were also obtained in pools P5 and P6. These constituted 355 and 821 bp stretches which were identical to the P4 sequence. The sca1 and ompA contigs in pool P4 were also highly similar to R. hoogstraalii, with 98.4–99.8% and 99.4% nucleotide and amino acid identity, respectively. The 428-nucleotide sca1 contig from the pool P6 was also identical to the sequence in pool P4. In addition, a longer section of ompA could be obtained from P6, which showed 1.4% divergence from P4 sequence and 98.6 and 99.2% nucleotide and amino acid identity, respectively, to R. hoogstraalii.

Overall, the obtained sequences enabled identification of the Rickettsia strain in pools P4 and P6 (Table 2), and the analysis of the 16S region could not provide data sufficient for strain discrimination in pools P5, P18, P21, P37 and P37 (Fig. 1). The available 355-nucleotide gltA contig from pool P5 revealed similar identity to several Rickettsia in BLASTn and MEGABLAST searches, therefore the precise identification of the strain in this pool also remained obscure.

Coxiella, Francisella and Borrelia findings

Bacterial 16S rRNA gene sequences other than Rickettsia were characterized in 26 pools (65%), 22 (84.6%) of which included questing/hunting ticks. In 16 tick pools (40%), comprising Rhipicephalus, Hyalomma, Haemaphysalis and Dermacentor specimens, 16S rRNA gene sequences with varying similarities to several Coxiella-like endosymbionts (CLE) were detected. CLEs were present in pools of questing (n = 12) as well as feeding (n = 4) ticks (Table 1). The sequences comprised 1088–1180 bp with up to 4.5% diversity. In the maximum-likelihood tree, three distinct clusters were observed (Fig. 2). The sequences from Rhipicephalus and Hyalomma species (MH645186–96) grouped with endosymbionts of Rhipicephalus spp., while the sequences from feeding Ha. parva (MH645183–5) remained distinct, sharing a common ancestor with CLE from Ixodes spp. The sequence originating from the questing D. marginatus tick pool (MH645197) also formed another clade with endosymbionts of the same tick species (Fig. 2).

Fig. 2
figure 2

Maximum-likelihood analysis of the Coxiella partial 16S rRNA gene sequences (1086 nucleotides). The tree was constructed using the Tamura-Nei model, with a bootstrap analysis of 1000 replicates. Sequences characterized in this study are given in bold and indicated with a symbol, GenBank accession number, pool code and host tick species. Bacterial strains are indicated by GenBank accession number, microorganism and strain/isolate name. Bootstrap values lower than 60 are not shown. Legionella pneumophila strain Philadelphia 1 was included as the outgroup

The 16S rRNA gene sequences related to Francisella species were identified in 9 pools (22.5%) with hunting/questing specimens. In contrast to CLE, these sequences belonged to more abundant Hyalomma ticks (8/9) (Table 1). They comprised contigs of 1290–1516 bp and formed two groups, namely P11, P20, P22, P25, P26 and P35, P37, P38, P40, with less than 1% intragroup divergence and 98.9% identity between groups. These groups, distinct from pathogenic Francisella and Wolbachia endosymbionts, could also be distinguished phylogenetically, as the pools with H. marginatum ticks (MH645186, MH645198-MH645200) formed high bootstrap supported clades with Francisella-like endosymbionts (FLE) of Hyalomma rufipes (Fig. 3). The other group, detected in H. aegyptium and R. sanguineus (s.l.) (MH645201-MH645205), clustered with sequences from H. aegyptium and Amblyomma spp., whereas FLE from Ixodes, Dermacentor and Haemaphysalis ticks remained distinct.

Fig. 3
figure 3

Maximum-likelihood analysis of the Francisella and Wolbachia partial 16S rRNA gene sequences (375 nucleotides). The tree was constructed using the Tamura-Nei model, with a bootstrap analysis of 1000 replicates. Sequences characterized in this study are given in bold and indicated with a symbol, GenBank accession number, pool code and host tick species. Bacterial strains are indicated by GenBank accession number, microorganism and strain/isolate name. Bootstrap values lower than 60 are not shown. Pasteurella multocida subsp. gallicida strain NCTC 10204 was included as the outgroup

The last group of 16S rRNA gene contigs constituted two sequences of 1361 and 1364 bp (MH628249-MH628250), identical except for 1–3 nucleotide terminal overhangs. They were detected in pools of hunting H. aegyptium ticks and were identical to B. turcica and Borrelia sp. recovered from H. aegyptium in Turkey and Amblyomma geoemydae in Japan. They grouped together along with several tick-associated Borrelia species in the maximum-likelihood tree, forming a separate clade distinct from relapsing fever and Lyme disease Borrelia (Fig. 4).

Fig. 4
figure 4

Maximum-likelihood analysis of the Borrelia partial 16S rRNA gene sequences (1352 nucleotides). The tree was constructed using the Tamura-Nei model, with a bootstrap analysis of 1000 replicates. Sequences characterized in this study are given in bold and indicated with a symbol, GenBank accession number, pool code and host tick species. Bacterial strains are indicated by GenBank accession number, microorganism and strain/isolate name. Bootstrap values lower than 60 are not shown

Babesia, Theilera and Hemolivia findings

Eukaryotic 18S rRNA gene contigs were obtained in four questing tick pools (10%). In three pools comprising Rhipicephalus spp., sequences related to Babesia and Theileria were detected. BLASTn analysis of the longest sequence (1454 bp) in the pool P34 (MH618772) revealed highest similarity rates of 99% to B. ovis. It further grouped phylogenetically with B. ovis with high bootstrap values, confirming the identification (Fig. 5). The 544-bp sequence in pool P23 (MH618773) displayed 97–98% identity to several Babesia sp. detected in ticks, but no definitive strain identification could be established. The recently described, presumably novel Babesia sequences from ticks and goats from Turkey [28, 29] also revealed 96.1–97% identity and were distantly related to this sequence (Fig. 5). The 1339-bp sequence obtained from pool P11 (MH618774) showed 88–89% identity to various Babesia and Theileria spp. and clustered with the Theileria spp. in the maximum-likelihood analysis (Fig. 5).

Fig. 5
figure 5

Maximum-likelihood analysis of the Babesia and Theilera partial 18S rRNA gene sequences (794 nucleotides). The tree was constructed using the Tamura-Nei model, with a bootstrap analysis of 1000 replicates. Sequences characterized in this study are given in bold and indicated with a symbol, GenBank accession number, pool code and host tick species. Protozoan strains are indicated by GenBank accession number, microorganism and strain/isolate name. Bootstrap values lower than 60 are not shown

Finally, a 472 bp sequence, with 99–100% identity to several Hemolivia mauritanica isolates, was obtained from a pool of hunting H. aegyptium ticks (Table 1). Despite the availability of a relatively short segment, the sequence (MH618775) grouped with He. mauritanica isolates in the maximum-likelihood tree (Fig. 6).

Fig. 6
figure 6

Maximum-likelihood analysis of the Hemolivia and Hepatozoon partial 18S rRNA gene sequences (465 nucleotides). The tree was constructed using the Tamura-Nei model, with a bootstrap analysis of 1000 replicates. Sequences characterized in this study are given in bold and indicated with a symbol, GenBank accession number, pool code and host tick species. Protozoan strains are indicated by GenBank accession number, microorganism and strain/isolate name. Bootstrap values lower than 50 are not shown. Babesia sp. isolate Kashi1 was included as the outgroup


We performed a cross-sectional screening for tick-associated bacteria and protozoans, using an NGS-based strategy in pools of field-collected ticks from various regions of Turkey. We adopted a straightforward approach for NGS, using standard and widely-used commercial assays for nucleic acid purification, cDNA and sequencing library preparation, performed without major modifications. We could detect Rickettsia in 17.5% of the pools, including questing and feeding ticks. The obtained sequences comprised a near-complete genome, partial conjugative plasmid as well as 16S rRNA, ompA, sca1 and gltA gene segments (Table 2). The strain could be identified as R. hoogstraalii in two pools comprising feeding Ha. parva ticks, while the available data were insufficient for precise strain characterization in the remaining specimens. Rickettsia (order Rickettsiales genus Rickettsia) are intracellular Gram-negative bacteria that infect eukaryotic cells [30]. Several species are recognized, currently organized within distinct groups, according to the genome-wide sequence data [31]. Rickettsia hoogstraalii is closely related to R. felis and both strains are classified within the spotted fever group that includes species causing tick-borne infections in humans [32]. Despite in vitro cytopathic effects on various cell lines, the pathogenesis of R. hoogstraalii in vertebrate hosts remains unknown [32]. The isolation was accomplished from Haemaphysalis sulcata in Croatia and it has been detected in several tick species from various countries, including Cyprus, Ethiopia, Japan, Spain, the Indian Ocean islands and the USA [32,33,34]. Rickettsia hoogstraalii was initially identified in Turkey in 2014 in Ha. parva and the follow-up efforts have detected this strain in Ha. parva and Haemaphysalis punctata ticks in Central Anatolia [11, 35,36,37]. In addition to the near-complete genome with 98.3% identity to the prototype genome, we characterized a segment of the rickettsial plasmid and identified protein motifs with conjugative transfer functions. Despite their strictly intracellular life-cycle and reductive genomic evolution, several Rickettsia spp. have been shown to possess plasmids, with the possibility of horizontal, plasmid-mediated DNA exchange in ticks [38, 39].

We could not characterize the detected Rickettsia in five tick pools, despite the availability of relatively long 16S rRNA gene sequences (Table 2). The 16S gene is highly conserved among Rickettsia, where the similarity level between two species exceeds 97.2% [31]. This constitutes an impediment for significant inferences of intragenus phylogeny and hampers strain identification, which can be overcome by sequencing citrate synthase or outer surface proteins and surface cell antigens [31, 40]. Such data for the tick pools in question could not be produced in this setting, due to the relatively limited number of target sequence reads obtained.

Ticks have been documented to harbor diverse bacterial strains engaged in facultative or obligate endosymbiotic interactions with their hosts [6, 7]. Many distinct genera of bacteria, including strains collectively named as CLE, FLE and Rickettsia-like endosymbionts, have been identified in ticks [41]. CLE and FLE were detected mostly in ticks, with varying infection rates in different species. CLE are ubiquitous, geographically widespread and detected in several tick species as well as in the spleen of wild mammals [7, 42, 43]. We characterized 16S rRNA gene sequences of CLE in 40% of the screened tick pools, which were the most frequently detected bacteria in the study cohort (Table 2). We detected CLEs in feeding as well as hunting/questing ticks, and observed a differential phylogenetic clustering of sequences according to the tick species (Fig. 2). The genus Coxiella is genetically divergent, with at least four highly divergent clades recognized, and CLE hosted by ticks are present in all clades [41, 44]. Interestingly, the phylogenetic patterns indicate that the well-known human pathogen Coxiella burnetii, the etiological agent of Q fever, has evolved from a tick-associated Coxiella [44]. This was also observed in our analysis where C. burnetii shared a common ancestor with Ixodes- and Haemaphysalis-associated sequences and formed a distinct clade among Coxiella (Fig. 2).

We further detected FLE in our cohort, with an incidence of 22.5%, occurring in hunting/questing ticks. FLEs are considered as an obligate symbiont alternate to CLE in some tick species and are, like CLE, genetically related to their pathogenic counterpart: Francisella tularensis, the etiological agent of tularemia [7, 41]. FLEs are widely distributed in Europe and identified in various tick species [45]. Interestingly, we observed a preferential detection of FLE in Hyalomma ticks (Table 1). Moreover, the FLE sequences formed phylogenetically distinct clusters associated with their tick hosts, suggesting differential evolutionary patterns in various hosts and ecological niches (Fig. 3). All FLE-related sequences remained distinct from pathogenic Francisella.

NGS provided Borrelia 16S rRNA gene sequences in 5% of the pools comprising hunting H. aegyptium ticks (Table 1). These sequences were identical to the previously characterized Borrelia turcica isolated from the same tick species [46, 47]. Borrelia turcica and closely related bacteria (Borrelia sp. tAG) are divergent from species associated with Lyme disease and relapsing fever, forming a third phylogenetic lineage within the genus Borrelia [48, 49], as observed in our analysis (Fig. 4). Also called reptile-associated Borrelia, members of this lineage are widely distributed, infecting various tick species [48,49,50,51]. Detected only in ticks or blood collected from tortoises so far, the consequences of human or animal exposure by these Borrelia are currently unknown [48, 49]. However, given the detection of several zoonotic agents and sporadic feeding on humans of H. aegyptium ticks, human infection by Borrelia turcica seems possible [51]. Therefore, characterization of the infecting strain in symptomatic individuals may provide information on the pathogenic potential of the members of this Borrelia lineage.

The outcomes of vertebrate infections with bacterial endosymbionts or apparently non-pathogenic bacteria in ticks remain obscure. No information regarding FLEs and Borrelia turcica pathogenicity is currently available. However, mild human infections caused by Coxiella-like bacteria were documented [52] as well as asymptomatic equine and severe avian pet infections [53,54,55]. These findings suggest that occasional human infections may occur. Therefore, they should be investigated in tick-associated infections in humans and animals without detectable pathogens. Another aspect is that non-pathogenic bacteria may interfere in the replication of tick-borne pathogens, influencing their abundance in vectors and transmission to vertebrate hosts. We could identify CLE and Babesia co-infection in a single questing R. sanguineus complex specimen (P23; Table 1, Additional file 1: Table S1), which indicates that co-infections are not extremely rare and can be detected by using appropriate methods in field-collected ticks.

The NGS-based approach further provided protozoan 18S rRNA gene sequences in 10% of the tick pools where Babesia, Theilera and Hemolivia spp. were identified (Table 1). The microorganisms could be characterized as B. ovis and He. mauritanica in R. bursa and H. aegyptium pools, respectively by pairwise comparisons and inferred phylogenies (Figs. 5 and 6). Babesiosis is prevalent in Turkey and B. ovis, the etiological agent of sheep babesiosis, was previously identified in R. bursa and Rhipicephalus turanicus ticks [35, 56, 57]. In addition to B. ovis, several other species were reported as well as a proposed novel Babesia in ticks and goats [9, 11, 58, 59]. Despite reliable identification of B. ovis, the Babesia sequence in the R. sanguineus complex pool (P23) remained unidentified due to insufficient sequence data.

The major shortcoming of this study is the relatively low number of total and target sequences obtained from tick pools. NGS-based approaches, when optimized for DNA/RNA deep sequencing, can produce up to 109 reads [60], which surpasses the overall sequencing efficiency observed in this study. The lack of tick-associated microorganisms in 30% of the pools can be attributed to this particular limitation, along with the overabundance of background signals from the host. Several steps and factors within the NGS workflow may affect sequencing efficiency and depth, producing potential biases in the representation of the original sequences [61]. Our strategy involved utilization of standardized specimen processing and library preparation, comparable to PCR-based pathogen screening. Targeted amplification and NGS of the bacterial and protozoan rRNA or pathogenic microorganisms can be alternate strategies, such as those we previously developed for viral hemorrhagic fever agents [62]. For a detailed investigation of the tick microbiome, individual ticks should be evaluated with a deeper sequencing strategy which we plan to employ for co-infected specimens in upcoming studies.


Using an NGS-based approach, we detected bacteria of the genera Rickettsia, Coxiella, Francisella, Borrelia and protozoans of the genera Babesia, Theileria and Hemolivia in questing and feeding ticks. A near-complete genome and the conjugative plasmid of R. hoogstraalii were assembled, along with several coding and non-coding Rickettsia genes in tick pools. Moreover, CLE and FLEs with the hosting tick species were documented.



Coxiella-like endosymbionts


Francisella-like endosymbionts


Next-generation sequencing


ribosomal RNA


  1. Dantas-Torres F, Chomel BB, Otranto D. Ticks and tick-borne diseases: a One Health perspective. Trends Parasitol. 2012;28:437–46.

    Article  Google Scholar 

  2. de la Fuente J, Estrada-Pena A, Venzal JM, Kocan KM, Sonenshine DE. Overview: ticks as vectors of pathogens that cause disease in humans and animals. Front Biosci. 2008;13:6938–46.

  3. Baneth G. Tick-borne infections of animals and humans: a common ground. Int J Parasitol. 2014;44:591–6.

    Article  Google Scholar 

  4. Dantas-Torres F. Climate change, biodiversity, ticks and tick-borne diseases: the butterfly effect. Int J Parasitol Parasites Wildl. 2015;4:452–61.

    Article  Google Scholar 

  5. Kilpatrick AM, Randolph SE. Drivers, dynamics, and control of emerging vector-borne zoonotic diseases. Lancet. 2012;380:1946–55.

    Article  Google Scholar 

  6. Ahantarig A, Trinachartvanit W, Baimai V, Grubhoffer L. Hard ticks and their bacterial endosymbionts (or would be pathogens). Folia Microbiol. 2013;58:419–28.

    Article  CAS  Google Scholar 

  7. Bonnet SI, Binetruy F, Hernandez-Jarguin AM, Duron O. The tick microbiome: why non-pathogenic microorganisms matter in tick biology and pathogen transmission. Front Cell Infect Microbiol. 2017;7:236.

    Article  Google Scholar 

  8. Greay TL, Gofton AW, Paparini A, Ryan UM, Oskam CL, Irwin PJ. Recent insights into the tick microbiome gained through next-generation sequencing. Parasit Vectors. 2018;11:12.

    Article  Google Scholar 

  9. Inci A, Yildirim A, Duzlu O, Doganay M, Aksoy S. Tick-borne diseases in Turkey: a review based on One Health perspective. PLoS Negl Trop Dis. 2016;10:e0005021.

    Article  Google Scholar 

  10. Bursali A, Keskin A, Tekin S. A review of the ticks (Acari: Ixodida) of Turkey: species diversity, hosts and geographical distribution. Exp Appl Acarol. 2012;57:91–104.

    Article  Google Scholar 

  11. Karasartova D, Gureser AS, Gokce T, Celebi B, Yapar D, Keskin A, et al. Bacterial and protozoal pathogens found in ticks collected from humans in Corum Province of Turkey. PLoS Negl Trop Dis. 2018;12:e0006395.

    Article  Google Scholar 

  12. Brinkmann A, Dinçer E, Polat C, Hekimoğlu O, Hacıoğlu S, Földes K, et al. A metagenomic survey identifies Tamdy orthonairovirus as well as divergent phlebo-, rhabdo-, chu- and flavi-like viruses in Anatolia, Turkey. Ticks Tick Borne Dis. 2018;9:1173–83.

    Article  Google Scholar 

  13. Apanaskevich DA, Horak IG. The genus Hyalomma Koch, 1844: V. Re evaluation of the taxonomic rank of taxa comprising the H. (Euhyalomma) marginatum Koch complex of species (Acari: Ixodidae) with redescription of all parasitic stages and notes on biology. Int J Acarol. 2008;34:13–42.

    Article  Google Scholar 

  14. Walker AR, Bouattour A, Camicas JL, Estrada-Peña A, Horak IG, Latif AA, et al. Ticks of Domestic Animals in Africa: A Guide to Identification of Species. 1st ed. Edinburgh: Bioscience Reports; 2003.

    Google Scholar 

  15. Filippova NA. Fauna of Russia and neighbouring countries. Ixodid ticks of subfamily Amblyomminae. St. Petersburg, Russia: Nauka Publishing House; 1997.

    Google Scholar 

  16. Estrada-Peña A, Bouattour A, Camicas JL, Walker AR. Ticks of Domestic Animals in the Mediterranean Region. 1st ed. Zaragoza: University of Zaragoza Press; 2004.

    Google Scholar 

  17. Walker JB, Keirans JE, Horak IG. The Genus Rhipicephalus (Acari, Ixodidae): A Guide to the Brown Ticks of the World. Revised ed. Cambridge: Cambridge University Press; 2000.

    Book  Google Scholar 

  18. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  Google Scholar 

  19. Huson DH, Beier S, Flade I, Gorska A, El-Hadidi M, Mitra S, et al. MEGAN Community Edition - interactive exploration and analysis of large-scale microbiome sequencing data. PLoS Comput Biol. 2016;12:e1004957.

    Article  Google Scholar 

  20. Herbig A, Maixner F, Bos KI, Zink A, Krause J, Huson DH. MALT: Fast alignment and analysis of metagenomic DNA sequence data applied to the Tyrolean Iceman. bioRxiv. 2016.

  21. Zerbino DR, Birney E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18:821–9.

  22. Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics. 2011;27:2194–200.

    Article  CAS  Google Scholar 

  23. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    Article  CAS  Google Scholar 

  24. Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994;22:4673–80.

    Article  CAS  Google Scholar 

  25. Marchler-Bauer A, Derbyshire MK, Gonzales NR, Lu S, Chitsaz F, Geer LY, et al. CDD: NCBI’s conserved domain database. Nucleic Acids Res. 2015;43:D222–6.

  26. Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Eddy SR, et al. The Pfam protein families database. Nucleic Acids Res. 2002;30:276–80.

    Article  CAS  Google Scholar 

  27. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.

    Article  CAS  Google Scholar 

  28. Ozubek S, Aktas M. Molecular evidence of a new Babesia sp. in goats. Vet Parasitol. 2017;233:1–8.

    Article  CAS  Google Scholar 

  29. Ozubek S, Aktas M. Molecular evidence for a novel species of Babesia in unfed Rhipicephalus sanguineus sensu lato (Acari: Ixodidae). J Med Entomol. 2018;55:1271–6.

    Article  Google Scholar 

  30. Kang YJ, Diao XN, Zhao GY, Chen MH, Xiong Y, Shi M, et al. Extensive diversity of Rickettsiales bacteria in two species of ticks from China and the evolution of the Rickettsiales. BMC Evol Biol. 2014;14:167.

    Article  Google Scholar 

  31. Merhej V, Angelakis E, Socolovschi C, Raoult D. Genotyping, evolution and epidemiological findings of Rickettsia species. Infect Genet Evol. 2014;25:122–37.

    Article  Google Scholar 

  32. Sentausa E, El Karkouri K, Nguyen TT, Caputo A, Raoult D, Fournier PE. Genome sequence of Rickettsia hoogstraalii, a geographically widely distributed tick-associated bacterium. Genome Announc. 2014;2:e01171–14.

    PubMed  PubMed Central  Google Scholar 

  33. Duh D, Punda-Polic V, Avsic-Zupanc T, Bouyer D, Walker DH, Popov VL, et al. Rickettsia hoogstraalii sp. nov., isolated from hard- and soft-bodied ticks. Int J Syst Evol Microbiol. 2010;60:977–84.

    Article  Google Scholar 

  34. Duh D, Punda-Polic V, Trilar T, Petrovec M, Bradaric N, Avsic-Zupanc T. Molecular identification of Rickettsia felis-like bacteria in Haemaphysalis sulcata ticks collected from domestic animals in southern Croatia. Ann N Y Acad Sci. 2006;1078:347–51.

    Article  CAS  Google Scholar 

  35. Orkun Ö, Karaer Z, Cakmak A, Nalbantoglu S. Spotted fever group rickettsiae in ticks in Turkey. Ticks Tick Borne Dis. 2014;5:213–8.

    Article  Google Scholar 

  36. Orkun Ö, Karaer Z, Cakmak A, Nalbantoglu S. Identification of tick-borne pathogens in ticks feeding on humans in Turkey. PLoS Negl Trop Dis. 2014;8:e3067.

    Article  Google Scholar 

  37. Keskin A, Bursali A, Keskin A, Tekin S. Molecular detection of spotted fever group Rickettsiae in ticks removed from humans in Turkey. Ticks Tick Borne Dis. 2016;7:951–3.

    Article  Google Scholar 

  38. El Karkouri K, Pontarotti P, Raoult D, Fournier PE. Origin and evolution of rickettsial plasmids. PLoS One. 2016;11:e0147492.

    Article  Google Scholar 

  39. Weinert LA, Welch JJ, Jiggins FM. Conjugation genes are common throughout the genus Rickettsia and are transmitted horizontally. Proc Biol Sci. 2009;276:3619–27.

    Article  CAS  Google Scholar 

  40. Guillemi EC, Tomassone L, Farber MD. Tick-borne Rickettsiales: molecular tools for the study of an emergent group of pathogens. J Microbiol Methods. 2015;119:87–97.

  41. Duron O, Binetruy F, Noel V, Cremaschi J, McCoy KD, Arnathau C, et al. Evolutionary changes in symbiont community structure in ticks. Mol Ecol. 2017;26:2905–21.

    Article  CAS  Google Scholar 

  42. Tsementzi D, Castro Gordillo J, Mahagna M, Gottlieb Y, Konstantinidis KT. Comparison of closely related, uncultivated Coxiella tick endosymbiont population genomes reveals clues about the mechanisms of symbiosis. Environ Microbiol. 2018;20:1751–64.

    Article  CAS  Google Scholar 

  43. Ge Y, Guo G, Ge B, Yin H, Yin H. The spleen microbiota of small wild mammals reveals distinct patterns with tick-borne bacteria. PLoS Negl Trop Dis. 2018;12:e0006499.

    Article  Google Scholar 

  44. Duron O, Noel V, McCoy KD, Bonazzi M, Sidi-Boumedine K, Morel O, et al. The recent evolution of a maternally-inherited endosymbiont of ticks led to the emergence of the Q fever pathogen, Coxiella burnetii. PLoS Pathog. 2015;11:e1004892.

    Article  Google Scholar 

  45. Spitalska E, Sparagano O, Stanko M, Schwarzova K, Spitalsky Z, Skultety L, Havlikova SF. Diversity of Coxiella-like and Francisella-like endosymbionts, and Rickettsia spp., Coxiella burnetii as pathogens in the tick populations of Slovakia, Central Europe. Ticks Tick Borne Dis. 2018;9:1207–11.

    Article  Google Scholar 

  46. Güner ES, Hashimoto N, Kadosaka T, Imai Y, Masuzawa T. A novel, fast-growing Borrelia sp. isolated from the hard tick Hyalomma aegyptium in Turkey. Microbiology. 2003;149:2539–44.

    Article  Google Scholar 

  47. Güner ES, Watanabe M, Hashimoto N, Kadosaka T, Kawamura Y, Ezaki T, et al. Borrelia turcica sp. nov., isolated from the hard tick Hyalomma aegyptium in Turkey. Int J Syst Evol Microbiol. 2004;54:1649–52.

    Article  Google Scholar 

  48. Takano A, Goka K, Une Y, Shimada Y, Fujita H, Shiino T, et al. Isolation and characterization of a novel Borrelia group of tick-borne borreliae from imported reptiles and their associated ticks. Environ Microbiol. 2010;12:134–46.

    Article  CAS  Google Scholar 

  49. Takano A, Fujita H, Kadosaka T, Konnai S, Tajima T, Watanabe H, et al. Characterization of reptile-associated Borrelia sp. in the vector tick, Amblyomma geoemydae, and its association with Lyme disease and relapsing fever Borrelia spp. Environ Microbiol Rep. 2011;3:632–7.

    Article  Google Scholar 

  50. Cutler SJ, Ruzic-Sabljic E, Potkonjak A. Emerging borreliae - expanding beyond Lyme borreliosis. Mol Cell Probes. 2017;31:22–7.

    Article  CAS  Google Scholar 

  51. Kalmar Z, Cozma V, Sprong H, Jahfari S, D’Amico G, Marcutan DI, et al. Transstadial transmission of Borrelia turcica in Hyalomma aegyptium ticks. PLoS One. 2015;10:e0115520.

  52. Angelakis E, Mediannikov O, Jos SL, Berenger JM, Parola P, Raoult D. Candidatus Coxiella massiliensis infection. Emerg Infect Dis. 2016;22:285–8.

    Article  CAS  Google Scholar 

  53. Seo MG, Lee SH, VanBik D, Ouh IO, Yun SH, Choi E, et al. Detection and genotyping of Coxiella burnetii and Coxiella-like bacteria in horses in South Korea. PLoS One. 2016;11:e0156710.

    Article  Google Scholar 

  54. Shivaprasad HL, Cadenas MB, Diab SS, Nordhausen R, Bradway D, Crespo R, et al. Coxiella-like infection in psittacines and a toucan. Avian Dis. 2008;52:426–32.

    Article  CAS  Google Scholar 

  55. Woc-Colburn AM, Garner MM, Bradway D, West G, D’Agostino J, Trupkiewicz J, et al. Fatal coxiellosis in Swainson’s Blue Mountain rainbow lorikeets (Trichoglossus haematodus moluccanus). Vet Pathol. 2008;45:247–54.

    Article  CAS  Google Scholar 

  56. Tanyel E, Guler N, Hokelek M, Ulger F, Sunbul M. A case of severe babesiosis treated successfully with exchange transfusion. Int J Infect Dis. 2015;38:83–5.

    Article  Google Scholar 

  57. Altay K, Aktas M, Dumanli N. Detection of Babesia ovis by PCR in Rhipicephalus bursa collected from naturally infested sheep and goats. Res Vet Sci. 2008;85:116–9.

    Article  CAS  Google Scholar 

  58. Aydin MF, Aktas M, Dumanli N. Molecular identification of Theileria and Babesia in ticks collected from sheep and goats in the Black Sea region of Turkey. Parasitol Res. 2015;114:65–9.

    Article  Google Scholar 

  59. Aktas M, Ozubek S. A survey of canine haemoprotozoan parasites from Turkey, including molecular evidence of an unnamed Babesia. Comp Immunol Microbiol Infect Dis. 2017;52:36–42.

    Article  Google Scholar 

  60. Liu Y, Zhou J, White KP. RNA-seq differential expression studies: more sequence or more replication? Bioinformatics. 2014;30:301–4.

    Article  CAS  Google Scholar 

  61. van Dijk EL, Jaszczyszyn Y, Thermes C. Library preparation methods for next-generation sequencing: tone down the bias. Exp Cell Res. 2014;322:12–20.

    Article  Google Scholar 

  62. Brinkmann A, Ergünay K, Radonic A, Kocak Tufan Z, Domingo C, Nitsche A. Development and preliminary evaluation of a multiplexed amplification and next generation sequencing method for viral hemorrhagic fever diagnostics. PLoS Negl Trop Dis. 2017;11:e0006075.

    Article  Google Scholar 

Download references


The authors are grateful to Ursula Erikli for meticulous copy-editing and N. Emin Güven for graphical support. Preliminary findings of this study have been presented at the IMED - International Meeting on Emerging Diseases and Surveillance, held November 9–12, 2018 in Vienna, Austria.


KE is a recipient of the Georg Forster Research Fellowship (HERMES) for Experienced Researchers by the Alexander von Humboldt Foundation, 2015. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

Availability of data and materials

Data supporting the conclusions of this article are included within the article and its additional files. The newly generated nucleotide sequences were deposited in the GenBank database under the accession numbers: MH645175-645181, MH645183-645205, MH618772-618775, MH618686, MH628249, MH628250, MH630144-630147, MH649268, MH649269, MH673722 and MH673723. The HTS data were deposited within the Sequence Read Archive (SRA) of the United States Library of Medicine National Center for Biotechnology Information website under the accession numbers SRX4641594-4641620.

Author information

Authors and Affiliations



Specimen collection and identification: ED and OH. Laboratory assays: AB, PH and ED. Data interpretation: AB and KE. Project planning, general overview and manuscript preparation: KE, AN and PH. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Koray Ergünay.

Ethics declarations

Ethics approval and consent to participate

The study involved testing of ticks, collected in the field during questing or from domesticated animals and no institutional or regional ethics committee approval was required. Removal of ticks from infested domestic animals was carried out with the informed consent and cooperation of the caretakers and/or owners.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Tick pools evaluated using high throughput sequencing. Table S2. Microbial detection rates in pools according to tick species. (XLS 38 kb)

Additional file 2:

Near-complete R. hoogstraalii genome sequence, assembled from the Ha. parva pool P4. The sequence is available in FASTA format and observed gaps following alignment to the R. hoogstraalii strain Croatica genome (CCXM01000001) are indicated. (TXT 1410 kb)

Additional file 3:

Alignment of the partial conjugative plasmid sequences of R. hoogstraalii p4 characterized in this study (GenBank: MH649269), with R. hoogstraalii strain Croatica (CCXM01000002), R. felis strain (CP000054) and R. australis strain Cutlack (CP003339). (PDF 63 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Brinkmann, A., Hekimoğlu, O., Dinçer, E. et al. A cross-sectional screening by next-generation sequencing reveals Rickettsia, Coxiella, Francisella, Borrelia, Babesia, Theileria and Hemolivia species in ticks from Anatolia. Parasites Vectors 12, 26 (2019).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: