Phylogenetic Position of Haemaphysalis kashmirensis and Haemaphysalis cornupunctata, with Notes on Rickettsia spp.

Despite high diversity in the Oriental region, ticks of the genus Haemaphysalis have been neglected regarding their genetic data and vector potential. This study aimed to genetically characterize three species of the genus Haemaphysalis: Haemaphysalis cornupunctata, Haemaphysalis kashmirensis, and Haemaphysalis montgomeryi infesting goats and sheep, and Rickettsia spp. associated with these tick species in the Hindu Kush Himalayan range of Pakistan. Altogether, 834 ticks were collected by examining 120 hosts including goats (64/120, 53.3%) and sheep (56/120, 46.6%), in which 86 (71.6%) hosts were found to be tick-infested. The morphologically identified ticks were subjected to DNA extraction and PCR for the amplification of partial 16S rDNA and cox fragments. Rickettsia spp. associated with the collected ticks were detected through the amplification of gltA, ompA and ompB partial fragments. The 16S rDNA of H. cornupunctata and H. montgomeryi showed a maximum identity of 100% with the sequences of the same species, whereas the 16S rDNA of H. kashmirensis showed the highest identity of 93–95% with Haemaphysalis sulcata. The cox sequence of H. montgomeryi displayed 100% identity with the same species. In comparison, the cox sequences of H. cornupunctata and H. kashmirensis showed maximum identities of 87.65–89.22% with Haemaphysalis punctata and 89.34% with H. sulcata, respectively. The gltA sequence of Rickettsia sp. from H. kashmirensis showed the highest identity of 97.89% with Rickettsia conorii subsp. raoultii, while the ompA and ompB fragments from the same DNA samples revealed 100% and 98.16% identity with Rickettsia sp. and “Candidatus Rickettsia longicornii”, respectively. Another gltA sequence amplified from H. montgomeryi ticks showed 100% identity with Rickettsia hoogstraalii, while the attempts to amplify ompA and ompB for R. hoogstraalii were unsuccessful. In the phylogenetic tree, the 16S rDNA of H. cornupunctata clustered with the corresponding species while its cox clustered with H. punctata. Both 16S rDNA and cox sequences of H. kashmirensis clustered with H. sulcata. The gltA sequence of Rickettsia sp. was clustered individually in the spotted fever (SF) group of Rickettsia, while the gltA sequence of R. hoogstraalii was clustered with the same species in the transition group of Rickettsia. In the SF group, the rickettsial ompA and ompB sequence clustered with undetermined Rickettsia sp. and “Candidatus Rickettsia longicornii”, respectively. This is the earliest study regarding the genetic characterization of H. kashmirensis. This study indicated that ticks belong to the genus Haemaphysalis have the potential of harboring and/or transmitting Rickettsia spp. in the region.


Introduction
Ticks are notorious ectoparasites infesting the majority of terrestrial and semiterrestrial vertebrates [1]. Ticks are cosmopolitan in distribution, particularly prevalent in tropical and subtropical regions [1]. In addition to anemia, a reduction in dairy products and meat, ticks can also transmit numerous infectious agents including viruses, bacteria, protozoa, and helminths to vertebrate hosts [2].
Haemaphysalis (Acari: Ixodidae) is the second largest genus of ixodid ticks comprising approximately 176 species [1]. Haemaphysalis ticks inhabit various landscapes mainly in the Oriental, Afrotropical, and Palearctic regions, while some of its members have been recorded in Australasia (Australia, New Zealand, and New Guinea), infesting a variety of free-roaming and domestic animals [1,3]. Some species of Haemaphysalis have been identified as potential vectors for various protozoan, viral, and bacterial pathogens to humans, domestic, and wild animals [4,5]. Haemaphysalis ticks such as Haemaphysalis longicornis, Haemaphysalis concinna, Haemaphysalis qinghaiensis, and Haemaphysalis flava have been associated with various Rickettsia spp. including Rickettsia japonica, R. conorii, Rickettsia monacensis, Rickettsia heilongjiangensis, and several undetermined Rickettsia spp. [6][7][8][9]. Moreover, the number of Rickettsia spp. associated with Haemaphysalis ticks are continuously increasing because of the advancement in molecular approaches [10].
The livestock hosts in Pakistan are known to be infested by a variety of Haemaphysalis ticks [16,19,20]. Genetic characterization of closely related Haemaphysalis ticks such as H. kashmirensis and H. cornupunctata is important for the surveillance and control of these ticks in the Oriental region including Pakistan, which is an epidemic hotspot for these parasites [14,16,18]. The capabilities of Haemaphysalis ticks as a vector for Rickettsia spp. have been neglected. These ticks need to be molecularly screened for associated Rickettsia spp. This study was aimed to the molecular characterization of H. kashmirensis and H. cornupunctata ticks infesting goats and sheep, and associated Rickettsia spp. in the HKH range of Pakistan.

Study Sites
The present study was conducted in four districts, namely: Dir Lower (34 • Figure 1).

Ethical Approval and Consent
The design of the present study was approved by the Advanced Study and Research Board (Dir/A&R/AWKUM/2022/9396) and the Ethical Committee of the Faculty of Chemical and Life Sciences, Abdul Wali Khan University Mardan, Pakistan. All animal owners were informed verbally, and approval was taken before observing their hosts.

Ethical Approval and Consent
The design of the present study was approved by the Advanced Study and Research Board (Dir/A&R/AWKUM/2022/9396) and the Ethical Committee of the Faculty of Chemical and Life Sciences, Abdul Wali Khan University Mardan, Pakistan. All animal owners were informed verbally, and approval was taken before observing their hosts.

Tick Collection and Morphological Identification
Villages of the four selected districts were randomly visited between March 2021 to February 2022 for tick collection. Small ruminants including goats and sheep were examined for ticks. When found, ticks from each host were separately collected in labeled tubes using a fine tweezer. The necessary information regarding collection sites, host type and gender, collection date, and environmental conditions (temperature and humidity) were noted. Collected ticks were rinsed with distilled water followed by 70% ethanol to remove contaminants. Subsequently, standard identification keys [12][13][14]21] were used to identify

Tick Collection and Morphological Identification
Villages of the four selected districts were randomly visited between March 2021 to February 2022 for tick collection. Small ruminants including goats and sheep were examined for ticks. When found, ticks from each host were separately collected in labeled tubes using a fine tweezer. The necessary information regarding collection sites, host type and gender, collection date, and environmental conditions (temperature and humidity) were noted. Collected ticks were rinsed with distilled water followed by 70% ethanol to remove contaminants. Subsequently, standard identification keys [12][13][14]21] were used to identify tick species based on observation of morphological features such as shape of capitulum, cornua and palp articles, cervical and lateral grooves, size and shape of coxa spurs, number of festoons, and shape of genital aperture. These morphological observations were performed under a stereomicroscope (SZ61, Olympus Corporation, Tokyo, Japan).

Statistical Analyses
Necessary tick data recorded from the four districts were compiled and arranged in spreadsheets using Microsoft Excel V. 2016 (Microsoft Office 365 ® ). Prevalence (no. of infested hosts×100/total no. of examined hosts), mean abundance (total no. of collected ticks/total no. of examined hosts), and mean intensity (total no. of collected ticks/no. of infested hosts) was determined. The necessary climate data for each month or seasonal were obtained via mean temperature ( • C), relative humidity (%), and total rainfall (mm) (climate-data.org; accessed on 15 November 2022).

DNA Extraction and PCR Amplification
Altogether, 108 tick specimens (three males, three females, and three nymphs per species per district) were randomly selected for molecular analysis. The morphologically identified ticks were individually dissected with a sterile blade and ground with a hygienic pestle in 1.5 mL Eppendorf tubes. The ground samples were individually subjected to the phenol-chloroform protocol [22] for DNA extraction. The extracted genomic DNA was quantified by using NanoDrop (Nano-Q, Optizen, South Korea).
The extracted DNA was subjected to conventional PCR to amplify 16S rDNA (460 bp) and cox (710 bp) fragments of tick species and citrate synthase (gltA), outer membrane protein subunit A (ompA), and outer membrane protein subunit B (ompB) fragment for any Rickettsia spp. Each PCR reaction comprised 25 µL volume: 8.5 µL of PCR water "nuclease free", 1 µL of each primer at a concentration of 10 pmol/µL, 2 µL of extracted DNA (50-100 ng/µL), and 12.5 µL DreamTaq green MasterMix (2X) ( Table 1). In each PCR reaction, PCR water was used as a negative control while Rhipicephalus microplus and Rickettsia massiliae DNA were taken as positive control for ticks and Rickettsia spp., respectively. PCR products were analyzed by horizontal electrophoresis in 2% agarose gel and examined under ultraviolet light of a Gel Documentation System (BioDoc-It™ Imaging Systems, Upland, CA, USA).

Sequencing and Phylogenetic Analysis
Amplicons were purified with GeneClean II Kit (Qbiogene, Illkirch, France) following the manufacturer's protocol and sequenced bi-directionally (Macrogen, Inc., Seoul, South Korea) via the Sanger-based sequencing method. The obtained sequences were trimmed by removing the poor-quality regions followed by obtaining a consensus sequence in SeqMan v 5.00 (DNASTAR, Inc., Madison, WI, USA). Maximum identity sequences were retrieved from GenBank using the Basic Local Alignment Search Tool (BLAST) [28] on the National Center for Biotechnology Information (NCBI) user interface. They were aligned with the obtained sequences in BioEdit Sequence Alignment Editor v. 7.0.5 [29] using CLUSTAL W multiple alignments [30]. The Neighbor-Joining method with the Kimura 2parameter model was applied to construct phylogenetic trees with 1000 bootstrap replicates in Molecular Evolutionary Genetic Analysis (MEGA-X) software [31]. The coding (cox, gltA, ompA, and ompB) nucleotide sequences were aligned by MUSCLE [32]. The final positions in the dataset comprised the obtained sequence of each fragment.  Table 2.

Ticks, Hosts, and Seasonal Data
Ticks were collected from 86 out of 120 examined hosts (goats and sheep) with an overall 71.3% prevalence. Goats were observed with a high prevalence (49/64, 76.6%) compared to sheep (37/56, 66.1%). An overall mean abundance of 7.0 ticks/examined host was noted, while the mean intensity was 9.7 ticks/infested host. Goats were found with higher mean abundance and mean intensity (8.1 and 10.6, respectively) compared to sheep (5.6 and 8.6, respectively). Seasonal parameters, including temperature, humidity, and rainfall with tick parameters such as the prevalence of infestation, mean abundance, and mean intensity are shown in Figure 2.
quences were retrieved from GenBank using the Basic Local Alignment Search Tool (BLAST) [28] on the National Center for Biotechnology Information (NCBI) user interface. They were aligned with the obtained sequences in BioEdit Sequence Alignment Editor v. 7.0.5 [29] using CLUSTAL W multiple alignments [30]. The Neighbor-Joining method with the Kimura 2-parameter model was applied to construct phylogenetic trees with 1000 bootstrap replicates in Molecular Evolutionary Genetic Analysis (MEGA-X) software [31]. The coding (cox, gltA, ompA, and ompB) nucleotide sequences were aligned by MUSCLE [32]. The final positions in the dataset comprised the obtained sequence of each fragment.  Table 2.

Ticks, Hosts, and Seasonal Data
Ticks were collected from 86 out of 120 examined hosts (goats and sheep) with an overall 71.3% prevalence. Goats were observed with a high prevalence (49/64, 76.6%) compared to sheep (37/56, 66.1%). An overall mean abundance of 7.0 ticks/examined host was noted, while the mean intensity was 9.7 ticks/infested host. Goats were found with higher mean abundance and mean intensity (8.1 and 10.6, respectively) compared to sheep (5.6 and 8.6, respectively). Seasonal parameters, including temperature, humidity, and rainfall with tick parameters such as the prevalence of infestation, mean abundance, and mean intensity are shown in Figure 2. The graph shows the variation in prevalence, mean abundance, and mean intensity related to mean temperature, mean relative humidity, and the total rainfall during this study. The graph shows the variation in prevalence, mean abundance, and mean intensity related to mean temperature, mean relative humidity, and the total rainfall during this study.

Molecular Analysis
Altogether, 72 sequences were obtained (one 16S rDNA and one cox for one female, one male, and one nymph of each tick species) for ticks. Trimmed sequences for H. cornupunctata, H. kashmirensis, and H. montgomeryi were 16S rDNA (395 bp) and cox (654 bp). The 16S rDNA of H. cornupunctata showed 100% identity with the same species from Pakistan. In contrast, its cox showed 87.65-89.22% identity with the same subgenus species: H. (Aboimisalis) punctata from Portugal, Iran, France, China, Hungary, and Romania. The BLAST outcome for the 16S rDNA of H. kashmirensis revealed 93-95% identity with the same subgenus species: H. (Herpetobia) sulcata reported from France, Turkey, China, and Pakistan, compared to cox which showed 89.34% identity with H. sulcata reported from China and Iran. The sequences obtained for H. montgomeryi were 100% identical to the sequences of the same species reported from Pakistan [16]. Thus, such sequences were excluded from further analysis.
For Rickettsia spp., bidirectional sequences were obtained for each amplified fragment. The consensus sequences of gltA (380 bp), ompA (504 bp), and ompB (449 bp) were amplified from H. kashmirensis, while gltA (340 bp) was amplified from H. montgomeryi. Rickettsial gltA sequence (380 bp) amplified from H. kashmirensis showed 97.89% identity with R. conorii subsp. raoultii from Brazil, China, and France, while another rickettsial gltA consensus sequence from H. montgomeryi displayed 100% identity with R. hoogstraalii reported from Italy. The ompA amplified from the same samples of H. kashmirensis showed 100% identity with Rickettsia sp. from China and Algeria. The ompB sequence of the corresponding sample showed 98.16% identity with "Candidatus Rickettsia longicornii" reported from China. The rate of infection was 15.74% (17/108) recorded for Rickettsia sp. (based on gltA, ompA, and ompB) followed by 5.56% (6/108) for R. hoogstraalii (based only on gltA). The details regarding the infection rate and the number of sequences obtained for each Rickettsia spp. are shown in Table 2.
The obtained 16S rDNA sequences were submitted to GenBank under the accession numbers: OQ024373 (H.

Phylogenetic Analysis
In a phylogenetic tree based on 16S rDNA, the H. cornupunctata sequence clustered with the same species reported from Pakistan (ON911369), whereas this species grouped in a sister clade with the species of the same subgenus: H. punctata reported from China (NC062064, MF002566, and MG021187) and Turkey (KR870978) (Figure 3 (Figure 4).
Phylogenetic tree based on rickettsial gltA, the Rickettsia sp. detected in H. kashmirensis clustered individually with the sequence of Rickettsia spp. of the spotted fever group. In contrast, the R. hoogstraalii detected in H. montgomeryi was clustered with the same species reported from Italy (KY418024-KY418025) ( Figure 5). The rickettsial ompA fragment detected in H. kashmirensis clustered with the Rickettsia sp. reported from China (MT361020, and grouped in a sister clade with "Candidatus Rickettsia longicornii" and "Candidatus Rickettsia jingxinensis", while the rickettsial ompB fragment clustered with "Candidatus Rickettsia longicornii" reported from China (MN026546, MK620854, MG906675, and MT511089) (Figure 7). (MT361020, MG228270, and MN644903) and Algeria (MZ064523-MZ064524 and JN943296) ( Figure 6) and grouped in a sister clade with "Candidatus Rickettsia longicornii" and "Candidatus Rickettsia jingxinensis", while the rickettsial ompB fragment clustered with "Candidatus Rickettsia longicornii" reported from China (MN026546, MK620854, MG906675, and MT511089) (Figure 7).     Neighbor-Joining phylogenetic tree based on gltA sequences of Rickettsia sp. and R. hoogstraalii. The gltA sequences of Rickettsia australis and Rickettsia akari were employed as an outgroup. All sequences have been denoted by their GenBank accession numbers, followed by species name and country name. The bootstrap values (1000-replications) are shown at each node. The sequences (OQ160793 and OQ160792) of the present study have been marked with bold and underlined fonts.

Discussion
As the geo-climatic conditions of the Oriental region including Pakistan suit the flourishment of Haemaphysalis ticks, the largest diversity of these ticks has been reported in this region [16]. The second most diverse genus of hard ticks (Ixodidae), Haemaphysalis comprises~173 tick species globally. In Pakistan, 13 Haemaphysalis spp. have been reported; however, genetic data of these ticks are limited [16]. Despite the huge diversity, genetic data regarding the genus Haemaphysalis ticks and associated Rickettsia spp. have been largely neglected. For this purpose, Haemaphysalis ticks were collected from goats and sheep in northern Pakistan, where several species of this genus are considered endemic [16,18,19]. Herein, the collected ticks were morpho-molecularly identified as H. cornupunctata, H. kashmirensis, and H. montgomeryi. The genetic characterization based on 16S rDNA and cox partial sequences of H. kashmirensis, and cox sequence for H. cornupunctata was achieved for the first time. An undetermined Rickettsia sp. based on gltA, ompA, and ompB sequence was molecularly characterized in H. kashmirensis, whereas R. hoogstraalii based on only gltA sequence was detected in H. montgomeryi.
The surveyed areas are part of the HKH mountain range, which has been considered as one of the most important biodiversity hotspots [18,33]. Previously, H. kashmirensis, H. cornupunctata, and H. montgomeryi have been reported from many locations of the HKH range, which spans different territories of the Indian subcontinent, including Pakistan [15,17,18], Kashmir and India [14], and Afghanistan [34]. These findings suggest that HKH mountain regions have a great diversity of the Haemaphysalis ticks owing to the abundance of suitable hosts and conducive climate conditions.
Haemaphysalis kashmirensis was found less in number than H. cornupunctata and H. montgomeryi, which could be due to the association of the later species with hosts other than goats and sheep in the family Bovidae [1]. Moreover, the adult ticks (female and male) were outnumbered by the immature ticks (nymphs), while no larval stage of any tick species was found on goats and sheep. Previous studies suggested that Agama tuberculata (Kashmir Rock Agama) is the main host of the nymphal and larval stages of H. kashmirensis, which is found in the HKH range [21]. Similarly, animals belonging to the Muridae, Herpestidae, Erinaceidae, Cricetidae, and Soricidae families have been recorded as the main hosts for the nymphal and larval ticks of H. cornupunctata and H. montgomeryi [1].
Haemaphysalis ticks have been reported as vectors for Rickettsia spp. including R. hoogstraalii and Rickettsia rhipicephali [5,35]. Herein, Rickettsia sp. was detected through gltA, ompA, and ompB, whereas R. hoogstraalii was detected only through gltA. The genetic characterization of R. hoogstraalii was also attempted through ompA and ompB; however, the amplifications of these fragments were unsuccessful. Amplification failures are common in the case of ompA, ompB that might be the lack of targeted genes, as shown in the transition group Rickettsia or due to primer mismatching [36][37][38]. Rickettsia hoogstraalii has been detected in ticks of the genus Haemaphysalis such as H. sulcata (Cyprus and Italy), H. punctata (Italy), and Haemaphysalis parva (Turkey) [39][40][41]. The pathogenicity of R. hoogstraalii is poorly known [42].
Genetic analyses based on molecular markers such as 16S rDNA and cox genes are extremely helpful in unveiling the systematics of ticks and phylogenetic positioning [43][44][45][46][47][48][49]. In the phylogenetic tree based on 16S rDNA and cox of H. cornupunctata, this species appeared in a monophyletic branch with H. punctata from Turkey, France, China, Iran, Romania, Portugal, and Hungary. Previously, these two tick species (H. cornupunctata and H. punctata) have been assigned to the same subgenus (Aboimisalis) on a morphological basis [11]. In a phylogenetic tree based on 16S rDNA and cox of H. kashmirensis, this specie clustered in a monophyletic clade with H. sulcata previously reported from Pakistan, China, Turkey, France, and Iran. On a morphological basis, H. sulcata and H. kashmirensis have been placed in the same subgenus, Herpetobia [11]. A phylogenetic clustering of these species with different closest species of the same subgenus could be associated with the missing genetic data of corresponding species in the GenBank. The phylogenetic analysis based on rickettsial gltA showed that the Rickettsia sp. detected in H. kashmirensis belonged to the SF group, whereas R. hoogstraalii detected in H. montgomeryi belonged to the transition group. The phylogenetic tree based on ompA and ompB, obtained from the same sample in which Rickettsia sp. was detected, validated the gltA-based phylogenetic analysis for Rickettsia sp.

Conclusions
This study contributes to the missing information regarding the genetic data of some Haemaphysalis ticks, especially H. kashmirensis, which was genetically characterized for the first time. The relationship of H. cornupunctata with subgenus Aboimisalis and H. kashmirensis with subgenus Herpetobia established on a morphological basis was confirmed through molecular-based phylogenetic analysis. Furthermore, a Rickettsia sp. was molecularly assessed in H. kashmirensis, whereas R. hoogstraalii was detected in H. montgomeryi. This study may assist in understanding the identification, evolutionary history, and molecular epidemiology of Haemaphysalis ticks and associated Rickettsia spp. Further studies should genetically characterize and evaluate the pathogenicity of these Rickettsia spp.