Identification of tick-borne pathogens by metagenomic next-generation sequencing in Dermacentor nuttalli and Ixodes persulcatus in Inner Mongolia, China

Hard ticks act as arthropod vectors in the transmission of human and animal pathogens and are widely distributed in northern China. The aim of this study is to screen the important tick-borne pathogens (TBPs) carried by hard ticks in Inner Mongolia using metagenomic next-generation sequencing (mNGS) and to estimate the risk of human infection imposed by tick bites. The adult Dermacentor nuttalli (n = 203) and Ixodes persulcatus (n = 36) ticks feeding on cattle were collected. The pooled DNA samples prepared from these ticks were sequenced as the templates for mNGS to survey the presence of TBPs at the genus level. Individual tick DNA samples were detected by genus--specific or group-specific nested polymerase chain reaction (PCR) of these TBPs and combined with DNA sequencing assay to confirm the results of mNGS. R. raoultii (45.32%, 92/203), Candidatus R. tarasevichiae (5.42%, 11/203), Anaplasma sp. Mongolia (26.60%, 54/203), Coxiella-like endosymbiont (CLE) (53.69%, 109/203), and Babesia venatorum (7.88%, 16/203) were detected in D. nuttalli, while R. raoultii (30.56%, 11/36), Anaplasma sp. Mongolia (27.80%, 10/36), and CLE (27.80%, 10/36) were detected in I. persulcatus. The double- and triple-pathogen/endosymbiont co-infections were detected in 40.39% of D. nuttalli and 13.89% of I. persulcatus, respectively. The dual co-infection with R. raoultii and CLE (14.29%, 29/203) and triple co-infection with R. raoultii, Anaplasma sp. Mongolia, and CLE (13.79%, 28/203) were most frequent in D. nuttalli. This study provides insight into the microbial diversity of D. nuttalli and I. persulcatus in Inner Mongolia, China, reporting for the first time that Candidatus R. tarasevichiae had been found in D. nuttalli in China, and for the first time in the world that Anaplasma sp. Mongolia has been detected in I. persulcatus. This study proves that various vertically transmitted pathogens co-inhabit D. nuttalli and I. persulcatus, and indicates that cattle in Inner Mongolia are exposed to several TBPs.


Background
Hard ticks (Acari: Ixodidae) are obligate blood-sucking parasitic arthropods which can infest mammals, birds, and reptiles, and act as arthropod vectors in the transmission of human and animal pathogens. A wide variety of pathogens can be maintained and transmitted by hard ticks, including Ehrlichia spp., Anaplasma spp., Rickettsia spp., Coxiella spp., Babesia spp., Borrelia spp., etc. [1]. In addition, a variety of endosymbionts, such as Coxiella-like, Rickettsia-like, and Arsenophonus-like endosymbionts, live inside hard ticks [2][3][4]. Therefore, hard ticks are usually considered to be the most important vectors of pathogens, and knowledge of the microbial communities within these ticks will be of benefit for risk assessment of tick-borne diseases.
Mongolia Hulunbuir League, one of the important pastoral regions in Inner Mongolia in China, contains significant amounts of pastures and is an important region for animal production. The eastern part of Hulunbuir, stretching across the primeval forest-covered area in the Daxing'anling Mountains, is an important habitat for hard ticks [15]. Many TBPs including C. burnetii [16], Rickettsia spp. [16][17][18], Anaplasma spp. [19], and tick-borne encephalitis virus [20] have been detected in hard ticks collected here. However, little attention has been given to co-infection with TBPs in ticks, and continued research is needed to fully comprehend the diversity of TBPs. In this study, we investigated the microbial communities in hard ticks collected from cattle in Hulunbuir League to reveal the coexistence of TBPs using metagenomic next-generation sequencing (mNGS) combined with nested polymerase chain reaction (PCR). The results of our study might provide broader knowledge of the microorganisms inside hard ticks in the region, thereby strengthening programs to prevent and control the potential infections caused by TBPs.

Tick washing, homogenization, and DNA extraction
To remove environmental contaminants, each tick was surface-sterilized twice with 75% ethanol, followed by phosphate-buffered saline (PBS) twice. Ticks were then individually homogenized in 300 μL of PBS using MagNA Lyser Green Beads (Roche, Mannheim, Germany), and DNA extraction was performed on 200 μL of each tick homogenate using a QIAamp ®  the manufacturer's instructions. The extracted genomic DNA was dissolved in 100 μL ultrapure water and stored at −20 °C for further analysis. For these previous steps, ultrapure water, sterile tubes, and filter tips were used, and all operations were carried out in a biological safety cabinet. Each time DNA extraction was performed, an extraction control (water) was added.
Individual DNA samples were mixed in an equal volume (20 μL) to prepare pooled DNA samples for full microbial genome sequencing using mNGS.

Metagenome assembly, gene prediction, and taxonomy prediction
All pooled DNA samples were paired-end sequenced on the Illumina HiSeq platform (insert size 350 bp, read length 150 bp) by the Beijing Genomics Institute (BGI) (Beijing, China). The reads with more than 40 nt lowquality bases (quality value ≤ 38) were removed. Meanwhile, the reads with more than 10 nt "N" bases were filtered out of the data sets. Lastly, the reads overlapping more than 15 nt bases with the adapters were removed. Reads that aligned to tick genes were also removed using Bowtie 2 (v2.2.4) with the parameters -end-to-end, -sensitive, -I 200, -X 400 [22,23]. Accordingly, the clean data were obtained.
Then the clean reads were mapped against scaffolds using SOAPdenovo (V2.04) with the parameters -d 1, -M 3, -R, -u, -F, -K 55 [24]. The unused reads from each sample were then assembled using the same parameters. The scaffolds were broken at N into the scaftigs [25], and the scaftigs with the length of ≥ 500 nt were used for further analysis [26]. Open reading frames (ORFs) in the scaftigs (≥ 500 bp) were predicted by MetaGeneMark (V2.10) [23,27]. A nonredundant gene catalog was obtained after processing by using CD-HIT (V4.5.8) with the parameters -c 0.95, -G 0, -aS 0.9, -g 1, -d 0 [28,29], and using a sequence identity threshold of 0.95 and a minimum coverage cutoff of 0.9. To determine the gene abundances, the reads were realigned with the gene catalog using Bowtie 2 and the following parameters: -m 200 -× 400 -s 119. Only genes with ≥ 2 mapped reads were deemed to be present in a sample [30]. Relative abundance of genes was calculated based on the number of reads mapped to the genes and the length of the genes as previously described [31][32][33].
To access the taxonomic assignments of genes, genes were aligned to the integrated NR database (version: 2018-01-02) of NCBI using DIAMOND (V0.9.9) and default parameters, with the exception of -k 50 -sensitive -e 0.00001 [34]. For each gene, the significant matches which were defined by e-values ≤ 10 * e-value of the top hit were retained to distinguish taxonomic groups [30]. Then the taxonomical level of each gene was determined by using the lowest common ancestor (LCA)-based algorithm implemented in MEGAN [35]. The results containing the number of genes and the abundance information for each sample, and the relative abundances of each taxonomic group were calculated by adding the relative abundances of genes annotated to the same feature [23,26,36].

Polymerase chain reaction (PCR)
Based on the results of mNGS, genus-/group-specific PCR was performed to confirm the presence of TBPs in individual ticks. PCR was performed using a PCR System 9700 (Applied Biosystems, GeneAmp ® , USA). For nested PCR, 1 μL of each individual DNA sample (150-330 ng) was used as template for the first round, and 1 μL of the primary PCR production was used as template for the second round. For the first round, a negative control (water) and an extraction control mentioned above were included in each PCR experiment. Tube strips with individual caps were used in amplification steps to prevent cross-contamination, and all PCR amplifications were carried out using PrimeSTAR ® HS (Premix) (TaKaRa, Beijing, China). All operations were carried out in a biological safety cabinet. Amplified products were then electrophoresed in 1.5% agarose gel, and the positive amplicons were sent to TSINGKE Biological Technology (Beijing, China) for sequencing. The PCR primers for the spotted fever group Rickettsia (SFGR) [37], Anaplasma spp. and Ehrlichia spp. [38], Coxiella spp. [39], and Babesia spp. [40] are presented in Table 1.

Phylogenetic analysis
The obtained DNA sequences were compared with those available in GenBank using the National Center for Biotechnology Information (NCBI; Bethesda, MD) Basic Local Alignment Search Tool (BLAST) search engine (http:// blast. ncbi. nlm. nih. gov/ blast. cgi), and multiple sequence alignment was performed using the ClustalW multiple alignment tool with the default parameters in MEGA 7.0. The phylogenetic analysis of gltA for SFGR, 16S rRNA for Anaplasma spp., 16S rRNA for Coxiella spp., or 18S rRNA for Babesia spp. was performed using the maximum likelihood method based in MEGA 7.0. Bootstrap values were estimated for 1000 replicates [41,42].

Taxonomic classification
A total of 239 adult hard ticks were identified as D. nuttalli (n = 203) (accession number: MK213083.1) and I. persulcatus (n = 36) (accession number: MH790201.1) based on morphological identifications confirmed by species-specific PCR and sequencing assay. Ten pools of D. nuttalli DNA samples were finally analyzed by mNGS on the Illumina HiSeq platform. Sequencing yielded between 5970 and 7475 million reads per pool, while all were of high quality (Clean_Q20 > 96%) (shown in Additional file 1: Table S1). The construction of a metagenomic library of I. persulcatus DNA samples failed.
The presence of the bacterial genera Rickettsia, Anaplasma, and Coxiella in the pooled tick samples was confirmed by the taxonomic profiles at genus level ( Fig. 2; Table 2). Rickettsia spp. were most abundant in sample pools 2-4 and 7-8 and also abundant in other pools. In pools 1, 9, and 10, Anaplasma spp. were abundant. However, Coxiella spp. were abundant only in pool 1. In addition, Pseudomonas spp. were most abundant in pools 5   and 6, and Psychrobacter spp. were most abundant in pool 10 ( Table 2).

Prevalence of tick-borne pathogens in individual ticks
By mNGS, the important pathogenic bacterial genera Rickettsia, Anaplasma, and Coxiella were found in the pooled tick samples, and thus each tick was detected by the genus-/group-specific PCR combined with sequencing in order to identify the TBPs carried by it. In addition, Babesia spp. were often detected in ticks, and thus each tick was detected by Babesia-specific PCR.  (Table 3).

Co-infection in individual ticks
In 190 TBP-positive ticks, 87 ticks (45.79%) were found to be co-infected with more than one species identified in the present study ( Table 3). The dual-and triplepathogen/endosymbiont co-infections were detected in 40.39% of D. nuttalli and 13.89% of I. persulcatus. The dual co-infection with R. raoultii and CLE and the triple co-infection with R. raoultii, Anaplasma sp. Mongolia and CLE were most frequent in D. nuttalli (Table 3).

Discussion
In recent years, much attention has been focused on ticks and TBPs in China. However, although a variety of pathogens have been identified, co-infection with multiple pathogens in hard ticks has rarely been investigated. In this study, we applied mNGS combined with nested PCR to survey TBPs in D. nuttalli and I. persulcatus feeding on cattle in Inner Mongolia, China. By mNGS, the endosymbionts including Coxiella spp., Rickettsia spp., Francisella spp., and "Candidatus Midichloria mitochondrii" have been recognized as the most abundant bacterial species identified frequently in entirely homogenized ticks [43]. Qiu et al. applied NGS to examine the microbiomes of salivary glands of ticks collected in Japan, revealing a large number of bacterial genera, including 71 I. ovatus, 127 I. persulcatus, and 59 H. flava, and detecting some of the medically important bacteria including Coxiella spp., Ehrlichia spp., and Rickettsia spp. [44].
In the present study, the pooled DNA samples of D. nuttalli and I. persulcatus collected were assayed by mNGS. The result revealed the presence of the bacterial genera Rickettsia, Anaplasma, and Coxiella in these ticks. In order to identify the bacteria at species level, each tick was detected by SFGR-, Anaplasma and Ehrlichia-, and Coxiella-specific PCR as well as Babesia-specific PCR, respectively. After sequencing of the DNA fragments amplified by PCR and the sequence comparison, two Rickettsia species (R. raoultii and Candidatus R. tarasevichiae), one Anaplasma species (Anaplasma sp. Mongolia), B. venatorum, and CLEs were found in D. nuttalli, while R. raoultii, Anaplasma sp. Mongolia, and CLEs were also found in I. persulcatus.
R. raoultii, a species of SFGR, was firstly detected in Dermacentor ticks collected in Russia in 1999 and isolated from Dermacentor ticks and named in 2008 [45]. It is one of the causative agents of tick-borne lymphadenopathy (TIBOLA), which is also known as Dermacentor-borne necrosis erythema and lymphadenopathy (DEBONEL) in humans [46]. R. raoultii has been found to be present in various ticks, including Dermacentor, Haemaphysalis, Rhipicephalus, Hyalomma, and Amblyomma. In the present study, R. raoultii was detected in 45.32% of D. nuttalli and 30.56% of I. persulcatus, suggesting that it was the dominant Rickettsia species prevalent in the hard ticks in Inner Mongolia, and this may  Candidatus R. tarasevichiae, an emerging tick-borne pathogen, is also a species of SFGR. It was first detected in I. persulcatus collected from the southern Urals and Siberia in 2003 [47] and then found in Haemaphysalis ticks in Far East regions in Russia [48]. Human cases caused by Candidatus R. tarasevichiae have been found in China and Russia [49,50]. In this study, Candidatus R. tarasevichiae was detected in 5.42% of D. nuttalli, which was most related to the Candidatus R. tarasevichiae Mulan-11 strain (MN450396.2) and Bayan-68 strain (MN450397.2) from I. persulcatus in China in phylogenetic analysis. This is the first time that Candidatus R. tarasevichiae has been detected in D. nuttalli.
Anaplasma sp. Mongolia was firstly detected in D. nuttalli [51] and bovine blood in Mongolia [51,52], demonstrating that the Anaplasma species is an important cattle pathogen. In this study, Anaplasma sp. Mongolia was detected in 26.6% of D. nuttalli and in 27.8% of I. persulcatus. This study is thus the first in the world to report the presence of Anaplasma sp. Mongolia in I. persulcatus.
CLEs are relatively common in the microbiota of various tick species around the world, forming multiple subclusters in the cluster of the genus Coxiella in phylogenetic analysis [39,53]. The presence of these symbiotic bacteria in ticks confers crucial and diverse benefits to the host, affecting its development, nutrition, chemical defense, or reproduction [53][54][55]. The prevalence of CLEs is from 6.25% in Rhipicephalus sanguineus to 100% in Amblyomma americanum in North America and Europe [53]. In the present study, CLEs were detected in 53.69% of D. nuttalli and 27.78% of I. persulcatus. Phylogenetic analysis suggested that the CLE strain of D. nuttalli, which was mostly related to that (KP994814.1) from D. silvarum in France, was different from that of I. persulcatus, which was mostly related to that (KM079619.1) of Haemaphysalis concinna from Russia.
Babesia spp. are the pathogenic agents of babesiosis in humans and animals. In the present study, B. venatorum was detected in 7.88% of D. nuttall. However, Babesia spp. was not found in the pooled tick samples by mNGS assay, which might be caused by the extremely low abundance of Babesia spp. in the pooled samples. According to phylogenetic analysis, B. venatorum of D. nuttalli  [56].
In the present study, the multiple pathogen/endosymbiont co-infections were detected in 40.39% of D. nuttalli and 13.89% of I. persulcatus. Ticks may acquire multiple pathogenic species during blood feeding on their vertebrate hosts, and the hosts may also be infected by the pathogens carried by ticks. Due to the development of molecular diagnostic methods, more and more cases with co-infection of multiple TBPs have been reported [57]. The co-infection may be the result of a single tick bite by the tick carrying more than one pathogen or the result of multiple bites by ticks carrying different pathogens. Therefore, the prevalence of co-infection with TBPs in people living in the area close to the natural focus of TBPs should be investigated in the future.

Conclusions
This study proves that various vertically transmitted pathogens co-inhabit D. nuttalli and I. persulcatus, and is the first to report that Candidatus R. tarasevichiae has been found in D. nuttalli in China, and the first in the world to report that Anaplasma sp. Mongolia has been detected in I. persulcatus. This study provides insight into the microbial diversity of D. nuttalli and I. persulcatus in Inner Mongolia, China, and indicates that cattle in Inner Mongolia are exposed to several TBPs.