Questing ticks were collected by dragging a 1-m² cotton cloth over vegetation in four tourist attraction sites in Khoa Yai National Park (14°26'19.5“N 101°22'20.1”E). Sampling was conducted by six people at each site for 1 h. Ticks were collected in one trip in November 2020 and immediately preserved in 90% ethanol before transporting to the laboratory for further processing. All sampling procedures and experimental manipulations were reviewed and approved as part of the animal collection protocol entitled “Surveillance of Tick- and Flea-Borne Diseases of Public Health Importance” (PN# 21-01). The project was also approved by Mahidol University – Institute animal care and use committee (MU-IACUC 2019/3). Research was conducted in compliance with the Animal Welfare Act and other federal statutes and regulations related to animals and experiments involving animals, and adhered to principles outlined in the “Guide for the Care and Use of Laboratory Animals,” NRC Publication, 2011 edition.

Questing ticks (unfed ticks) were morphologically identified using taxonomic keys (36, 37). Each tick was vortexed for 1 min in 3% sodium hypochlorite and then transferred to a new tube and vortexed for 1 min in 70% alcohol followed by three washes in sterile PBS in the same manner. Ticks were air dried for 10 min on Whatman® filter paper (Sigma-Aldrich, Saint Louis, MO, USA) before DNA extraction. Whole ticks in 250 μl of ATL buffer (lysis buffer, component of the extraction kit) 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 of 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 (elution buffer, component of extraction kit) and stored at−20°C until use. Ultrapure DNA/RNA-free distilled water as well as PBS buffer used for tick surface sterilization were also included as an extraction control.

Following DNA extraction, the bacterial-specific 16S rDNA (V3–V4, a 550-bp fragment) was amplified as previously described (15). Negative control PCR reactions were included in every experimental run using Ultrapure DNA/RNA-free distilled water in place of DNA template. PCR reactions were also performed with eluates from mock DNA extractions as well as from PBS buffer used for tick surface sterilization. PCR product was cleaned using AMPure magnetic bead-based purification system (Beckman Coulter, UK) following the manufacturer's instructions. Purified PCR products were eluted and quantified using the Qubit dsDNA HS Assay Kit (Invitrogen Life Technologies, MA).

The library was prepared with dual indices and Illumina sequencing adapters attached to purified PCR products using the Nextera XT Index Kit following the manufacturer's protocol (Illumina Inc., San Diego, CA). For index control reaction, a combination of index primers that were not used with samples was also included with PCR grade water as template. The number of reads recovered from these particular index combinations were used to filter the cross-contaminations between indexed PCR primers and to identify errors in an Illumina sample sheet. Libraries were cleaned using Agencourt AMPure XP beads. The purity of the libraries was checked on the QIAxcel Advanced System (Qiagen) with a QIAxcel DNA High Resolution Cartridge. Purified amplicon libraries were quantified using the Qubit dsDNA HS Assay Kit (Invitrogen). DNA concentration was calculated and normalized to reach 4.0 nM for each library. Five microliters of DNA from each library was pooled for a NGS run. Pooled libraries were denatured and diluted to a final concentration of 8 pM with a 10% PhiX (Illumina) control. Sequencing was performed using the MiSeq Reagent Kit V3 on the Illumina MiSeq System.

The sequence reads generated by the 16S rRNA on MiSeq sequencers were processed on CLC Genomics workbench v 12.0.3 (Qiagen, Aarhus A/S, http://www.clcbio.com). High-throughput sequences were imported into CLC Genomics Workbench according to quality scores of Illumina pipeline 1.8. In order to achieve the highest quality sequences for clustering, paired reads were merged in CLC microbial genomics module v 4.8 using default settings (mismatch cost = 1; minimum score = 40; gap cost = 4 and maximum unaligned end mismatch = 5). Primer sequences were trimmed from merged reads using parameters (trim using quality scores = 0.01, trim ambiguous nucleotides = 2, and discard read length shorter than 150 bp). Samples were removed from analysis if the number of reads was <100 or <25% from the median (the median number of reads across all samples). Chimeric sequences were detected and removed. Only filtered and merged sequences were clustered into operational taxonomic units (OTUs) according to a threshold of 97% sequence identity. All such processes were performed using CLC microbial genomics module v 4.8. Reference OTU data used in the present study were downloaded from the Greengenes database V13.8 (38). OTUs with combined abundance <10 reads were removed from the downstream analysis. Alpha diversity estimates (Observed OTUs, Simpson's index, and Shannon entropy) were analyzed on the quality-filtered OTU table at the genus level using CLC Microbial Genomics Module v 4.8. Beta diversity analysis for microbiome compositional difference between groups was calculated using a distance-based non-parametric test, the generalized UniFrac distances (39). The statistical significance of microbiome compositional difference between groups (tick genera, tick developmental stages, Francisella persica infection, Coxiella-like endosymbiont infection, and Rickettsia spp. infection statuses) was then compared using PERMANOVA using 99,999 permutations of the distance values. All analyses mentioned here were performed with CLC Microbial Genomics Module v 4.8.

Real-time PCR and PCR assays were performed on 433 individual ticks for the detection of bacteria (Rickettsia spp., Borrelia spp., Anaplasma spp., Ehrlichia spp., and Coxiella-like endosymbiont) and the taxonomic species assignment. Other potential pathogenic bacteria (Francisella spp. and Coxiella spp.) detected by NGS analysis (read count > 1) were also confirmed by real-time PCR and PCR assays. Detailed methods for assays and target gene(s) for selected pathogens are provided as online Supplementary data (Supplementary Table 1). For all real-time PCR, the reaction consisted of 1X Platinum quantitative PCR SuperMix-UDG (Invitrogen) using standard real-time PCR conditions with primer/probe concentrations and annealing temperatures as indicated in Supplementary Table 1. For conventional PCR, the assay was carried out in a 50-μl reaction volume containing 0.5 U of iProof High-Fidelity DNA Polymerase, 200 μM dNTPs, MgCl₂, and primer concentration as indicated (Supplementary Table 1). The PCR conditions consisted of 98°C for 3 min, followed by 40 cycles of 98°C for 10 s, annealing temperature for 30 s (indicated in Supplementary Table 1 for each pathogen), and 72°C for 45 s.

PCR amplicons were purified using the QIAquick® PCR Purification Kit (Qiagen) according to the manufacturer's instructions. The PCR products were cycle-sequenced using an ABI BigDye™ Terminator v3.1 Cycle Sequencing Kit, ethanol precipitated, and run on a SeqStudio Genetic Analyzer (Applied Biosystems ThermoFisher, Thailand). Sequences of each sample and pathogen were assembled using Sequencher™ ver. 5.4.6 (GeneCodes Corp., Ann Arbor, MI). The pathogen sequences were aligned with reference sequences retrieved from the GenBank database using the MUSCLE codon alignment algorithm (40). A maximum likelihood phylogenetic tree was then constructed from bacterial target genes (Supplementary Table 1) using the best fit model of nucleotide substitution with bootstrapping (1,000 replicates) in MEGA 6 (41).

Differences in alpha diversity indices of the bacterial community composition, based on metagenomics data, were determined by Mann–Whitney U test (between two groups) or Kruskal–Wallis test (across all groups) and the critical range (p < 0.05) was determined. Statistical analyses (two-way ANOVA tests and mean, 95% confidence interval) and scatter plots were performed using GraphPad Prism version 5.04 for Windows (GraphPad Software, San Diego, CA. www.graphpad.com). Some graphical illustrations presented in this study as well as the Chi-square independence test and Fisher's exact test were performed in the R environment for statistical computing (42, 43). A nucleotide distance matrix was generated using “Compute Pairwise Distance” in MEGA 6 (41).

A total of 21,229 questing ticks were collected in November 2020 in Khao Yai National Park. The majority of ticks collected were larvae (n = 20,916) accounting for 98.5%, followed by females (n = 145, 0.7%), males (n = 85, 0.4%), and nymphs (n = 83, 0.4%). Species identification was done for the adult stage only and Haemaphysalis lagrangei, Haemaphysalis obesa, and Haemaphysalis longicornis were found at the highest rate (22%−24%), followed by Haemaphysalis shimoga (15%), Dermacentor auratus (9%), Amblyomma testudinarium (4%), Haemaphysalis papuana (1%), and Dermacentor steini (1%). All adult ticks (n = 230) and nymphs (n = 83) and 0.6% of larvae (n = 120) were selected for pathogen detection individually (n = 433). Identification of larval stage was performed on 120 selected samples and Amblyomma spp. (n = 35), Dermacentor spp. (n = 35), and Haemaphysalis spp. (n = 50) were found. In total, there were Haemaphysalis spp. (331, 76.4%), Dermacentor spp. (58, 13.4%), and Amblyomma spp. (44, 10.2%) included in this study. Additionally, a subset of ticks (n = 82) were selected for studying the bacterial community profile using 16S rRNA Next-Generation Sequencing (Table 1). Five to six ticks per species and life stage were selected for NGS, with the exception of immature Haemaphysalis spp., where 10 larvae and 20 nymphs were included. The number of ticks selected for NGS and the final number of ticks that passed the quality filter are described in detail in the online Supplementary data (Supplementary Table 2).

Overall, pathogenic bacteria were detected in 62 out of the total number of examined ticks collected from Khao Yai National Park, for a 14.3% prevalence rate. Pathogens were found in 48 of the tested Haemaphysalis spp. (14.5%), nine in Amblyomma spp. (20.5%), and five in Dermacentor spp. (8.6%). Rickettsia spp. was the most common pathogenic bacteria (7.9%), followed by Ehrlichia spp. (3.2%), and Anaplasma spp. and Borrelia spp. with a similar rate of 1.6% (Table 2). Haemaphysalis spp. (n = 331) were infected by all pathogens mentioned with Rickettsia spp. at the highest rate (7.6%), followed by Ehrlichia spp. (3.2%), Anaplasma spp. (1.8%), and Borrelia spp. (0.9%). In Amblyomma spp. (n = 44) and Dermacentor spp. (n = 58), two bacterial pathogens were found in each genus. Both were infected with Rickettsia (18.2% and 1.7% prevalence rates, respectively). Anaplasma spp. was detected only in Amblyomma spp. (2.3%), while Borrelia spp. was detected only in Dermacentor spp. (6.9%). When examining pathogen prevalence in tick life stages, the infection rate was found to increase from 10.0 to 19.3% from larvae to adults, respectively (Figure 1). The trend of increasing pathogen prevalence by life stage is clearly observed for Rickettsia spp. and Ehrlichia spp.

Sequence and phylogenetic analyses showed that among Rickettsia species detected in ticks, R. montana (n = 21) was the predominant species found in Haemaphysalis spp., while R. raoultii (n = 6) was detected mostly in Amblyomma spp. (n = 5), and three other Rickettsia species (2 = R. heilongjiangensis, 3 = R. monacensis, and 2 = Rickettsia sp.) were detected in Haemaphysalis spp. and Amblyomma spp. (Table 3). Borrelia theileri (n = 7) was the only Borrelia species found, and was found in nearly equal numbers in Dermacentor spp. and Haemaphysalis spp. ticks. The majority of Anaplasmataceae bacteria (7 = Anaplasma spp. and 14 = Ehrlichia spp.) were mostly detected in Haemaphysalis spp. with E. ewingii (n = 9) being the dominant species found. Phylogenetic analyses for all bacteria can be found in online Supplementary Data (Supplementary Figure 1).

Co-infection between pathogenic bacteria examined in this study only occurred in a small number of ticks, primarily female Haemaphysalis spp. (n = 3). All co-infections occurred with Rickettsia spp. and Anaplasmataceae (Ehrlichia spp. or Anaplasma spp.) (Figures 2A,B) and accounted for 4.8% of infected ticks and 0.7% of all ticks examined. All three co-infected ticks were also positive for the Coxiella-like endosymbiont.

The Coxiella-like endosymbiont was detected in 65.8% of all ticks examined, with the highest infection rate in nymphs (86.7%) and females (83.4%) (Table 4). Of all genera, Haemaphysalis spp. exhibited the highest rate of infection with the Coxiella-like endosymbiont, with an 81.0% infection rate. The majority of Haemaphysalis spp. infected were females (95.1%), followed by nymphs (88.9%), males (87.0%), and, to a lesser extent, larvae (24.0%). The second highest infection rate belonged to Amblyomma spp. with a 38.6% infection rate and only detected in females (66.7%) and larvae (37.1%). Surprisingly, no Dermacentor spp. were positive for the Coxiella-like endosymbiont.

Figures 2A,B are scatter plots of individual ticks depicting the pathogenic bacterial infection status compared between two tick populations—those harboring the Coxiella-like endosymbiont and those without. There was no significant difference in the number of ticks infected with Borrelia spp. or Anaplasma spp. between the Coxiella-like endosymbiont-positive and -negative groups. However, ticks harboring the Coxiella-like endosymbiont had greater rates of Ehrlichia spp. or Rickettsia spp. infection than those without the Coxiella-like endosymbiont. Chi-square and Fisher's exact tests were used to test for significant differences of Rickettsia spp. or Ehrlichia spp. infection between endosymbiont-positive and -negative ticks. Results show that the proportion of Rickettsia spp. infection in ticks harboring the Coxiella-like endosymbiont [8.8%, CI (5.5, 12.1%)] was not significantly different from the proportion in ticks without the Coxiella-like endosymbiont [6.1%%, CI (2.2, 9.9%)] with statistical values; Chi-square = 0.638, df = 1, p-value = 0.424. Likewise, the proportion of all pathogenic bacterial infection between ticks harboring the Coxiella-like endosymbiont [16.1%, CI (11.9, 20.4%)] and ticks without Coxiella-like endosymbiont [10.8%, CI (5.8, 15.8%)] was not significantly different (Chi-square = 0.7374, df = 1, p-value = 0.3905). Ehrlichia infection was observed only in ticks harboring Coxiella-like endosymbiont with the proportion of 4.9%, CI [2.4%, 7.4%]. Fisher's Exact test was used to determine the association between Ehrlichia infection and infection with the Coxiella-like endosymbiont. Results indicate that infection by these two bacteria in ticks is dependent (alpha = 0.05, p-value = 0.003), suggesting that the Ehrlichia infection is significantly associated with ticks harboring the Coxiella-like endosymbiont.

Of 433 samples tested in this study, 82 samples were selected for metagenomics NGS analysis (Table 1). After performing sequence quality filters and removing samples with low reads, only 58 samples (range: 1,085–472,236, mean number of reads ± SD = 37,780 ± 98,225) passed the criteria and were subjected to further OTU clustering and alpha and beta diversity analyses. A total of 2,222,970 reads passed quality filters and 268 OTUs were found across all samples. Comparison of the number of passed-filter reads being used for OTU clustering between tick stages showed that larvae (n = 19; 96,713 ± 157,880) generated more reads than nymphs (n = 17; 10,161 ± 9,845) and adults (n = 22; 8,226 ± 8,603). Eight Amblyomma spp., 10 Dermacentor spp., and 40 Haemaphysalis spp. were included in metagenomics analysis (Table 1).

The classification of OTUs from each sample was made against the Greengenes reference database and the similarity threshold was set at 0.97 in the CLC microbial genomics module. There were nine recorded phyla found among all tick samples studied: Proteobacteria (91%), Bacteroidetes (4%), Firmicutes (3%), Cyanobacteria (1%), and Actinobacteria (1%) (Figure 3A). There were five major phyla found in controls for DNA extraction, PCR, and Indexing; however, the majority of Proteobacteria phylum (73%) detected was genus Enterobacteria (89%), and very small amount of genus Rickettsia (0.0092%) (Supplementary Figure 2). Only 1% of reads could not be classified using Greengenes reference and were removed before performing downstream process of abundance analysis. Figure 3B shows the abundance of bacterial taxa in tick genera at the family level. There were one or two predominant bacterial taxa in each tick genus. Rickettsiaceae (83%) was the predominant bacterial taxon in Amblyomma spp., while Francisellaceae (50%) and Methylobacteriaceae (26%) were the major taxa found in Dermacentor spp. Haemaphysalis spp. harbored a more diverse bacterial spectrum than the other two tick genera. Coxiellaceae (43%) and Methylobacteriaceae (24%) were the predominant bacterial taxa found, with a lesser abundance of Rickettsiaceae (4.6%). Beta analysis using distance-based non-parametric test indicated that microbiome composition differed significantly across tick genera by the generalized UniFrac distance (df = 2, pseudo-F = 5.812, p-value = 0.00001).

Alpha diversity estimates of bacterial communities in ticks and statistical significance (p-value) for group comparison are summarized in Table 5. The alpha diversity analyses showed that the bacterial community of Haemaphysalis spp. (n = 40) was significantly different from Amblyomma spp. (n = 8) with p-value < 0.05 for all tests (Observed OTUs, Simpson's index, Shannon entropy) when measured at the genus level as shown in Figure 4A. While there were few ticks positive for Francisella persica (n = 4), alpha diversity estimates indicated that their bacterial profiles were quite different from ticks that were not infected with F. persica with statistical significance (p-value < 0.001) (Figure 4B).

As three Haemaphysalis spp. life stages (Larva = 9, Nymph = 15, Male = 6, Female = 10) were included in NGS, this genus was selected for testing the bacterial community difference among tick life stages (Figures 5A,B). Alpha diversity estimates (Simpson's index and Shannon entropy) showed that while the overall bacterial profiles among life stages were not significantly different with regard to bacterial taxa present, the number of observed OTUs among life stages was significantly different (Table 5; Figure 5B). Likewise, the beta diversity analysis using distance-based non-parametric test showed no significant difference across Haemaphysalis spp. life stages by the generalized UniFrac distance (df = 3, pseudo-F = 0.905, p-value = 0.5638). Comparing Coxiellaceae abundance across life stages, this group comprised 57–59% abundance of the total bacteria community in male and female ticks, 40% in nymphs, and 28% in larvae (Figure 5A).