A novel Rickettsia, Candidatus Rickettsia takensis, and the first record of Candidatus Rickettsia laoensis in Dermacentor from Northwestern Thailand

Three hundred and forty-four tick samples were collected from vegetation at Taksin Maharat National Park, Tak province, northwestern Thailand. They were morphologically identified and molecularly confirmed by 16S rRNA and COI genes as Dermacentor laothaiensis (n = 105), D. steini (n = 139), and D. auratus (n = 100). These ticks were examined for the spotted fever group rickettsiae (SFGRs) using PCR and DNA sequencing of six genes; 17-kDa, gltA, 16S rRNA, ompA, ompB, and sca4. Of these ticks, 6.10% (21/344) gave positive results for the presence of SFGRs. Phylogenetic analyses of the SFGRs clearly indicated that a novel genotype assigned as Candidatus Rickettsia takensis was detected in D. laothaiensis (19/105) and at lesser frequency in D. steini (1/139). Furthermore, Candidatus Rickettsia laoensis was also found at a low frequency in D. auratus (1/100), the first record in Thailand. Although, the pathogenicities of these SFGRs remain unknown, our findings suggest potential risks of SFGRs being transmitted via ticks near the border between Thailand and Myanmar, a gateway of daily migrations of local people and visitors both legal and illegal.


Results
Tick species identification. A total of 344 tick samples from the study area comprised three morphological species of Dermacentor ticks. These were D. laothaiensis (n = 105: 36 males, 69 females), D. steini (n = 139: 68 males, 71 females), and D. auratus (n = 100: 42 males, 58 females) (Fig. 1a). The morphology of D. laothaiensis was observed to be very similar to D. steini, and both were clearly different from D. auratus (Fig. 1b-j). D. steini showed closely spaced external and internal spurs on coxa-I, a distinct narrow dark border stripe on the densely punctuated scutum or pseudoscutum, and a narrowly to broadly V-shaped aperture of the female genital structure. D. laothaiensis, otherwise similar in morphology to D. steini, differed in having an indistinct broad hazy central brown stripe on the scutum or pseudoscutum, moderately dense large punctations, and very narrowly V-shaped genital aperture. D. auratus exhibited widely spaced external and internal spurs on coxa-I, a distinct narrow line through the center of scutum/ pseudoscutum, and a U-shaped genital aperture; Fig. 1b-j. Phylogenetic analysis results based on 16S rRNA and COI genes confirmed the existence of these 3 related species (Figs. 2a, c and S1). The pairwise diagram for these tick species based on 16S rRNA gene (Fig. 2b) demonstrated that D. laothaiensis had 91%, 91.8%, and 92% similarity compared with D. steini from Thailand (ON705029), D. auratus from Thailand (ON705025), and D. steini from Malaysia (MK296405), respectively. Their corresponding total numbers of differentiation positions of nucleotide base sequences alignments were 36 (27 bases/ 9 gaps), 33 (31 bases/ 2 gaps), and 32 (29 bases/ 3 gaps), respectively. Thailand D. steini was 91% similar to D. auratus from Thailand (ON705025) and 95.8% similar to D. steini from Malaysia (MK296405). It differed from both taxa in positions of nucleotide base sequences alignments by 36 (31 bases/ 5 gaps), and 17 (15 bases/ 2 gaps), respectively. Our D. auratus tick sample, showed a 99.5% match with both a D. auratus specimen originating from Malaysia (MT914184) and one other from Thailand (KC170746). Their total numbers of differentiation positions of nucleotide base sequences alignments were 2 (2 bases/ 0 gap) and 2 (1 base/ 1 gap), respectively. Likewise, the pairwise diagram for the tick species using the COI gene (Fig. 2d)  Infection rate and molecular identification of SFGRs. The prevalence of rickettsiae from tick samples was investigated. The results using three housekeeping genes commonly used for SFGR screening, namely 17-kDa (434 bp), gltA (380 bp), and 16S rRNA (1,328 bp), all showed a positive rate of 6.10% (21/344). Further analyses of the amplified DNA samples employing three genes commonly used in Rickettsia genotyping, i.e., ompA (532 bp), ompB (1,295 bp), and sca4 (1,100 bp), clearly confirmed the presence of SFGRs in these positive ticks. The three species of ticks and their Rickettsia infection rates were as follows: 18.10% (19/105) in D. laothaiensis, 0.72% (1/139) in D. steini, and 1% (1/100) in D. auratus (Fig. 1). All positive amplicons were then sequenced and phylogenetic trees were constructed for the six individual targeted genes (17-kDa, gltA, 16S rRNA, ompA, ompB, and sca4). They showed more or less similar patterns of phylogenetic relationships with the closely related SFGRs from the NCBI database (Fig. 3). Our findings showed that there were two detected SFGRs in two separate clades as detailed below.
Clade 1, the isolated clade, included only our novel genotype, Ca. R. takensis. Based on 17-kDa, the novel genotype showed phylogenetic relationships with R. japonica from China and Japan, Ca. R. vini from the Czech Republic, Rickettsia sp. TCM1 from Thailand and R. fournieri from Australia (Fig. 3a). However, its gltA was closely related to Rickettsia sp. ARRL2016-159 from Australia, R. japonica from China and Japan, and Ca. R. vini from the Czech Republic (Fig. 3c). The 16S rRNA gene, ompA, ompB, and sca4 all demonstrated their closest relationships with R. fournieri from Australia ( Fig. 3e,g,i,k).
Clade 2 included another SFGR detected in the present study, Ca. R. laoensis (DAT10M3). This may be the first Thailand record of Ca. R. laoensis. Based on the 17-kDa gene, the Thailand specimen grouped with Ca. R. laoensis from Laos and Rickettsia sp. from India (Fig. 3a). Its gltA gene grouped with Rickettsia sp. from India, www.nature.com/scientificreports/ Rickettsia sp. RDla440 from the Thailand-Myanmar border and Ca. R. laoensis from Taiwan (Fig. 3c). On 16S rRNA, it grouped with Rickettsia sp. from Pakistan (Fig. 3e). Close phylogenetic relationships with Ca. R. laoensis were notably demonstrated for ompA from India and Laos (Fig. 3g), ompB and sca4 from Laos (Fig. 3i,k). The results of these phylogenetic analyses supports the conclusion that there may be two unknown groups of Rickettsia spp. occurring in these tick vectors, which are tentatively referred to Ca. R. takensis and Ca. R. laoensis as mentioned above. However, some other SFGRs also showed phylogenetic relationships close to our discovered SFGRs. Therefore, to validate the species classification of our two putative detected SFGRs, their sequence alignments, and pairwise comparisons were performed against reference SFGRs. The results of the nucleotide similarity percentages determined were then applied to pose criteria for a novel Rickettsia determination using gene sequences 20  The highest percentage, indicating the closest similarity in each target gene is shown in bold font in the diagram. The 17-kDa gene which was used to screen Rickettsia in the population, showed the highest value (99.5%; 2 bases/ 0 gap/ 1 amino acid) for several SFGRs (R. fournieri from Australia, Ca. R. vini from the Czech Republic, and Rickettsia sp. TCM1 from Thailand; Fig. 3b). The gltA gene showed the highest value of similarity (99.7%, 1 base/ 0 gap/ 0 amino acid) for R. japonica from Thailand (PMK94), R. heilongjiangensis from Japan, and Rickettsia sp. ARRL2016-159 from Australia (Fig. 3d). For 16S rRNA the highest values of 99.4% (6 bases/ 1 gap) were obtained for R. fournieri from Australia, R. japonica from China (LA16/2015) and Japan (YHM and YM) (Fig. 3f). The ompA, ompB, and sca4 genes demonstrated similar results, finding that our Ca. R. takensis was closest to the R. fournieri from Australia. Their highest percentages and detail of variables were as follows: ompA 96.4% (12 bases/ 3 gaps/ 13 amino acids) (Fig. 3h), ompB 98.1% (15 bases/ 6 gaps/ 15 amino acids) (Fig. 3j), and sca4 98.9% (11 bases/ 0 gap/ 7 amino acids) (Fig. 3l). In addition, the pairwise comparisons between Ca. R. takensis and either Ca. R. laoensis (DAT10M3) from Thailand or Rickettsia sp. RDla440 from the Thailand-Myanmar border, showed exact similarity of 99.3% (2 bases/ 0 gap/ 1 amino acid) based on the only data available for the gltA gene (Fig. 3d). Based on the criteria for a novel Rickettsia determination using gene sequences of 16S rRNA, gltA, ompA, ompB, and sca4, an isolate can be classified as a new Rickettsia genotype when it exhibits no more than one of the following degrees of nucleotide similarity with the most homologous validated species (cutoff percentages ≥ 99.8, ≥ 99.9, ≥ 98.8, ≥ 99.2, and ≥ 99.3) 20 . The result of Ca. R. takensis fulfilled this criterion with the degrees of similarity 99.4, 99.7, 96.4, 98.1, and 98.9, respectively. None of the degrees were equal to or over the cutoff percentages of the standard criteria. This indicated the occurrence of a new Rickettsia genotype. Moreover, this Rickettsia was detected in ticks, not from rickettsial isolates: therefore Candidatus status could be applied. Therefore, this evidence showed that we found a novel genotype of SFGR in D. laothaiensis and D. steini tick samples. In addition, this novel Rickettsia was assigned as Candidatus Rickettsia takensis after the name of the collecting site, Tak province.
According to the evidence, we also found Ca. R. laoensis (DAT10M3) infection in one male in the D. auratus tick samples. The highest percentages of similarities based on the targeted genes were shown in pairwise comparisons. (1) DAT10M3 and Ca. R. laoensis from Laos ompA 100%, ompB 100%, and sca4 99.9% (1 base/ 0 gap/ 0 amino acid) (Fig. 3h, 3j, and 3l).  (Fig. 3f). Overall, phylogenetic analyses and sequence alignments indicated that DAT10M3 and Ca. R. laoensis from Laos shared almost identical genotypes. This shows that Ca. R. laoensis is present on the northwestern border of Thailand. However, there have not yet been any reports on whether Ca. R. laoensis is a pathogen although it was detected in ticks (Haemaphysalis sp. 16 and D. auratus in the present study).

Discussion
More than 40 species of Dermacentor ticks are known worldwide [23][24][25] . So far, eight species (D. auratus, D. compactus, D. steini, D. filippovae, D. pasteuri, D. laothaiensis, D. falsosteini, and D. tricuspis) have been found in Thailand [24][25][26][27][28] . Of these, three species were morphologically and molecularly identified as D. laothaiensis, D. steini, and D. auratus in the present study. D. laothaiensis was recently described, and it was reported as morphologically closely related to D. steini 26 . However, we found that D. laothaiensis was quite distinct from D. steini as these two species resolving into separate clades during phylogenetic analyses. The distinction was further supported as data was harvested from nucleotide sequence alignments (Fig. S2) demonstrating differentiation during pairwise comparisons among the tick species (Fig. 2). The results of our molecular analyses and phylogenetic relationships based on 16S rRNA and COI genes constitute the first molecular evidence indicating the difference between the D. laothaiensis and D. steini, and also showed that both were distinct from D. auratus.
In Thailand, several studies have examined the detection of Rickettsia in Dermacentor ticks collected from vegetation and mammal hosts. Malaisri Tak (16°N 99°E), another border province with the same longitude comparing with Kanchanaburi province as mentioned above that there were some reports of SFGRs. This may indicate the distribution of SFGRs along the Thailand-Myanmar border.
In this study, our phylogenetic classification of Ca. R. takensis lacked clarity because its genotype showed close relationships to (1) R. fournieri from Australia and R. japonica from China and Japan (16S rRNA) (2) Rickettsia sp. from Australia, R. japonica from Thailand, and R. heilongjiangensis from Japan (gltA) and (3) R. fournieri from Australia (ompA, ompB, and sca4) (Fig. 3d,f,h,j,l). However, individually targeted gene trees demonstrated a unique Rickettsia genotype resolving separately from the other reference clades (Fig. 3a,c,e,g,i,k). This was likely supported by the results summarized in the diagrams of Fig. 3. The criteria for classifying Rickettsia species 20 were applied, and it was confirmed that Ca. R. takensis was a member of the genus Rickettsia according to its 16S rRNA and gltA homology against some other validated Rickettsia spp. by the degree of homology of ≥ 98.1% and ≥ 86.5%, respectively for these two genes. Furthermore, the presence of the ompA gene also confirmed that it belonged to the SFGRs. Again, when we examined Rickettsia species identification by MLST, we found that in no case was the degree of nucleotide similarity (99.4, 99.7, 96.4, 98.1, and 98.9), compared to the most homologous species from the NCBI database ( Fig. 3 and Table 1), equal to or over the respective standard criteria cutoff percentages (≥ 99.8, ≥ 99.9, ≥ 98.8, ≥ 99.2, and ≥ 99.3). Therefore, these results indicated the uniqueness of a novel genotype of Ca. R. takensis reported in this study. Moreover, current investigations of rickettsiae associated with Dermacentor ticks in the vicinity of the Thailand-Myanmar border showed the presence of a novel genotype, Ca. R. takensis, predominantly detected in D. laothaiensis and to a lesser extent in D. steini. Ca. R. takensis found in all 20 ticks exhibited similar nucleotide sequences for the six genes (17-kDa, gltA, 16S rRNA, ompA, ompB, and sca4) examined. Surprisingly, only one male of D. auratus was infected with Ca. R. laoensis (DAT10M3). Pairwise comparisons of nucleotide sequences of these two genotypes of rickettsiae based on the targeted genes showed 97.4%, 99.3%, 98.6%, 89.9%, 95.2%, and 97.5% similarities, respectively. These data confirmed that our detected SFGRs, either Ca. R. takensis or Ca. R. laoensis, were two distinct Rickettsia species that infected the population of our tick samples.
Another differentiated SFGRs genotype, Ca. R. laoensis (DAT10M3), was classified and this yielded the first molecular evidence for its presence in Thailand. Importantly, based on the gltA gene, our sample DAT10M3 exhibited a close relationship with the Rickettsia sp. RDla440, the SFGRs reported in Dermacentor larvae from Kanchanaburi 12 . However, the gltA gene sequence analysis was the only available data for RDla440, although it yielded 100% similarity compared to our DAT10M3 (299 bp) sequence alignment. Subsequently it was reported that RDla440 appeared to be more closely related to the sequences of Rickettsia sp. strain DnS14 and Rickettsia sp. strain RpA4, differing in only 2 bp (99.7% similarity) based on gltA (1,108 bp) 12 . However, both DnS14 and RpA4 were later identified as R. raoultii based on 16S rRNA, gltA, ompA, ompB, and sca4 genes 35 . Therefore, although the results with the gltA gene indicated DAT10M3 had 99.7% similarity (variable of 1 base) with either DnS14 or RpA4 strains detected in ticks from the former Soviet Union 36 , data from a single gene might be insufficient for SFGRs classification. Likewise, based on our results, our sample DAT10M3 was most likely Ca. R. laoensis, related to samples of Ca. R. laoensis reported in Laos 16 . However, it was not certain whether the RDla440 detected earlier at the same border longitude of this was in fact Ca. R. laoensis due to the limited molecular evidence. Although gltA has frequently been used as a target for generic diagnostics based on PCR because this approach can easily identify a number of Rickettsia [37][38][39] , this gene is also highly conserved. Therefore, gltA on its own might not be enough to differentiate between closely related Rickettsia species. Consequently, the MLST of 17-kDa, gltA, 16S rRNA, ompA, ompB, and sca4 were suggested as potentially helpful and appropriate for the reconstruction of the evolutionary relationship of diverse but closely related Rickettsia species, as documented in our results and previous reports 1,2,16,38,40,41 .
Ca. R. laoensis infecting Haemaphysalis nymphs (4.5%) was first reported in Khammouan province, Laos 16 . We found that Ca. R. laoensis in Thailand and Laos showed the highest percentage values during pairwise comparison of the five genes, 17-kDa, gltA, ompA, ompB, and sca4, with 99.5%, 99.6%, 100%, 100%, and 99.9% similarity, respectively. Khammouan province in Laos is located near the northeastern border of Thailand, bordering the Mekong River, about 700 km east of Tak province which lies on the Myanmar border (Fig. 1). This geographical evidence seems to suggest that Ca. R. laoensis is probably a common species, and widely distributed in Thailand, at least in the upper part of the country.  41    www.nature.com/scientificreports/ around Southeast Asia, including Thailand. As well as seeking to detect Ca. R. laoensis in different species of tick vectors in this region, more information about its pathogenicity is needed. Therefore, more fieldwork investigating Ca. R. laoensis infection in ticks and humans in this region is required. Although the information on species diversity of rickettsiae related to the emerging rickettsioses via tick vectors is still poorly understood in Thailand, our study has enhanced knowledge of SFGRs harbored in these 3 species of Dermacentor ticks in this particular area of northwestern Thailand. Thus, the prevalence of these SFGRs infected tick vectors along the border of Thailand and Myanmar, extending for approximately 2,400 km from the north to the south, will be very useful and informative because these areas are gateways for the movement of local people between the two countries. The border areas also involve legal and illegal livestock trading, increasing the potential risk of SFGRs exposure in this region. Furthermore, the polymorphic genotypes of Ca. R. takensis and Ca. R. laoensis discovered in this study indicated the genetic diversity of SFGRs in this area. In addition, understanding the movements of wildlife, especially mammalian hosts of ticks, may contribute to the knowledge of the distribution of rickettsial agents in this particular forested area. The data seem to suggest a potential risk for the distribution of this SFGRs in neighboring countries, including Thailand, Laos, and Myanmar. The diagnosis and management of rickettsial infections in this region remains challenging, for instance, the results of clinical differential diagnoses for SFGRs are usually disregarded in resource-poor settings because of the limitations of immunofluorescence serologic tests caused by antigenic cross-reactions 43 . While the identification of phylogenetically distinct SFGRs to species by molecular methods is well known to have practical applications, public health facilities and attention related to rickettsiae needs to be improved to enable greater accessibility in Southeast Asia.

Materials and methods
Tick collection and species identification. A total of 344 adults of Dermacentor ticks were collected from vegetation in the forest of a tourist trail (with shrubs and trees along both sides) at Taksin Maharat National Park, Mae Sot district, Tak province, about 20 km from the Thailand-Myanmar border (Fig. 1a). The collecting of tick samples was conducted upon approval from the Department of National Parks, Wildlife and Plant Conservation of Thailand (permission no. TS 0907.4/10,680). Collections were made at 754 m asl, and ticks were searched for and collected, with blunt forceps or hand picking from vegetation in 2016 with an average temperature of 23 °C, and an average humidity of 80% for the environmental conditions. All tick samples were kept in 70% ethanol and stored at −20 °C for laboratory study. Each tick was morphologically identified using available descriptions 26,44,45 . The samples of these ticks were subsequently confirmed by molecular techniques using DNA sequencing of two genes; 16S ribosomal RNA (16S rRNA: 16S + 1, 16S-1) and cytochrome c oxidase I (COI: RON, TCOIR), as described by Black and Piesman and Simon et al. 46,47 . The references for primers used are listed in Table 2. DNA extraction. DNA extractions of all tick samples were performed after they had been morphologically identified. An individual tick was serially rinsed and cleaned thoroughly in 70% ethanol, 10% sodium hypochlorite (NaOCl), and sterilized distilled water three times (1 min each) to avoid external bacterial contamination. Rickettsia detection, identification and sequencing analyses. PCR was performed to amplify rickettsial DNA from extracted products of each Dermacentor tick. Overall, six targeted partial fragments of genes were used (Table 2). Initially, detection and screening to determine the bacteria of the genus Rickettsia were carried out employing the genus-common 17-kDa surface antigen gene (17-kDa), the pan bacterial gene encoding 16S rRNA (16S rRNA) and the citrate synthase gene (gltA). Rickettsia-positive rates were then calculated. Additionally, we selected positive samples to test with other PCR amplifications that targeted the SFGRs-specific 190-kDa outer membrane protein A gene (ompA), the 120-kDa outer membrane protein B gene (ompB), and the PS120 or D protein-encoding gene (sca4) to determine their SFGRs features. The PCR amplicons of 17-kDa 434 bp, 16S rRNA 1,328 bp, ompA 532 bp, ompB 1,295 bp, and sca4 1,100 bp were purified using the NucleoSpin Gel and PCR Clean-up Mini kit (MACHEREY-NAGEL, Germany). The purified amplicons were sequenced in both directions at Macrogen Co., LTD (Seoul, Korea). Furthermore, any ambiguous positive PCR products were gel purified, inserted into the pGEM® -T easy vector system (Promega, USA), and transformed into E. coli (strain JM109) using standard cloning protocols (or according to the manufacturer's protocol). A positive cloning colony was grown on an LB-agar plate (Luria-Bertani) containing 20 mg/ml IPTG, 20 mg/ml X-Gal and 100 μg/ml ampicillin. The recombinant plasmids were also extracted and subjected to DNA sequencing.
Data analysis. Data analyses of the obtained nucleotide sequences were primarily performed to determine the highest similarity using the Basic Local Alignment Search Tool (BLAST) 52 . Then, multiple sequence alignments with the related DNA sequences (Table S1 and S2) retrieved from the GenBank database were aligned by the CLUSTAL method, BioEdit v.2.0.0 program, while the amino acid polymorphisms were determined using the Amino Acid Explorer of the NCBI server 52 . Phylogenetic trees were constructed using Neighbor-joining (NJ) method based on 16S rRNA and COI genes for tick species identification and Maximum-parsimony (MP) method based on 17-kDa, gltA, 16S rRNA, ompA, ompB, and sca4 genes for rickettsial genotype identification. These analyses were performed in the Molecular Evolutionary Genetics Analysis (MEGA) 7.0 software 53 . Distance matrix analysis was generated by the Tamura 3-parameters with gamma distribution (T92 + G) and Subtree-Pruning-Regrafting (SPR) for multiple substitutions. Bootstrap values were based on 1,000 replicates to estimate the support for nodes within the phylogenetic trees.

Data availability
The nucleotide sequences of the Rickettsia obtained in this study are deposited to GenBank database (https:// www. ncbi. nlm. nih. gov/ genba nk/).