Spatio-Temporal Patterns of Ticks and Molecular Survey of Anaplasma marginale, with Notes on Their Phylogeny

Hard ticks (Ixodida: Ixodidae) are medically important ectoparasites that feed on all classes of terrestrial vertebrates. Recently, we molecularly characterized hard ticks and associated Anaplasma spp. in the northern and central regions of Khyber Pakhtunkhwa (KP), Pakistan; however, this knowledge was missing in the southern regions. This study aimed to investigate tick prevalence, host range, genetic diversity, and molecular survey of Anaplasma spp. in a wide range of tick species in two distinct physiographic regions of southern KP. A total of 1873 hard ticks were randomly collected from 443/837 hosts (cattle, Asian water buffaloes, horses, goats, sheep, dogs, and camels) in Lakki Marwat, Bannu, and Orakzai districts of KP. Overall, 12 tick species were morphologically identified, among which Hyalomma dromedarii was the most prevalent species (390/1873, 20.9%), followed by Hy. anatolicum (294, 15.7%), Rhipicephalus microplus (262, 14%), Hy. scupense (207, 11.1%), R. sanguineus (136, 7.3%), R. turanicus (121, 6.5%), Haemaphysalis cornupunctata (107, 5.7%), R. haemaphysaloides (110, 5.9%), Ha. montgomeryi (87, 4.6%), Hy. isaaci (58, 3.1%), Ha. bispinosa (54, 2.9%), and Ha. sulcata (47, 2.5%). The extracted DNA from a subset of each tick species was subjected to PCR to amplify cox1 or 16S rRNA sequences of ticks and 16S rRNA sequences of Anaplasma spp. The tick cox1 sequences showed 99–100% identities with the sequences of the same species, whereas 16S rRNA sequences of R. turanicus, Ha. montgomeryi and Ha. sulcata showed 97–100% identities with the corresponding species. The 16S rRNA sequence of Ha. cornupunctata showed 92% identity with the species from the same subgenus, such as Ha. punctata. The 16S rRNA sequence of Anaplasma spp. showed 100% identity with Anaplasma marginale. Moreover, 54 ticks were found positive for A. marginale with a total infection rate of 17.2%. The highest infection rate was recorded in Hy. dromedarii (31.1%) and the lowest in each R. haemaphysaloides and R. sanguineus (20%). All the cox1 or 16S rRNA sequences in phylogenetic trees clustered with the same species, except Ha. cornupunctata, which clustered with the Ha. (Aboimisalis) punctata. In this study, Ha. cornupunctata was reported for the first time at the molecular level. The genetic characterization of ixodid ticks and molecular detection of associated A. marginale will assist in the epidemiological surveillance of these parasites in the region.

The pathogenic agents of anaplasmosis are highly prevalent worldwide, particularly in tropical and subtropical regions [19]. These pathogens have a wide genetic range and adversely affect the livestock industry [20,21]. Knowing that this field of research has attracted attention, there are still very few available studies restricted to a few areas in Pakistan about the molecular data of ticks and Anaplasma spp. [2,8,13,18,[22][23][24]. Our previous study has demonstrated the molecular assessment of hard ticks and associated A. marginale collected from livestock hosts in the northern and central regions of Khyber Pakhtunkhwa (KP), Pakistan [24]. Still, similar studies are missing from the southern regions. This study aimed to investigate tick prevalence, genetic diversity, and molecular survey of associated Anaplasma spp. in a wide range of tick species in two distinct physiographic regions of southern KP.

Ethical Approval
Before this study, ethical approval was taken from the Advanced Studies and Research Board (Dir/A&R/AWKUM/2020/4871) of the Faculty of Chemical and Life Sciences, Abdul Wali Khan University Mardan, KP, Pakistan. Furthermore, written and/or oral consents were obtained from the animals' owners for tick collection.

Study Area
The current study investigated three districts of southern KP, including Lakki Marwat (32. . These districts belong to two distinct physiographic regions, one with a "hot semi-arid climate" (Bannu and Lakki Marwat) and the other with a "humid subtropical climate" (Orakzai). Based on the ecological zones, the former is mainly a "desert plain," and the latter is mainly a semi-arid piedmont. The geographic coordinates of each collection site were obtained using Global Positioning System (GPS) and loaded into a Microsoft Excel sheet to design a map using ArcGIS 10.

Tick Collection and Preservation
Tick collection was carried out from March 2019 to February 2020 with a regular visit to the study area once a month. Ticks were randomly collected using forceps from different vertebrate hosts, including cattle, Asian water buffaloes, horses, goats, sheep, dogs, and camels ( Figure 2). Tick specimens were rinsed with distilled water followed by 70% ethanol and were stored in 100% ethanol in properly labeled tubes for onward molecular experiments. During tick collection, the relevant information regarding collection date, host type, and place of collection of the ticks were noted.

Tick Collection and Preservation
Tick collection was carried out from March 2019 to February 2020 with a regular visit to the study area once a month. Ticks were randomly collected using forceps from different vertebrate hosts, including cattle, Asian water buffaloes, horses, goats, sheep, dogs, and camels ( Figure 2). Tick specimens were rinsed with distilled water followed by 70% ethanol and were stored in 100% ethanol in properly labeled tubes for onward molecular experiments. During tick collection, the relevant information regarding collection date, host type, and place of collection of the ticks were noted.

DNA Extraction and PCR
Before the genomic DNA extraction, ticks were washed with distilled water and dried on sterile filter paper. The ticks were crushed with sterilized pestles in 1.5 mL sterile Eppendorf tubes. Genomic DNA was extracted individually from each tick using the phenol-chloroform method according to the standard protocol. The DNA pellet was hydrated by adding 30 µL of nuclease-free water [34]. The quality and quantity of genomic DNA were determined through Nano-Q (Optizen, Daejeon, Korea).
By using reference primers and PCR conditions (Table 1), the extracted DNA was subjected to amplifying partial fragments of ticks cox1 and 16S rRNA genes and screened for 16S rRNA of Anaplasma spp. in Table 2 through a PCR. Each PCR reaction was prepared in a 20 µL reaction mixture and contained: 12 µL of DreamTaq MasterMix (Thermo Fisher Scientific, Inc., Waltham, MA, USA), 1 µL of each forward and reverse primers (10 µM), 2 µL (50 ng/µL) genomic DNA template and 4 µL PCR water (nuclease-free). The DNA of R. microplus and Rickettsia massiliae were used as positive controls for ticks and Anaplasma spp., respectively, while PCR water (nuclease-free) was used as a negative control. The amplified DNA was run on a 1.5% agarose gel, dyed with 2 µL ethidium bromide, and observed by a Gel documentation system (BioDoc-It™ Imaging Systems UVP, LLC, Upland, CA, USA).

DNA Extraction and PCR
Before the genomic DNA extraction, ticks were washed with distilled water and dried on sterile filter paper. The ticks were crushed with sterilized pestles in 1.5 mL sterile Eppendorf tubes. Genomic DNA was extracted individually from each tick using the phenol-chloroform method according to the standard protocol. The DNA pellet was hydrated by adding 30 µL of nuclease-free water [34]. The quality and quantity of genomic DNA were determined through Nano-Q (Optizen, Daejeon, Korea).
By using reference primers and PCR conditions (Table 1), the extracted DNA was subjected to amplifying partial fragments of ticks cox1 and 16S rRNA genes and screened for 16S rRNA of Anaplasma spp. in Table 2 through a PCR. Each PCR reaction was prepared in a 20 µL reaction mixture and contained: 12 µL of DreamTaq MasterMix (Thermo Fisher Scientific, Inc., Waltham, MA, USA), 1 µL of each forward and reverse primers (10 µM), 2 µL (50 ng/µL) genomic DNA template and 4 µL PCR water (nuclease-free). The DNA of R. microplus and Rickettsia massiliae were used as positive controls for ticks and Anaplasma spp., respectively, while PCR water (nuclease-free) was used as a negative control. The amplified DNA was run on a 1.5% agarose gel, dyed with 2 µL ethidium bromide, and observed by a Gel documentation system (BioDoc-It™ Imaging Systems UVP, LLC, Upland, CA, USA). • Count for fully, partially and unengorged.

DNA Sequencing and Phylogenetic Analysis
Purification of PCR products was performed using GeneClean II Kit (Qbiogene, Illkirch, France) following the manufacturer's protocol. A total of 64 (cox1 40 and 16S rRNA 24) amplified PCR products for ticks and 18 (3 from each Anaplasma positive tick species amplicons) for 16S rRNA Anaplasma spp. were submitted for bidirectional DNA sequencing (Macrogen, Inc., Seoul, South Korea). The sequences were cropped to remove the primers and poor reading regions through SeqMan V. 5 (DNASTAR). The obtained purified sequences were subjected to the Basic Local Alignment Search Tool (BLAST) [38] at National Center for Biotechnology Information (NCBI), and the homologous sequences were downloaded. These sequences were aligned with obtained sequences along with an outgroup in BioEdit Sequence Alignment Editor V. 7.0.5 (Raleigh, NC, USA) [39]. The phylogenetic trees were constructed by using the Maximum-Likelihood model (1000 bootstrap replicons) in Molecular Evolutionary Genetics Analysis (MEGA-X) [40].

Morphologically Identified Ticks
The morphological identification confirmed 12 tick species belonging to the three genera of hard ticks. Haemaphysalis species were only found in the Orakzai district, while we could not collect these species in the other two districts. The details of each tick species reported from the study area are provided (Figure 3).

Sequencing Analysis
From the extracted tick DNA of 12 tick species, the partial fragments of cox1 were amplified for eight tick species, whereas 16S rRNA was amplified for four tick species. Clean cox1 sequences were obtained from eight tick species: Hy. dromedarii (743

Phylogenetic Analysis
The phylogenetic tree for the cox1 sequences of Hy. dromedarii, Hy. anatolicum, Hy. scupense, Hy. isaaci, Ha. bispinosa, R. microplus, R. sanguineus, and R. haemaphysaloides were constructed combinedly with 49 sequences downloaded from NCBI based on the maximum identity. In the phylogenetic tree, the obtained cox1 sequences were clustered to the corresponding species reported from different countries, such as Hy. dromedarii from Egypt and Tunisia, Hy. anatolicum from India, Egypt, and China, Hy. scupense from France, Spain, China, and Turkey, Hy. isaaci from Pakistan, Ha. bispinosa from India and Bangladesh, R. microplus from Pakistan, India, and China, R. sanguineus from Iran and India, and R. haemaphysaloides from Pakistan, India, and China. In the case of 16S rRNA, the phylogenetic tree of R. turanicus, Ha. cornupunctata, Ha. montgomeryi and Ha. sulcata was constructed with 27 sequences downloaded from NCBI based on the maximum identity. In the phylogenetic tree, 16S rRNA sequences of R. turanicus, Ha. montgomeryi and Ha. sulcata clustered with the same species reported from Afghanistan, Pakistan, and China, while Ha. cornupunctata clustered with the species of the same subgenus Ha. (Aboimisalis) punctata reported from China, Turkey, Italy, Spain, and Portugal.
The obtained partial 16S rRNA sequence of A. marginale was uploaded to the GenBank (ON528757).

Discussion
Pakistan has an agrarian economy where agriculture contributes approximately 21% to gross domestic product (GDP) and 45% to the labor force [41]. Ticks pose severe threats to the livestock and economy of the country. Knowledge regarding molecular surveillance of ticks and A. marginale and their host range in different physiographic is essential for implementing adequate measures against these parasites in Pakistan. The present study

Discussion
Pakistan has an agrarian economy where agriculture contributes approximately 21% to gross domestic product (GDP) and 45% to the labor force [41]. Ticks pose severe threats to the livestock and economy of the country. Knowledge regarding molecular surveillance of ticks and A. marginale and their host range in different physiographic is essential for implementing adequate measures against these parasites in Pakistan. The present study was executed in two distinct physiographic regions in southern KP, Pakistan. The targeted areas were selected because ticks and tick-borne diseases are common in these regions but mainly remained unexplored, and to compare tick diversity in two regions that are geographically close but physiographically and climatically different. Herein, 12 tick species were morphologically and molecularly identified. Four tick species, including Ha. bispinosa, Ha. cornupunctata, Hy. dromedarii and Hy. isaaci were genetically characterized for the first time from Pakistan. Furthermore, the molecular survey was conducted to screen a subset of the collected 12 species for A. marginale, in which this pathogen was detected in six species. Among these species, A. marginale was detected for the first time in Hy. dromedarii, Hy. scupense, R. sanguineus and R. haemaphysaloides from Pakistan.
Environmental and climatic factors influence the distribution and prevalence of ticks within a specific region [42]. Previous studies considered Hyalomma spp. as successful ticks in harsh desert regions [43,44]. Similarly, as a larger proportion of the current study area was a desert plain, the genus Hyalomma was the most prevalent, followed by genus Rhipicephalus and Haemaphysalis. Herein, unlike [2,3,8], Hy. dromedarii was the most prevalent in the region owing to the screening of a larger number of camels compared to other hosts. According to the studies performed in the region [2,3,8], R. microplus and Ha. bispinosa were the most prevalent tick species in the genus Rhipicephalus and genus Haemaphysalis, respectively.
The highest prevalence of A. marginale occurs in those regions where R. microplus is endemic [45]. This implies that R. microplus is one of the most competent vectors for A. marginale. Comparatively, A. marginale was highly detected in Hy. dromedarii, followed by R. microplus in the present study. This pathogen was also detected in four other tick species, including Hy. anatolicum, Hy. dromedarii, R. sanguineus and Hy. scupense. To the best of our knowledge, the detection of A. marginale in R. sanguineus is exceptionally rare [46], and this pathogen has not been detected in Hy. scupense. However, experimentally it has been demonstrated that this pathogen can be successfully transmitted by R. sanguineus [47,48]. Therefore, such unexpected outcomes need to be further evaluated because the presence of a pathogen DNA in a tick species does not ensure it as a biological vector. Moreover, a global increase in moments of infected/carrier livestock and/or tick-infested livestock across international borders can further worsen the situation regarding this pathogen [49].
For the host range of ticks, the resemblance among hosts' ecology might be more significant than evolutionary similarity [50]. A wide host range was recorded for Hy. anatolicum that could be attributed to its two or three host life cycle with the infestation on different ungulates [51]. A comparatively wide host range was also noted for one host tick species such as R. microplus. This might be due to common practices in the study area, such as placing different hosts in the same shelter, overcrowded herds, and combined grazing.
Research has been focused on understanding the evolutionary history and taxonomy of ticks and TBPs using standard genetic markers [36,52,53]. The mitochondrial gene cox1 has been considered an appropriate genetic marker for understanding tick phylogenetic relationships, especially at the species level [52]. The 16S rRNA gene has also been considered a reliable marker for tick identification [36,52] and is of prime importance in evaluating bacterial phylogeny and taxonomy [54,55]. When taking these into account, the cox1 sequences were obtained for eight tick species (Hy. dromedarii, Hy. anatolicum, Hy. scupense, Hy. isaaci, Ha. bispinosa, R. microplus, R. sanguineus, and R. haemaphysaloides). For the remaining four tick species (R. turanicus, Ha. cornupunctata, Ha. montgomeryi, and Ha. sulcata), we were able to obtain only 16S rRNA sequences. The A. marginale associated with these ticks was molecularly assessed by targeting the partial 16S rRNA gene. Except for Ha. cornupunctata, all other Haemaphysalis species were clustered with related species reported from the Oriental and neighboring Palearctic zoogeographical regions (Figures 5 and 6). In the cox1-based phylogenetic tree, the monophyletic clade containing Ha. bispinosa was basal to the remaining ixodid tick species. In tick 16S rRNA-based tree,