Metagenomic profiles of Dermacentor tick pathogens from across Mongolia, using next generation sequencing.

Tick-borne diseases are a major public health concern in Mongolia. Nomadic pastoralists, which make up ~ 26% of Mongolia's population, are at an increased risk of both tick bite exposure and economic loss associated with clinical disease in herds. This study sought to further characterize tick-borne pathogens present in Dermacentor ticks (n = 1,773) sampled in 2019 from 15 of Mongolia's 21 aimags (provinces). The ticks were morphologically identified and sorted into 377 pools which were then screened using Next-Generation Sequencing paired with confirmatory PCR and DNA sequence analysis. Rickettsia spp. were detected in 88.33% of pools, while Anaplasma spp. and Bartonella spp. were detected in 3.18 and 0.79% of pools, respectively. Khentii had the highest infection rate for Rickettsia spp. (76.61%; CI: 34.65-94.79%), while Arkhangai had the highest infection rate for Anaplasma spp. (7.79%; CI:4.04-13.72%). The exclusive detection of Anaplasma spp. in tick pools collected from livestock supports previous work in this area that suggests livestock play a significant role in disease maintenance. The detection of Anaplasma, Bartonella, and Rickettsia demonstrates a heightened risk for infection throughout Mongolia, with this study, to our knowledge, documenting the first detection of Bartonella melophagi in ticks collected in Mongolia. Further research deploying NGS methods is needed to characterize tick-borne pathogens in other endemic tick species found in Mongolia, including Hyalomma asiaticum and Ixodes persulcatus.

Introduction

Ticks and the pathogens they carry pose a significant threat to both human and animal health. This holds true in Mongolia, where an estimated 26% of the population continues to live a nomadic pastoral lifestyle and 37% of households own livestock (Odontsetseg et al., 2009; Boldbaatar et al., 2017; Barnes et al., 2020). These populations spend prolonged periods of time moving herds through tick habitats, resulting in a heightened risk for exposure to ticks and tick-borne diseases (TBDs). The Mongolian economy is also likely impacted by the effects of TBDs, where an estimated 67 million heads of livestock are present within the country and roughly 18% of the nation’s GDP comes from animal-related products (Odontsetseg et al., 2009). In neighboring China, an estimated $70 million every year is lost due to the impact of tick-borne disease impacts on small mammal production (Yin and Luo, 2007).

Ticks gathered in Mongolia have previously tested positive for various TBDs, including Anaplasma spp., Borrelia spp., Crimean-Congo hemorrhagic fever, Ehrlichia spp., Rickettsia spp., and tick-borne encephalitis virus (Moore et al., 2018; Voorhees et al., 2018; Černý et al., 2019; von Fricken et al., 2020a). Rickettsial diseases are of particular concern due to high rates of severe illness and death in previously healthy individuals, (Aung et al., 2014; Biggs, 2016; von Fricken et al., 2018). A previous study by our team found that 20% of humans and livestock animals in Mongolia have had past exposure to Rickettsia spp., with variations observed by geographic location (von Fricken et al., 2018). This also held true when examining previous exposure to Anaplasma spp., which was detected in 37% of nomadic herders and over 40% of livestock (von Fricken et al., 2018). We also have detected Anaplasma ovis infection rates as high as 80% in sheep and 69% in goats, which aligns with what has previously been detected in ticks from the same region (Ochirkhuu et al., 2017; von Fricken et al., 2018, 2020a; Enkhtaivan et al., 2019; Fischer et al., 2020; Chaorattanakawee et al., 2022). Anaplasmosis in livestock can result in anoxia, abortions, infertility, significant weight loss, and even death, all of which can impact economic security in pastoralist communities.

Dermacentor ticks are the most common and one of the more important ticks of medical and veterinary concern within Mongolia due to their wide geographic range and the pathogens they carry (Černý et al., 2019). Ticks collected from southern and central aimags have previously had high pool positivity rates (> 80%) for Rickettsia spp., with molecular detections of R. raoultii, R. sibirica mongolitimonae, and R. sibirica reported (Fischer et al., 2020; von Fricken et al., 2020b). In contrast, a study of pathogens within ticks collected from aimags of central Mongolia found lower overall levels of Anaplasma spp. within Dermacentor ticks, although the infectivity rates increased substantially when specifically examining ticks removed from livestock (von Fricken et al., 2020a). Additional pathogens have been detected within Dermacentor ticks collected from Mongolian aimags include Babesia caballi, B. equi, Borrelia afzelii, Candidatus Midichloria sp., Candidatus Neoehrlichia mukurensis, Theileria equi, and T. orientalis (Battsetseg et al., 2001; Javkhlan et al., 2014; Fischer et al., 2020). In neighboring countries, pathogens reported from Dermacentor spp. ticks include Babesia venatorum, Borrelia miyamotoi, Brucella spp., Francisella tularensis subsp. holarctica, Rickettsia aeschlimannii, and the Far Eastern genotype of tick-borne encephalitis virus (Zhang et al., 2008; Wei et al., 2016; He et al., 2018; Yin et al., 2018; Huang et al., 2020; Gao et al., 2021; Jiao et al., 2021).

The potential threat tick-borne diseases present to both the Mongolian population and its growing ecotourism industry is substantial, given the high rates of various pathogens reported in previous tick survey studies (Černý et al., 2019; Fischer et al., 2020; von Fricken et al., 2020b). Improved molecular characterization of TBDs within Mongolia may help inform future preventative measures for locals and visitors, while also establishing baseline sequence data to monitor evolution over time. The variety of tick-borne pathogens found within Mongolia complicates attempts to fully characterize pathogens found in samples collected within the country. Our research group has recently used an analytical workflow on livestock blood samples from three aimags in Mongolia, initially applying next-generation sequencing (NGS) to obtain a snapshot of pathogen groups present, followed by conventional PCR and Sanger sequencing for confirmation and species characterization (Chaorattanakawee et al., 2022). In this study, we deploy Next-Generation Sequencing for on Dermacentor ticks collected from a wide geographic range of Mongolia to further our understanding of tick-borne pathogens in this region.

Materials and methods

Dermacentor ticks were collected from the environment (questing) and off domestic animals from 15 aimags across Mongolia in 2019 (Uvs, Khovd, Govi-Altai, Zavhan, Khuvsgul, Arkhangai, Bayankhongor, Arkhangai, Uvurkhangai, Bulgan, Tuv, Dundgovi, Khentii, Dornogovi, Sukhbaatar, and Dornod; Figure 1). Adult ticks were morphologically identified as D. nuttalli or to the genus level as Dermacentor spp. by entomologists using local keys (Boldbaatar and Byambaa, 2015). In total, 7,275 ticks were collected and sorted into 1,489 pools according to location and collection source (environment vs. animal). Of these pools, 377 pools of adult stage ticks, representing pools from all sampled provinces, were selected for analysis by next-generation sequencing, including 51 pools collected from livestock (Tables 1, 2). Whole ticks in 250 μl of ATL buffer were punctured with a fine tip under a stereomicroscope to release the tissue from the hard chitin exoskeleton prior to adding 2 mg/ml of Proteinase K solution. Samples were then incubated at 55°C overnight. A total volume of 250 μl homogenized solution was then used for DNA extraction on the QIAsymphony® SP instrument with QIAsymphony® DSP DNA Mini Kit using Tissue LC 200 DSP protocol (Qiagen, Hombrechtikon, Switzerland). The DNA was eluted in 50 μl of ATE buffer and stored at −20°C until use.

Figure 1

Distribution of tick collection events symbolized according to collection source. A map of Mongolia representing the location of tick pools that were chosen for further analysis. Tick pools are symbolized according to their collection source off animals (brown circle) or from the environment (green square). Individual aimags are colored according to the number of pools analyzed per aimag to demonstrate sampling intensity.

Table 1

Samples selected for pathogen screening by next-generation sequencing (NGS) in this study.

Provinces	Number of tick pools	Number of ticks	Number of tick pools selected for NGS	Number of ticks in pools selected for NGS
Arkhangai	61	305	27	135
Bayankhongor	101	496	29	143
Bulgan	14	56	14	56
Dornod	122	599	34	159
Dornogovi	5	18	5	18
Dundgovi	47	226	27	126
Govi-Altai	135	675	24	120
Khentii	8	17	8	17
Khovd	162	810	26	130
Khuvsgul	223	1,050	51	232
Sukhbaatar	140	693	27	133
Tuv	4	21	4	21
Uvs	162	810	29	145
Uvurkhangai	124	594	42	188
Zavhan	181	905	30	150
Total	1,489	7,275	377	1773

Table 2

Pool positivity rate by collection source [% and (95% CI)].

Pathogens	Sources (% Infection, 95% CI)
	Animals (N = 51)	Environment (N = 307)	Rock and bush (N = 19)	Total
Rickettsia	41 [80.4% (69.5, 91.3%)]	275 [89.6% (86.2, 93.0%)]	17 [89.5% (75.7, 103.0%)]	333
Anaplasma	12 [23.5% (11.9, 35.2%)]	0	0	12
Bartonella	1 [1.96% (−1.8, 5.8%)]	2 [0.7% (−0.3, 1.6%)]	0	3

Bacterial 16S DNA amplification

Nested PCR was performed as described in Chaorattanakawee et al., 2022 to amplify both the V1-V6 region and V3-V4 region of the bacterial 16S rDNA. Each round of PCR included both an ultrapure DNA/RNA-free water negative control and a mock DNA extraction control. The nested PCR amplicon products were isolated using AMPure XP magnetic beads and the quality of the products was assessed as previously described (Chaorattanakawee et al., 2022). Amplicon products were stored at −20°C until further analysis.

Library preparation and sequencing

The Nextera XT Index Kit v2 (Illumina) was used for index PCR to attach the dual indices and Illumina sequencing adapters to purified 16S amplicons as previously described. Each batch of indexing reactions included a DNA/RNA-free water as a negative control. The index PCR products were cleaned using AMPure XP beads, followed by library purity analysis using the QIAxcel Advanced System (Qiagen). The index libraries were then quantified using the Qubit dsDNA HS Assay Kit (Invitrogen). Libraries were denatured with NaOH according to the manufacturer’s protocol (Ilumina). Sequencing was performed using the MiSeq Reagent Kit V3 with the Ilumina Miseq System. A 10% PhiX internal control (Ilumina) was included in each low-diversity library run.

NGS data analysis

Sequence reads produced by the Ilumina MiSeq system were processed using the CLC Genomics workbench (v 11.0.1) and CLC microbial genomics module (v 3.0; Qiagen, Aarhus A/S1), which included merging paired reads, primer sequence removal, low read sample removal, and chimeric sequence removal. The filtered sequences were then clustered into operational taxonomic unites (OTUs) using a threshold of 97% sequence identity and the reference OTU database downloaded from the Greengenes database (v 13.8) and SILVA 16S (v 132). Pathogen reads detected in the negative controls represented cross-contamination and were used to subtract respective reads detected in samples.

Pathogen characterization by PCR and sanger sequencing

To confirm the detection of pathogens and the taxonomic assignment as indicated by NGS analysis, PCR and DNA sequencing were conducted on NGS samples with read counts above a set threshold. The assays and gene targets for selected pathogens (Anaplasma, Bartonella, Rickettsia, Coxiella) were detailed previously in Takhampunya et al. (2019). PCR amplification products were cleaned using the ExoSAP-IT kit (Applied Biosystem), followed by cycle-sequencing and sequencing using the SeqStudio Genetic Analyzer (Applied Biosystems), as previously described (Chaorattanakawee et al., 2022). Sample sequences were assembled using Sequencher (v 5.1, Gene Codes Corp.) and aligned with GenBank reference sequences using the MUSCLE codon alignment program. Maximum likelihood phylogenetic trees were constructed for each bacterial target gene using MEGA 6.

Mapping

ArcGIS Pro (v 2.8.0, ESRI) was used for spatial visualizations of data, including tick collection events, tick collection source, and pathogen detection. The map layer of Mongolia and its delineated aimags was accessed from ESRI.

Statistical analysis

To estimate the probability of pathogen detection within the pooled samples, prevalence rates, maximum likelihood estimates (MLE) and minimum infection rates (MIR) were calculated, which is standard when analyzing pooled tick data. The MLE and MIR estimates were conducted in Excel with the use of the CDC’s Mosquito Surveillance Software tool which calculate point and confidence intervals using pooled data that take into account individual pool sample sizes to estimate infection rates.

Results

Detection of Rickettsia, Anaplasma, and Bartonella

The summary results for the Rickettsia spp., confirmed through qPCR analysis and DNA sequencing, are presented in Table 3 and Figure 2 Overall, Rickettsia spp. were detected in tick pools from all aimags sampled, with 88% of pools testing positive (333/377). The highest Rickettsia spp. pool detection rate was seen in Tuv (100%) followed by: Dornod (97%) Dundgovi (96%) and Sukhbaatar (96%), while the Bulgan aimag showed the lowest pool positivity rate (57%). Maximum likelihood estimates (MLE) found an average prevalence of 37.30% (95%CI: 33.50–41.01%), where Dornod aimag had the highest MLE of 55.40% (95%CI: 34.38–70.67%) and a MIR of 20.75% (95%CI: 14.45–27.06%) and Bulgan had the lowest MLE of 18.34% (95%CI: 9.21–30.23%) with a MIR of 14.29% (95%CI: 5.12–23.45%). In general, higher MLEs were found in tick pools collected from eastern and western aimags of Mongolia, with lower MLEs seen in central aimags (Figure 2).

Table 3

Maximum likelihood estimates of Rickettsia spp. by region based on confirmatory results including 95% confidence intervals.

Rickettsia spp.
Province	Positive pools	Total number of ticks	MLE
			Point	Low	High
Arkhangai	20/27 (74%)	135	23.7	14.9	33.3
Bayankhongor	24/29 (83%)	143	31.2	20.0	42.4
Bulgan	8/14 (57%)	56	18.3	9.2	30.2
Dornod	33/34 (97%)	159	55.4	34.4	70.7
Dornogovi	4/5 (80%)	18	33.4	12.1	54.9
Dundgovi	26/27 (96%)	126	51.5	31.1	67.1
Govi-Altai	22/24 (92%)	120	39.2	23.7	52.9
Khentii	7/8 (88%)	17	76.6	34.7	94.8
Khovd	24/26 (92%)	130	40.1	24.7	53.7
Khuvsgul	41/51 (80%)	232	31	22.6	39.4
Sukhbaatar	26/27 (96%)	133	48.7	29.2	63.9
Tuv	4/4 (100%)	21	N/A	N/A	N/A
Uvs	27/29 (93%)	145	41.4	26.2	54.7
Uvurkhangai	39/42 (93%)	188	43.9	31	55.1
Zavhan	28/30 (93%)	150	41.8	26.6	55.0
Total	333/377 (88%)	1773	37.3	33.5	41.0

Figure 2

Pool MLE for Rickettsia spp. A map showing the distribution of identified Rickettsia species, with an aimag color gradient representing the Rickettsia spp. MLE of sampled pools within the aimag. MLE calculation for Tuv is N/A because detection rate was 100%.

Summary results for the Anaplasma spp., confirmed through PCR and DNA sequencing of the tick pools, are presented in Table 4. Pools were found to have an overall positivity rate of 3.18% for Anaplasma spp. (12/377), with only ticks sampled from Arkhangai (33% of pools), Uvs (7% of pools) and Uvurkhangai (2% of pools) testing positive. MLE found an average prevalence of 0.69% (95%CI: 0.39–1.19%), with Arkhangai having the highest MLE of 7.79% (95%CI: 4.04–13.72%) and MIR of 6.67% (95%CI: 2.46–10.87%). In contrast, Uvurkhangai had an MLE of 0.53% (95%CI: 0.09–2.87%) with a MIR of 0.53% (95%CI: 0–1.57%).

Table 4

Maximum likelihood estimates of Anaplasma spp. by region based on confirmatory results including 95% confidence intervals.

Anaplasma spp.
Province	Positive pools	Total number of ticks	MLE
			Point	Low	High
Arkhangai	9/27 (33%)	135	7.8	4.0	13.7
Bayankhongor	0/29 (0%)	143	–	–	–
Bulgan	0/14 (0%)	56	–	–	–
Dornod	0/34 (0%)	159	–	–	–
Dornogovi	0/5 (0%)	18	–	–	–
Dundgovi	0/27 (0%)	126	–	–	–
Govi-Altai	0/24 (0%)	120	–	–	–
Khentii	0/8 (0%)	17	–	–	–
Khovd	0/26 (0%)	130	–	–	–
Khuvsgul	0/51 (0%)	232	–	–	–
Sukhbaatar	0/27 (0%)	133	–	–	–
Tuv	0/4 (0%)	21	–	–	–
Uvs	2/29 (7%)	145	1.4	0.4	4.8
Uvurkhangai	1/42 (2%)	188	0.5	0.1	2.9
Zavhan	0/30 (0%)	150	–	–	–
Total	12/377 (3%)	1773	0.7	0.4	1.2

Table 5 summarizes the results for the Bartonella spp., confirmed through PCR and DNA sequencing of the tick pools. The overall Bartonella spp. pool positivity rate was found to be 0.79% (3/377), with Arkhangai being the only region with positive pools (3/27 pools). Maximum likelihood estimates (MLE) found an overall prevalence of 0.17% (95% CI: 0–0.06–0.50%). Arkhangai had a MLE of 2.33% (95%CI: 0.78–6.37%) and a MIR of 2.22% (95%CI: 0–4.71%). A full list of sequence accession numbers by gene target and microorganism can be found in Table 6.

Table 5

Maximum likelihood estimates of Bartonella spp. by region based on confirmatory results including 95% confidence intervals.

Bartonella spp.
Province	Positive pools	Total number of ticks	MLE
			Point	Low	High
Arkhangai	3/27 (11%)	135	2.3	0.8	6.4
Bayankhongor	0/29 (0%)	143	–	–	–
Bulgan	0/14 (0%)	56	–	–	–
Dornod	0/34 (0%)	159	–	–	–
Dornogovi	0/5 (0%)	18	–	–	–
Dundgovi	0/27 (0%)	126	–	–	–
Govi-Altai	0/24 (0%)	120	–	–	–
Khentii	0/8 (0%)	17	–	–	–
Khovd	0/26 (0%)	130	–	–	–
Khuvsgul	0/51 (0%)	232	–	–	–
Sukhbaatar	0/27 (0%)	133	–	–	–
Tuv	0/4 (0%)	21	–	–	–
Uvs	0/29 (0%)	145	–	–	–
Uvurkhangai	0/42 (0%)	188	–	–	–
Zavhan	0/30 (0%)	150	–	–	–
Total	3/377 (0.8%)	1773	0.2	0.0	0.5

Table 6

Identity and Genbank accession numbers for Anaplasma, Bartonella and Rickettsia spp. from pooled samples Dermacentor spp.

Target Gene	Organism	Location	GenBank#	Identity
gltA	Bartonella schoenbuchii	Arkhangai	OM281134-OM281135	99.21% AJ564635.1
	B. melophagi	Arkhangai	OM281136	100% AY692475.1
	Rickettsia raoultii	Arkhangai, Bayankhongor, Dornogovi, Dornod, Dundgovi, Govi-Altai, Khentii, Khovd, Khuvsgul, Sukhbaatar, Tuv, Uvs, Uvurkhangai, Zavhan	OM28112; OM281162; OM281137-OM281146; OM281148-OM281155; OM281157; OM281160-OM281168; OM281170-OM281171; OM2811173-OM281177; OM281179-OM281185; OM281187-OM281192; KU961538;	100% MT178337.1
	R. raoultii	Dornod, Khovd, Khuvsgul, Uvs	OM281156; OM281158; OM281159; OM281169; OM281178; OM281186	100% OK638145.1
	R. sibirica/R. slovaca	Govi-Altai	OM281147	100% MG811709.1;
ompA	R. raoultii	Bayankhongor, Dornod, Govi-Altai, Khentii, Khovd, Khuvsgul, Sukhbaatar, Uvs, Uvurkhangai	OM281193-OM281217	100% MK726326.1
16S rRNA	Anaplasma ovis A. capra A. centrale A. marginale	Arkhangai	OM320148-OM320155	100% MN266936.1
	A. ovis	Uvs	OM320157	100% MN266936.1
groEL	A. ovis	Arkhangai, Uvs	OM281118-OM281120; OM281122-OM281128	99.69% MT268377.1
	A. ovis, A. centrale, A. marginale	Arkhangai	OM28121	99.39% MT268375.1; 92.05 KY305559.1
	A. ovis	Uvs	OM281232	100% MH292916.1
	A. ovis	Arkhangai	OM281229-OM281231	100% MH292916.1
16S rRNA	Coxiella endosymbiont of Dermacentor marginatus	Arkhangai, Dornod, Dornogovi, Khentii, Khuvsgul, Sukhbaatar, Uvurkhangai	OM333168-OM333184	99.44% MZ047981.1

Pathogen detection by tick source

The pathogen pool positivity rate by tick collection source is detailed in Table 2. Rickettsia spp. was detected in 80.4% (95% CI 69.5, 91.3) of tick pools removed from livestock animals, with Tuv having the highest pool positivity rate (100%) and Khuvsgul having the lowest pool positivity rate (33.3%; Table 3). The Rickettsia spp. infection rate in ticks collected from different sources (animal vs. environment including from rock and bush) was compared using Chi-square test and no significant difference was found (Chi-square = 2.7685, df = 1, value of p = 0.09614). Of note, Anaplasma spp. was only detected in tick pools collected from animals, with a pool positivity rate of 23.5% (95% CI 11.9, 35.2). Arkhangai had the highest level of pool positivity, with 47.4% of tick pools collected from animals having Anaplasma spp. DNA present (Table 1). Bartonella spp. was detected in 1.96% (95% CI -1.8, 5.8) of pools of ticks removed from animals and 0.7% (95% CI –0.3, 1.6) of tick pools collected from the environment. As discussed above, all three pools testing positive for Bartonella spp. came from the Arkhangai aimag, with one pool representing ticks collected from animals (5.26% of animal tick pools from Arkhangai (Table 2), and the other two pools being ticks collected from the environment (25% of environmental tick pools from Arkhangai).

Pathogen species confirmation

DNA sequencing allowed for pathogen species confirmation of pools testing positive for the various bacterial groups. Within Rickettsia spp. positive tick pools, the gltA and ompA sequences were analyzed, with the summarizing maximum likelihood (ML) tree presented in Figure 3. As shown in the ML tree, most pools were identified as having R. raoultii (n = 332), although one environmental tick pool from Govi-Altai had a pathogen identified as R. sibirica/Rickettsia slovaca (100% sequence identity). Anaplasma species were identified by analyses of both the 16S rDNA and groEL, placing Anaplasma-positive pools within the A. capra, A. centrale, A. marginale, and A. ovis group, with all Anaplasma-positive pools eventually being grouped within the A. ovis group (n = 12, Figure 4). Finally, ML gltA gene analysis of the three pools that tested positive for Bartonella spp. identified the species as Bartonella melophagi (Figure 5). Figure 6 summarizes the geographic distribution of identified microbial species as well as the proportion of pools within an aimag that tested positive for Anaplasma, Bartonella, or Rickettsia species. A higher proportion of tick pools tested positive for a pathogen in the eastern and western part of Mongolia, which was largely driven by high detection rates of R. raoultii.

Figure 3

Maximum likelihood (ML) tree was constructed from gltA gene (A) and ompA gene (B) of Rickettsia spp. using T92 + G model with 1,000 bootstrap replicates (> 50% are shown on each node). Sequences of tick samples in this study are shown in red letters.

Figure 4

ML tree was constructed from 16S rRNA (A) and groEL genes using K9 + G (B) and TN93 + G + I (C) models, respectively, with 1,000 bootstrap and value over 50% are indicated on each node.

Figure 5

ML tree constructed from gltA gene of Bartonella spp. using T92 + G model with 1,000 bootstrap replicates (> 50% value are shown on each node). Sequence of tick samples in this study are indicated in red letters.

Figure 6

Distribution of identified microbial species and pool positivity rate by aimag. A map demonstrating the geographic distribution of microbial species identified within the sampled tick pools. Each identified microbial species has a different symbol. The aimags are colored according to the proportion of pools that tested positive for Anaplasma, Bartonella, or Rickettsia species.

Discussion

This study continues previous work describing the microbial diversity found within Dermacentor ticks of Mongolia, applying next-generation sequencing to ticks collected from a wider geographic range. Findings from this study reiterate that Rickettsia spp., specifically R. raoultii, are highly prevalent across Mongolia, with aimag pools having a positivity rate of above 50%. While previous work has documented high pool positivity levels for Rickettsia spp. in Dermacentor ticks from southern and central aimags of Mongolia (Fischer et al., 2020; von Fricken et al., 2020b), this study has a much wider geographic range and represents, to our knowledge, the first time NGS methods have been applied to testing Dermacentor ticks in Mongolia. The finding of higher MLEs and pool-positivity rates in the eastern part of Mongolia may indicate a potential hotspot for Dermacentor tick-related Rickettsia spp. exposure, warranting future human and animal serological studies in these areas. Of note, the total MLE in this study (37.30%; 95% CI: 33.50–41.01) is similar to a previous work that sampled ticks from five southern aimags, where the MLE for Dermacentor spp. was 33.2% (95% CI: 30.1–36.2; von Fricken et al., 2020b). Additionally, zero larvae and few nymphs were found across 191 geographically distinct collection events spread out through 15 aimags, which we believe is suggestive of Dermacentor ticks spending earlier life cycle stages underground in rodent burrows, given harsh dry winter seasons common in Mongolia. This theory is also supported by high levels of Rickettsia detection in rodent reservoirs, where 17/18 Meriones unguiculatus (Mongolian Gerbil) tested positive for Rickettsia DNA, with this rodent commonly found across Mongolia (Pulscher et al., 2018). When paired with evidence of transovarial transmission of Rickettsia in Dermacentor ticks from Mongolia (Moore et al., 2018), it is not surprising to observe such high detection rates across this wide geographic range. The Rickettsia species identified through NGS analysis include R. raoultii and R. sibirica/R. slovaca, which aligns with previous reports of R. raoultii and R. sibirica from Dermacentor spp. collected from the Omnogovi, Dornogovi, Govi-Altai, Khovd, Khentii and Bayankhongor aimags (Fischer et al., 2020; von Fricken et al., 2020b). Infections with R. raoultii typically manifest with eschars and lymphadenopathy, although severe cases have been reported with pulmonary edema (Li et al., 2018). Similarly, R. sibirica subspecies present with non-specific flu-like symptoms accompanied by a rash and eschar, with more severe complications such as disseminated intravascular coagulation, renal failure, and neurological symptoms (Nouchi et al., 2018). Therefore, the detection of these pathogens across a wide geographic distribution of ticks within Mongolia represent a major public health threat that is likely under reported in pastoral communities, due to limited access to healthcare in rural regions and low treatment-seeking behaviors within this population (Lkhagvatseren et al., 2019).

While only seen in three aimags (Arkhangai, Uvs, and Uvurkhangai), Anaplasma spp. was still detected in 3.18% of pools overall, eventually being identified as A. ovis. This Anaplasma species causes anaplasmosis in sheep, goats, and wildlife ruminants, often characterized as a subclinical disease which can lead to reduced milk production and spontaneous abortions (Cabezas-Cruz et al., 2019). Of note, all positive pools for Anaplasma spp. came from ticks removed from livestock, primarily sheep and goat, with zero detections occurring in environmental samples. A previous study of Anaplasma spp. and Ehrlichia spp. within ticks collected from central Mongolia reported a similar pattern of high MLE rates when ticks collected from animals were considered separately from ticks collected from the environment (Shao et al., 2020; von Fricken et al., 2020a). This pattern of pathogen distribution within ticks may result from ticks taking partial blood meals from infected livestock hosts and then detaching and reattaching to other livestock hosts, spreading the disease within herds in the process. We do not believe rodents or transovarial transmission plays a significant role in A. ovis transmission cycles given the absence of detection in such a large sample. These findings highlight an important component of Anaplasma spp. disease ecology where livestock act as amplifying hosts. Given these observed patterns of Anaplasma spp. transmission within Mongolia, the lack of Anaplasma spp. detection in many of the aimags in this study should not be interpreted as an actual absence of this microbe group, as many of the pools within these Anaplasma-negative aimags were primarily pools collected from the environment. Further investigation into co-feeding transmission between ticks and potential vectors is warranted, when paired with the high seroprevalence and pathogen detection found in previous studies screening livestock (Zhang et al., 2008; Jiao et al., 2021).

Bartonella spp. are typically transmitted by fleas and lice, however there is an ongoing larger discussion about what potential role ticks play within transmission cycles (Angelakis et al., 2010; Cheslock and Embers, 2019). Studies that have suggested the possible role of ticks as vectors of Bartonella spp. include the reported presence of Bartonella spp. DNA within various tick groups collected from around the world, including Dermacentor ticks, and epidemiological studies have often noted tick bites preceding Bartonella spp.-related illnesses (Wikswo et al., 2007; Angelakis et al., 2010; Zając et al., 2015). Experimental data supporting tick-mediated transmission of Bartonella spp. is limited, but includes the ability of Bartonella spp. to replicate within multiple tick species cell lines (Billeter et al., 2009), the ability of Bartonella-infected ticks to transmit Bartonella infection between animal models while feeding (Noguchi, 1926; Reis et al., 2011), and the detection of Bartonella spp. DNA or bacilli in infected tick midgut, salivary glands, and feces (Cotté et al., 2008; Reis et al., 2011; Wechtaisong et al., 2021). Despite these findings, there is still a lack of consensus regarding the ability of ticks to vector Bartonella species. Regardless of whether ticks play a role in transmitting Bartonella in Mongolia, here we provide further evidence that Bartonella melophagi is present in Mongolia. Within this study, Bartonella melophagi was found in 11.1% of tick pools from the Arkhangai aimag, including a tick pool collected from sheep. Sheep are considered the reservoir host for B. melophagi, with sheep keds being the common insect vector of this pathogen (Maggi et al., 2009). Human infections with this microbe have been reported, resulting in flu-like symptoms, bite-site lesion, neurological symptoms, and heart irregularities (Maggi et al., 2009). Although our group has previously reported B. melophagi within Mongolian sheep, to our knowledge, this is the first report of B. melophagi within ticks collected from Mongolia (Chaorattanakawee et al., 2022). Of note, this Bartonella species has been reported in Dermacentor, Hyalomma, and Rhipicephalus ticks collected from Xinjiang, China (Ni et al., 2021). The detection of B. melophagi in both domestic animals and ticks of Mongolia emphasizes a need for future studies to characterize the disease ecology of this pathogen, including determining the role of ticks in disease transmission and the possibility of transboundary disease movement between Mongolia and China.

All tick pools analyzed by qPCR and conventional PCR targeting the Coxiella burnetii transposase gene were negative, suggesting the absence of this pathogen from the ticks sampled. Of note, 35.8% of pools tested positive for a Coxiella-like bacteria, warranting further investigation.

This study illustrates the utility of NGS in characterizing the diversity of tick-borne pathogens found in Dermacentor ticks collected from geographically distinct locations. We applied NGS to receive a “snapshot” of the various bacterial groups present within tick pools, which then guided confirmatory assays to allow for accurate identification of tick-borne pathogen species. Given the large number of tick-borne pathogens present within Mongolia, the use of a nontargeting analytical method is appropriate to avoid unintentionally excluding the detection of certain microbial species which might be missed relying on other detection processes. Importantly, certain results of this study corroborate what has previously been reported concerning the epidemiology of tick-borne pathogens within Mongolia. This includes the high prevalence of Rickettsia spp., particularly R. raoultii, among Dermacentor ticks collected across the country (Fischer et al., 2020; von Fricken et al., 2020b). We also report the detection of A. ovis exclusively in ticks removed from livestock, which is in agreement with a previous study demonstrating higher levels of Anaplasma spp. in ticks collected from animals (Shao et al., 2020; von Fricken et al., 2020a). While this study expands on the knowledge concerning the geographical distribution of Rickettsia spp. and Anaplasma spp. within Dermacentor ticks, the observation that many of the results are in agreement with previous studies indicates that NGS offers a valid, novel approach for the characterization of tick-borne pathogens. Of note, the novel detection of Bartonella spp. DNA within Dermacentor spp. ticks collected within Mongolia also demonstrates the ability of NGS to discover new pathogen-vector relationships, which may be more difficult to detect using other molecular processes. Future use of NGS to describe the microbial diversity found in the various tick species from Mongolia will further contribute to a more complete characterization of the tick-borne pathogens in circulation within the country. The results from this study will contribute to more detailed risk mapping for tick-borne pathogens, which will help inform disease prevention interventions that benefit populations at increased risk of disease exposure. The abilities of NGS to identify novel vector-pathogen associations will also prove to be vital for local health care practitioners by informing them on the various tick-borne diseases they should include as part of differential diagnoses strategies.

The use of NGS in epidemiological surveys of vector species such as ticks has many advantages over pathogen detection methods that are typically used in such studies, such as PCR and immunofluorescence. While this study focused primarily on clinically relevant pathogens, sequencing data resulting from NGS allows for the creation of a library of microbial sequences, which can promote the tracking of microbial evolution overtime (Deurenberg et al., 2017). Importantly, by surpassing the need for selecting pathogen-specific molecular probes, the use of NGS streamlines the ability of groups to rapidly identify uncommon or previously uncharacterized microbial agents, which may represent emergent diseases (Wu et al., 2021).

Limitations

In this study, ticks were only identified morphologically, which limits our ability to infer findings beyond the genus level. The reliance on morphological identification of ticks may have led to misidentification of the tick species analyzed in this study, hence the decision to keep most of our discussion at the genus level. The decision to pool ticks within this study also introduces some limitations, including difficulties in determining the true prevalence of the various microbial agents that were detected. For example, in instances in which 100% of tick pools are positive, it is not possible to calculate a maximum likelihood estimate, which is what occurred for the Tuv aimag pools for our Rickettsia results. Although there was one pool in which all three pathogens were detected and 12 pools in which two pathogens were detected, discussion of co-infection status of ticks is also complicated by pooling of ticks. While the use of NGS may prove useful for the characterization of pathogenic microbes within ticks and other vector insect species, it is important to note that the detection of microbial DNA does not necessarily indicate the presence of viable microbial organisms within the tick sample. The detection of microbial DNA may also represent remnant DNA from a recent bloodmeal (Tokarz et al., 2019). Therefore, caution must be taken when interpreting NGS results from blood-feeding arthropods.

Conclusion

Here we report the use of NGS to assess the diversity of pathogens within Dermacentor ticks collected from 15 different aimags of Mongolia. The results of this study highlight a high level of Rickettsia detected across all sampled aimags, including the presence of R. sibirica/R. slovaca in Govi-Altai, as well as detections of A. ovis in samples removed from livestock. These findings also highlight the first reported detection of B. melophagi in ticks from this region. Future studies should make use of NGS analysis to further characterize the diversity of pathogens found in other medically relevant tick species within Mongolia, as this method allows for the detection of multiple pathogens simultaneously.

Author’s note

The material in this manuscript has been reviewed by the Walter Reed Army Institute of Research. There is no objection to its presentation and/or publication. The opinions or assertions contained herein are the private views of the author, and are not to be construed as official, or as reflecting true views of the Department of the Army, Department of the Navy, the Department of Defense, or the US Government. Multiple authors are military service members and/or federal/contracted employees of the United States government. The work prepared by DoD personnel in this manuscript was prepared as part of their official duties. Title 17 U.S.C. 105 provides that `copyright protection under this title is not available for any work of the United States Government.’ Title 17 U.S.C. 101 defines a United States Government work as work prepared by a military service member or employee of the United States. Government as part of that person’s official duties.

Data availability statement

The data presented in the study are deposited in the NCBI GenBank repository and can be found at: https://www.ncbi.nlm.nih.gov/genbank/, OM281134-OM281136; OM28112; OM281162; OM281137-OM281146; OM281148 -OM281155; OM281157; OM281160-OM281168; OM281170-OM281171; OM2811173 -OM281177; OM281179-OM281185; OM281187 -OM281192; OM281156; OM281158; OM281159; OM281169; OM281178; OM281186; OM281147; OM281193-OM281217; OM320148 -OM320155; OM320157; OM281118-OM281120; OM281122 -OM281128; OM28121; OM281232; OM281229 -OM281231; OM333168-OM333184.

Author contributions

DA: methodology, formal analysis, sample collection, investigation, and writing—original draft. AL and GM: formal analysis, data visualization, and writing—original draft. RT: methodology, formal analysis, investigation, data curation, writing—review and editing, and visualization. JS and NC: methodology, formal analysis, and investigation. BP-S, JF, and EL: conceptualization, writing—review and editing, supervision, project administration, and funding acquisition. BBd: investigation, writing—review and editing, and project administration. SD and JH: sample collection, conceptualization, and funding acquisition. BBr: morphological identification and sample collection. NT: writing—review and editing, supervision, and project administration. MF: conceptualization, formal analysis, investigation, data curation, writing—original draft, visualization, supervision, and funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

This work was supported/funded by work unit number D0016 with funding provided by the Armed Forces Health Surveillance Division (AFHSD) Global Emerging Infections Surveillance (GEIS) Branch, ProMIS ID # PO133_19_AF_N2. Lab analysis was funded by the AFHSD-GEIS under study #P0128_20_AF_14.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.