Sequence Analysis of the TSA 56 kDa and GroEL Gene Reveals the Circulation of a Novel Genotype of Orientia tsutsugamushi in Human Cases of Scrub Typhus in Thanjavur District, Tamil Nadu, Southeastern India

1. Introduction

Scrub typhus (St), a re-merging mite-borne infection of humans, has a widespread distribution in Southeast Asia (SEA). Orientia tsutsugamushi (Ot), the causative agent of St, has been found to exist as multiple genotypes and remain in circulation among rodents and mite vectors in nature. The incriminated primary mite vectors are classified under the genus Leptotrombidium sp, Schoengastiella sp, and Ornithonyssus sp. [1,2,3]. The disease is known to be endemic mainly in the “Tsutsugamushi Triangle” i.e., in the Southeast Asian countries, it is estimated that around one billion people are suffering from scrub typhus fever and, one million people are affected. It is considered one of the leading infectious diseases in the Asia-Pacific region [4,5]. In India, the case fatality ratio of scrub typhus ranged from 12 to 623 during the years 2002-2014 and St cases have been widely distributed in many states in the country particularly, in Himachal Pradesh, Assam, and Tamil Nadu [6,7]. Scrub typhus is misdiagnosed and underreported because of its overlapping signs and symptoms with the other commonly occurring leptospirosis, dengue, brucellosis, and typhoid etc [5,8,9]. Though the causative agent was believed to exist as a single species (Orientia tsutsugamushi), the recent discovery of two new candidatus species (Candidatus O. chuto and Candidatus O. chiloensis) has indicated its rich diversity [10,11]. Globally, more than 40 antigenic subtypes of Ot have been documented [12]. Karp, Kato, Gilliam, Boryong, Ikeda, Kawasaki, Kuroki, Saitama and Shimokoshi are key strains, and TA678, TA686, TA716, TA763, and TH1817 have been considered prototype strains [9].

Though ST is a re-emerging zoonotic problem in India, its rich genetic diversity of the circulating pathogen has not been determined in endemic areas. Even though scrub typhus is considered a seriously emerging public health problem in this area (proportion of cases- 37.6%) [13], the information available about the circulating genotypes especially in the southeastern part of the country is scanty. The Kato, Karp-like, Gilliam, Ikeda, Neimeng-65, TA716, TA763 genotypes have been identified from different states in India [14,15,16,17]. However, routine clinical cases are treated based on the ELISA tests, which include the 56-kD outer membrane recombinant protein of Karp, Kato, Gilliam and TA716 genotypes in the diagnostic panel [18].

O. tsutsugamushi has adapted to survive in the mite population through vertical and horizontal transmission routes in the environment [19]. Since this organism infects and circulates in a simple mite to the rodents and humans, undergo genetic recombination, which could result in the emergence of variants [20]. These variants act as virulent and develop clinical complications in humans [16]. Hence, an in-depth understanding of circulating genotypic study with the appropriate markers is essential. Among the different molecular markers, the 56 kDA and GroEL markers are more accurate in the detection of genotypes [19,21]. The available data on the genetic diversity of the pathogen primarily depend on the analysis of the 56 kDa gene [19,22]. However, the inclusion of additional genes-based analysis could elucidate the in-depth dynamics of the genetic diversity of the pathogen. Hence, in the present attempt, the GroEL and 56 kDa genes were used for genotyping of the Ot in this area.

Scrub typhus fever cases have been routinely diagnosed in patients seeking medical assistance at Government Thanjavur Medical College Hospital in recent years. Due to the lack of a genome-based investigation, no data is available on the circulating genotypes in this area. The GroEL gene based earlier study with a limited number of human and rodent samples in Theni district, South India, revealed the circulation of three genotypic variants. Therefore, the present study aims to explore the circulating Ot genotypes by analyzing the GroEL and 56 kDa TSA protein genes in and around the Thanjavur district (Figure 1).

2. Materials and Methods

2.1. Study Area

Thanjavur districts, which lie on the East Coast of Tamil Nadu between 78° 45’ and 70° 25’ E; 9° 50’ and 11° 25’ N in the Cauvery delta, also called the “rice bowl” of Tamil Nadu. The inhabitants are involved in agriculture activities mainly paddy cultivation [23], the rodent population is abundant and different species of rodents were recorded with mite infestation in domestic and peri-domestic areas adjacent to rice fields [2].

2.2. Sample Collection

The blood sample collected from the patient clinically suspected of St infection, seeking medical assistance from the Government Medical College & Hospital, at Thanjavur has been utilized for the detection of genotype. After serological diagnosis by commercial IgM-ELISA kits (Scrub Typhus Detect™ IgM ELISA, IndBios International, Inc, USA), the positive samples (n=59) were subjected to the genome analysis study. A molecular genotypic assay was carried out at the ICMR-Vector Control Research Centre, Field Unit. The necessary ethical clearance was obtained from the Institutional Human Ethical Committee (TMC/IHEC: 1118/2023).

2.3. Detection of O. tsutsugamushi by Nested PCR

The Genomic DNA was extracted from 59 serum samples using the commercially available QIAamp DNA Blood Mini Kit, (Cat No. 51104, Qiagen, Germany). The DNA extracted were subjected to a nested polymerase chain reaction (nPCR) assay targeting the Ot GroEL gene as per the protocol described by Weihong Li et al. [24]. This gene encodes a 60-kDa heat shock protein that is more diverse and highly conserved among O. tsutsugamushi but sufficiently variable to differentiate the genus Orientia and Rickettsia with 99.5% identity (21). Two µl of DNA was used as a template for the first round of PCR with the primers Gro-1, 5’-AAGAAGGA/CGTGATAAC-3’ and Gro-2, 5’-ACTTCA/CGTAGCACC-3’. Subsequently, one microliter of the first round PCR product was further used as the template for nested PCR with the primers TF1, 5’-ATATATCACAGTACTTTGCAAC-3’ and TR2, and 5’-GTTCCTAACTTAGATGTATCAT-3’ to amplify a product of 365 bp of the GroEL gene fragment in a final volume of 25 µl for both PCR reactions. In the post-PCR processing, the amplified DNA fragment was resolved in a 2% agarose gel by loading 5 µl of the PCR product with the 100 bp DNA ladder and a negative (water) control loaded alongside. The remaining amplified gene fragment was further subjected to Sanger DNA sequencing in both forward and reverse directions.

2.4. PCR Amplification of 56 kDA VDI, VII and VIII

All the fourteen samples amplified with the primer specific for GroEl gene were subjected to nested PCR amplification of the final fragment size of ~900 base pair of the 56 kDa gene, with the following outer and inner primers, JG-OtF584 (5’-CAA TGT CTG CGT TGT CGT TGC) and RTS9 (5’-ACAGAT GCA CTA TTA GGC AA), and the inner primers were F (5’-AGC GCTAGG TTT ATT AGC AT) and RTS8 (5’-AGG ATT AGA GTG TGG TCCTT) respectively as described by Toon Ruang-areerate et al. [22]. Five µl of the PCR product loaded in 1% agarose gel and after the confirmation of bands, the remaining PCR amplified products were purified and custom DNA sequenced.

2.5. DNA Sequencing and Phylogenetic Analysis

The DNA sequencing chromatograms were manually examined and edited by comparing both forward and reverse DNA strands. The final DNA sequences thus obtained were subjected to the nucleotide BLAST sequence analysis tool in the NCBI database to determine the sequence similarity. Fourteen partial GroEL and twelve 56 kDa DNA sequences generated in this study were deposited in the NCBI GenBank database and accession numbers were obtained (GenBank Acc. nos. GroEL: OR887439-OR887452; 56 kDA: PQ381691-PQ381702). The nearest neighbour sequences exhibited with low e-value and higher the nucleotide identity in BLAST analysis against the respective GroEL and 56 kDa study sequences were grouped. In addition, the Ot genotype-specific reference sequences for the two genes belonging to Karp, Kato, Gilliam-types, and TA716 strains were also included to calculate the nucleotide diversity (d) and check the phylogenetic relationship among them. Phylogenetic analysis was conducted using the Molecular Evolutionary Genetics Analysis [MEGA 11] software [25]. The mean nucleotide diversity of the entire GroEL study sequences, between and within the genotypes was calculated after grouping the DNA sequences based on the phylogenetic clades.

Similarity profile and bootscan analysis for checking the similarity of the study sequences with neighboring reference sequences and to detect the event of recombination among the strain was carried out by the Simplot software 3.5.1 [26]. Only three prototype Karp, Kato, and TA716 sequences were included to avoid the noise during the Simplot analysis. A neighbor joining similarity as well as recombinant distance plot was drawn with nucleotide substitution model of Kimura 2 parameter (K2P), window, step size of 200 and 20, respectively. The consensus output was obtained with computing 1000 bootstrap replicate and the similarity and recombination was detected among the study sequences.

3. Results

3.1. GroEL Gene Sequence Analysis

A total of 59 samples were processed by Nested PCR, and the GroEL gene-specific product (361 bp) was amplified in 14 (23.7%) samples. The analysis of the fourteen DNA sequences showed three distinct phylogenetic clades viz, new Ot-TJTN, Karp-like and Kuroki-Gilliam types in the Neighbour-Joining (NJ) dendrogram (Figure 2a). Further, to check the consistency of the phylogenetic inference obtained in the NJ tree, the Maximum Likelihood tree (ML) and Maximum Parsimony tree (MP) were constructed with the same DNA sequences, which showed consistent genotypic clades as obtained in the NJ tree. Among the 14 O. tsutsugamushi sequences generated in this study, 9 (64%) sequences formed a distinct new genotypic cluster related to TA716-like Thai strain, 4 (28.5%), 1 (7%) sequence formed Karp-like and Kuroki-Gilliam types, respectively.

The percentage of nucleotide identity for the 9 sequences was 100% with the DNA sequences submitted from south India, GenBank Acc.nos. ON156002 and CP166954, followed by 97.71% (EF551288) from Thailand, 97.43% (CP044031; LS398552) from China and UK, respectively. The percentage of BLAST identity scores for four other Karp-like sequences, with the ON15600, JX188400, was 99.71% and with the Karp reference sequences M31887 and LS398548, it was 99.43%. Two out of four sequences in the Karp-like clade were found identical and the other two (OR887446 and OR887441) differed by 2 and 1 nucleotide substitution, respectively. One out of 14 sequences were 100% identical to the Kawasaki-related ON156004 and the other Kuroki, Boryong, Gilliam and Kawasaki sequences showed 96.86% identity (Table 1).

3.2. TSA 56 kDA Gene VDI-III Sequence Analysis

All the 14 GroEL PCR positive samples were subjected to amplification of TSA 56 kDA variable domain-I to III. During the post-PCR processing, two different sized fragment patterns of the PCR amplified 56 kDA gene were recorded. Out of 14, a total of 12 TSA 56 kDa DNA sequences were generated and subjected to multiple sequence analysis in MEGA.11 software. DNA sequences having the query coverage of 99-100% with 98.3% to 99.9% nucleotide identity were retrieved from GenBank database and a conventional as well as the radiation distance phylogenetic tree was constructed. The resulting dendrogram with the 56 kDa sequences retrieved from GenBank showed a distinct cluster of Ot-TJTN and another Karp-like genotype. The one Gilliam type Ot sequence obtained in GroEL PCR did not amplify in the 56 kDa PCR assay. Among the twelve 56 kDA sequences, eight sequences of the new Ot-TJTN formed as a separate monophyletic cluster with the Ot whole genome study, complete and partial TSA 56 kDa sequences from Vellore, South India and the complete 56kDA coding sequence from Taiwan (CP142420) [27]. The TA716 and Kato genotypes are a common ancestral group to the newly emerged lineage. The other four Karp-like sequences (PQ381691, PQ381693, PQ381697 and PQ381700) formed as a monophyletic group along with the reference strain and sequences from India, Taiwan, China and Cambodia. Individual DNA sequences were analyzed in BLAST and information such as the GenBank accession numbers, length of query coverage, nucleotide identity, source of isolation of Ot and the country from where the sequence was deposited were collated and presented in Table 2.

3.3. TSA 56 kDA Gene VDI-III Similarity Plot Analysis

Only the eight O. tsutsugamushi 56 kDa gene sequences identified as a new cluster in this study were subjected to Simplot analysis with the 56 kDa strain specific reference sequences. The Similarity plot values -2.18, -3.06 and -5.06 indicated that the strain TA716-Thai isolate is related to 56 kDa Ot sequences generated in this study, followed by distantly related the Karp and Kato, respectively (Figure. 3a). The result of the bootscan plot demonstrated the event of possible recombination between TA716 like strain and Kato reference sequences at the nucleotide position 373 and 655 of the 56 kDa variable domain I-III (VDI to VDIII) (Figure. 3b).

Figure 3. (A) Similarity analysis of Simplot using the new Ot-TJTN TSA 56 kDA sequences as query. TA716, Karp and Kato reference sequences shows a overlapping curves in D-Ι against the new genotype (Ot-TJTN); (B). SimPlot bootscan graphs with reference sequences. The genotype of Ot-TJTN was assigned as query. The recombination breakpoint is clear between TA716 like at the nucleotide position 373 and 655 of the Ot TSA variable domain I to III.

Figure 3. (A) Similarity analysis of Simplot using the new Ot-TJTN TSA 56 kDA sequences as query. TA716, Karp and Kato reference sequences shows a overlapping curves in D-Ι against the new genotype (Ot-TJTN); (B). SimPlot bootscan graphs with reference sequences. The genotype of Ot-TJTN was assigned as query. The recombination breakpoint is clear between TA716 like at the nucleotide position 373 and 655 of the Ot TSA variable domain I to III.

4. Discussions

Scrub typhus (St), one of the neglected tropical diseases, has recently been reported to be escalating case incidences in endemic areas. It is considered to be an emerging public health problem globally [5,7]. In India, even though the serological evidence of the circulation of O. tsutsugamushi genotypes has been well documented [28], the available information on the genotypes of Ot in South India is very scanty. Generating the genotypic information from the particular region is important since it will help to design the diagnostic panel as well as the development of a vaccine. Different molecular markers are being used for the identification of the Ot genotypes [9]. However, to our limited knowledge, this is the first report on exploiting the GroEL and 56 kDa genes from south-eastern part of India, to screen, analyze and determine the phylogenetic relationship of the Ot recovered from human St cases. Among the two markers used in this study, the former encodes 60-kda heat-sock protein, which is more diverse, and highly conserved among O. tsutsugamushi but sufficiently variable to differentiate the genus Orientia and Rickettsia, and is suitable for the reconstruction of a phylogenetic tree [21] and, the latter is suitable and widely used for genotyping of Ot since it has four variable domains [14,18,22]. The 56 kDA TSA gene is mainly responsible for the elicitation of host immune response and neutralizing antibodies. [29]. The point mutation and genetic recombination in the TSA gene play a major role in determining the genotypic variation more accurately than the multilocus sequence typing (MLST) method (18).

In our earlier study, the GroEL gene sequence-based phylogenetic analysis has demonstrated the circulation of three distinct genotypes of Ot from this region [30]. Interestingly, nine out of 14 sequences generated in this study were found to have higher the nucleotide identity with the Ot DNA sequences available in the NCBI, GenBank database. Among the identical Ot DNA sequences retrieved, the following GroEL sequences with GB Ac Nos. ON156002, CP166954, CP166956 and PP754991 generated from human and rodent (unpublished) samples respectively, from Tamil Nadu, south India [30] are 100% identical, followed by EF551288 (TA716 Thailand) 97.71%, UT76 China) and CP044031 (Wji/2014 China) 97.41%. Six out of nine sequences in the new Ot-TJTN strain generated were found to be 100% identical, which probably indicates the close ancestral relationship of this isolate (Figure 4). However, this needs further confirmation.

In the other GroEL sequence cluster, four out of 14 sequences generated were found to be 99.14% to 99.71% of nucleotide identities with Karp-like sequences from India, Thailand, Japan, the United Kingdom, and the USA. Among the closely related sequences retrieved from the genome database, the flowing sequences, viz ON156000 VCRC/India, EF551309 Thailand, and JX188400 Japan were found to have the highest nucleotide identity of 99.71% (Table 1). Interestingly, out of the four Karp-like sequences, one (GB No. OR887446) formed a distinct sub-clade, due to the presence of two nucleotide substitutions. It indicates the circulation of multiple clones in this area may be due to recombination events occurring between strains of the same genotypes [31]. In the other clade, a single sequence, out of 14 sequences was found to have a nucleotide identity of 96.29% to 100% with the Kuroki-Gilliam-Boryong-Kawasaki DNA sequences retrieved from the genome database (Figure 2a). Among them, the Kuroki-Gilliam (ON156004) strain detected from human cases in the Theni district, South India belonging to the earlier study was 100% identical [30]. The remaining three closely identical sequences, namely, the Kuroki (LX188394), Boryong (AM494475) and Kawasaki (AY191587) from Japan and South Korea, had 96.86% nucleotide identity (Table 1). Only the Gilliam-type Ot sequences generated from this area were 100% identical. This result further indicates the introduction and circulation of a strain derived from common ancestry that is circulating in the south-western district of Theni, which is situated about 250 km apart from present study area in south-eastern district of Thanjavur.

The intra-strain mean diversity (d) was found to be higher (0.004) in the Karp-related group when compared with the new Ot-TJTN (0.001) and Kuroki (0.000) groups. The intra-mean diversity indicates that the Ot-TJTN strain is less diverse than the Karp-related strain. A single Kuroki-Gilliam-type sequence from this study (OR887445) was grouped with a Kawasaki-related sequence (ON156004 VCRC/India); no nucleotide diversity was detected within this group. Further, the inter-strain mean diversity was calculated after including reference sequences [21] and closely related sequences retrieved from the genome database in the appropriate group. The results showed a mean diversity of d=0.031 between Ot TJTN-Karp-like & Ot-TJTN-Kato-related strains and the highest diversity of d=0.052 was calculated between TJTN and Kuroki-Gilliam types. These d values indicate the independent circulation of the TJTN strain in this region.

Since the TSA 56 kDA protein is immunodominant, the intragenic recombination of this gene results in the genotypic variation in O. tsutsugamushi . The 56 kDA TSA based genotyping has been applied for determination of circulating strains in a particular geographic region [22]. In this study, analyses of hyper variable domain I-III (VD-I to VD-III) demonstrated the circulation of a new strain. The length of the DNA fragments was 650 bp and 692 bp for the Ot-TJTN and Karp-like sequences, respectively. Deletion of 42 nucleotides in the new genotypic clade resulted in amino acid variation when compared with the related reference strains included for analysis and Karp-like sequences detected in this study (Figure 5).

TSA 56 kDA amino acid similarity analysis in BLASTP showed 98% to 100% with the Ot strain sequences submitted from north Tamil Nadu (XDE91401- XDE91404) followed by 86% identity with the TSA sequences from Japan (AA072147). The neighbor-joining conventional and radiantion dendrograms clearly showed the two genotypic clusters (Figure 6).

The one with eight sequences belonged to a new genotype and the other was a Karp-like clade. The retrieved TSA 56 kDA DNA sequences (CP166954, CP166956) from north Tamil Nadu, India, which were isolated from humans were matched with 100% nucleotide identity with the designated new isolate Ot-TJTN. The other Ot 56 kDa sequences generated from mites (Chigger) collected in Vellore (MW604717), Hosur (MW604716), India; Taiwan (GU120144) were matched at 99.8% and 92.0%, respectively. Four sequences, especially from Taiwan, SEA, showed close identity with query coverage and nucleotide identity of 91.0%-99.0% and 92.0-98.5%, respectively. The eight retrieved sequences available without genotypic information from India showed nucleotide identity of 99.6%-100% with Ot-TJTN isolate. The other four Karp-like sequences showed 98.2%-99.8% nucleotide similarity with Ot 56 kDa sequences generated from chiggers and humans in India, China, Taiwan, Cambodia. The Karp-like clades exhibit a higher percentage of query coverage and similarity from different countries in Southeast Asia. On the other hand, the new genotype showed a 100% similarity to the isolates from this geographical region and 98.5%, 92.6% identity with Ot isolated from humans and chiggers from Taiwan, respectively. The above analysis confirms the generation of different clones due to genetic recombination in the chigger mites and maintenance and circulation of region-specific strains [22,32].

Simplot analysis for testing of the nucleotide similarity and bootscan analysis for detection of a recombination event between the study sequences (designated as query sequences) and prototype reference sequences showed a partial overlapping curve in the TA716 and Kato strain at the VD-I. The bootscan analysis showed a possible recombination event between TA716 and Kato at the nucleotide position 373 and 655 of the TSA 56 kDA gene. The Simplot test supports the result of phylogenetic analysis (Figure 3b) in concordance and the detected new cluster is genetically related to the Kato and TA716 strains but divergent to the extent of the new monophyletic group.

The results of the earlier studies have demonstrated the dominant occurrence of Ot belonging to the Kato-like strain in many parts of India [14,17]. However, the evidence of the present study shows, that the predominant distribution of that new Ot-TJTN genotype might have evolved from the Kato/TA716-genogroup [15,16]. In India, the results of 56 kDa TSA gene sequence-based genotype studies have been analyzed and found the TA678 strains with close genetic similarity to the Koto and TA716 genotypes [14,16]. Interestingly, the Ot-TJTN strain has been isolated in rodents (P754990, P754991 unpublished) and human cases (ON156002, ON156003) from Theni districts matched 100% nucleotide identity with GroEL sequences in this study [30]. Notably, the nine GroEL and eight the 56 kDa TSA sequences generated from this study formed a separate cluster, and they did not cluster with either the Koto genotypes or Koto reference sequences available in the GenBank database. The higher GroEL nucleotide percentage (97.71%) identity was obtained with EPW1038 of the TA716 (EF551288) Thailand reference strain next to the ON156002 (100%) generated from this geographic region. Similarly, the 56 kDA sequences too did not cluster with the sequences belonging to different genotypes except the 56 kDA sequences generated from this geographic region, followed by Taiwan isolates with 92.6% to 92.9% (CP142420, MW495758 and MW495734). The result of these analyses indicates that all the 9 GroEL and eight 56 kDA Ot-TJTH sequences are undescribed but closely related to the TA716 Thai isolate (Figure 4 & Figure 5). Further, this implies that the Ot-TJTN strain has been reported for the first time in this area, and it is worth understanding their possible regional-specific distribution, if any.

Worldwide, in recent days, molecular genotype studies have revealed the distribution of different genotypes. Even though information on the serological-based identification of the Ot genotype from different states of India is available, due to the lack of genome-based studies across the different geographical regions, the strain identification results were concluded, and the close identity was drawn mainly in comparison with major prototypes of Karp, Kato, Gilliam, Boryong, Ikeda, Kawasaki, Kuroki, Saitama, and Shimokoshi available in genome database. In this study, the predominant occurrence of Ot-TJTN was confirmed with two reliable molecular markers. Interestingly, Ot genotypes identified in this study show close identity with strains in the SEA region (Figure 2), which confirms the circulation of ancestral haplotypes throughout the SEA continent [22]. Therefore, these strains may be considered to be specific to this geographical region because the antigenic heterogeneity is restricted to the geographic epicenter of St endemicity [33].

This laboratory-based study indicates the rich diversity of Ot in these areas, causing significant public health problems. Similarly, our recent observation [30] in the Theni district, Tamil Nadu, showed the circulation of multiple genotypes of Ot and causing clinical cases in humans. These results further indicate the necessity of conducting in-depth studies to understand the Ot distribution in different parts of the state, since the same strains are circulating in human cases, and rodent mites in Thanjavur and adjacent districts in this region.

5. Conclusions

Scrub typhus is a re-emerging public health problem in India. The results of this preliminary study demonstrate a rich genetic diversity in O. tsutsugamushi circulating in and around Thanjavur district, Tamil Nadu, India. The phylogenetic analysis based on the partial genome sequence of the GroEL and 56 kDA genes shows the predominant circulation of a novel strain Ot-TJTN along with Karp-like and Kuroki-Gilliam-type strains. The result of this study urges the need to conduct a genome-based in-depth study on the scrub typhus pathogen, considering the increased scrub typhus cases [4,5], neurological complications, emergence of new species [11,12,34] seasonal distribution of the strain [35] and the evolving nature of the pathogen. These findings would eventually help to include the new antigenic types in the immunodiagnostic panel. Also, the outcome of this study warrants further studies on understanding the mechanisms of the ongoing transmission dynamics of Ot stains in vector mites, rodent hosts, and humans for implementing appropriate control measures, particularly one health approach in the future.