Risk factors associated with tick infestations on equids in Khyber Pakhtunkhwa, Pakistan, with notes on Rickettsia massiliae detection

Studies on ticks infesting equids are lacking in various parts of the world, including Khyber Pakhtunkhwa (KP), Pakistan. The aim of this study was to investigate the diversity of ticks infesting equids, associated risk factors and rickettsial detection in ticks from equids in KP. Inspection of 404 equid hosts from November 2018 to October 2019 resulted in the collection of 550 ticks. Data on tick-associated risk factors were collected from equid owners by means of a questionnaire. After morphological identification, partial DNA sequences of the tick mitochondrial 16S rRNA gene were used for taxonomic confirmation of species. Partial sequences of the gltA and ompA genes were used for Rickettsia detection in ticks. A total of 550 tick specimens were collected on 324 (80.2%) of the equids inspected, of which 161 were horses (50%), 145 (45%) were donkeys and 18 were mules (5%). The ticks were identified as belonging to the following five species: Rhipicephalus microplus (341 specimens, 62% of the total ticks), Rh. haemaphysaloides (126, 23%), Rh. turanicus (39, 7%), Rh. sanguineus (s.l.) (33, 6%) and Hyalomma anatolicum (11, 2%). The most prevalent tick life stage was adult females (279, 51%) followed by adult males (186, 34%) and nymphs (85, 15%). Higher tick infestations were observed on male equids (relative risk [RR] 0.7432, P < 0.0005) and adult equids (RR 1.268, P < 0.0020). Ticks were frequently attached to the axial region of horses (55, 21%), sternum of donkeys (44, 21%) and belly of mules (19, 23%) (P < 0.04). Temporal patterns of tick infestation in association with temperature and humidity were highly significant (P < 0.05). Risk factors, such as animal housing (P < 0.0003), living management (P < 0.006), grazing type (P < 0.01) and location in hilly areas (P < 0.02), significantly enhanced the chances for tick infestation. Tick species analyzed in this study were phylogenetically related to species from Afghanistan, China, South Africa and Taiwan. Partial sequences of the gltA and ompA genes obtained from Rh. microplus and Rh. haemaphysaloides were 100% identical to the spotted fever group pathogen Rickettsia massiliae. Equids exposed to significant risk factors were infected by one or more of at least five tick species in KP, Pakistan, and some of the ticks harbored the human pathogen R. massiliae.

regions, where approximately 80% of the world's cattle population is at risk of infestation [1]. These hematophagous ectoparasites play a major role in the transmission of pathogens, including bacteria, such as intracellular coccobacilli of the genus Rickettsia, and several protozoans and viruses that cause disease and are a threat to human and veterinary health [2][3][4]. Only 10% of tick species have been identified as carriers of pathogens; the remaining species require further research to determine whether they are vectors of disease.
Both the persistence of various life stages of ticks and their distribution depend mainly on such factors as habitat, climatic conditions, host availability and anthropogenic activities. Both biotic (vegetation structure, host availability and management) and abiotic (rainfall regimen, temperature range) factors are most likely to influence tick distributions on a global scale [5][6][7]. The potential impacts of climate change and transboundary movements of humans and animals on the distribution of ticks have been widely studied [8,9]. Understanding the niche of a tick species is a complex exercise due to the multiple variables operating on the different stages of the tick life-cycle, including regulation of the patterns of tick and host abundance and their spatial encounters [10]. Some tick species exhibit ecological flexibility and adapt easily to changing climate and novel habitats [10,11]. Efforts have been made at global, regional and local scales to characterize the major environmental factors influencing tick distributions using statistical models and correlative approaches [12]. The documentation of risk factors associated with tick distribution contributes to the design of cost-effective control measures.
Livestock activities are an integral part of Pakistan's economy and the backbone of rural income, as more than 70% of the population live in rural areas [13]. According to the Agricultural Census Organization, the total equid population in Pakistan was 4.8 million in 2006, which rose to 5.5 million as per the census report in 2012-2013. The Pakistan Economic Survey 2018-2019 reported that the horse population was 0.4 million, donkey population was 5.4 million and mule population was 0.2 million (https:// www. thene ws. com. pk/ print/ 482965-donke ys-popul ation-incre ased). Khyber Pakhtunkhwa (KP) is a mountainous region of Pakistan that lacks sufficient access to transport infrastructure. Consequently, equids play important roles in providing transportation, but also in farming, pumping water and milling activities in elevated areas. Equids are also used in fairs and games, such as horse races and polo festivals. However, despite their economic importance, the welfare of equids has received far less attention than that given to other animals in Pakistan.
Equids are frequently employed in outdoor activities and are therefore exposed to tick infestations [14,15]. Previous studies conducted in Pakistan were mostly confined to domestic animals but excluded equids, thereby limiting our understanding of the prevalence of and associated risk factors for tick infestations on these animals [13,16,17]. In addition, although Pakistan's climatic conditions are suitable for the survival of tick species that parasitize domestic animals, investigations on the frequency and distribution of ticks infesting equids are lacking, as are data on Rickettsia associated to equid ticks that are pathogenic to humans. The present study was designed to investigate the diversity of ticks, assess the presence of Rickettsia and estimate the risk factors associated with ticks infesting equids in KP.

Study area
We selected nine districts in northern, southern and central KP ( Based on access, nine regions of each of these districts were selected for tick sampling. A global positioning system (GPS) was used to obtain the geographic coordinates data and loaded onto a Microsoft Excel (Microsoft Corp., Redmond, WA, USA) spreadsheet to develop a distribution map for the study area using ArcGIS 10.3.1 (Fig. 1).

Ethical consent and information collection
Ethical consent was obtained from the advanced studies and research board of the Abdul Wali Khan University Mardan. Oral/written consent was obtained from the owner of equids for the collection of ticks. A questionnaire was provided to equid holders that assisted in collection of data such as: host-related information (gender and age), eco-geography and herd management.

Tick collection and morphological identification
Ticks were collected between November 2018 and October 2019 from horses, donkeys and mules at different sites (villages, towns) in tehsils (i. e. administrative divisions) within the nine districts (Fig. 1). Two visits per month were made to each area, and information related to climate was retrieved. In infested areas, ticks of all stages were removed from the base of the tail, ears, perianal area, sternum, scrotal and belly, using forceps without damaging the morphological features of the ticks. The collected ticks were placed in vials containing 100% ethanol and brought to the laboratory for identification.
Collected ticks were identified morphologically under a StereoZoom microscope (HT StereoZoom) and following taxonomic keys [18,19]. The specimens were preserved in 100% ethanol for later molecular analysis.

DNA extraction and PCR
Tick specimens were rinsed in distilled water and 70% ethanol and dried prior to DNA extraction [17]. Specimens were cut into pieces with a sterilized scalpel inside Eppendorf tubes. Genomic DNA of 91 randomly selected ticks was extracted using the phenol-chloroform method [20,21]. Sterile conditions were maintained throughout the process. The concentration of DNA in each sample was quantified using a spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). The extracted DNA was stored at − 20 °C.
Conventional PCR targeting a fragment (approx. 460 bp) of the tick mitochondrial 16S rRNA gene was carried out as previously described [22] to verify successful extractions (Table 1). A real-time PCR (StepOne; Applied Biosystems, Thermo Fisher Scientific, Waltham, MA USA) was performed to detect a short fragment (147 bp) of the Rickettsia citrate synthase-encoding gene (gltA) using primers CS-5 and CS-6 and an internal fluorogenic probe (6-carboxyfluorescein [6-FAM]) [2]. Positive samples were subjected to two conventional PCR assays, one to obtain a larger sequence (401 bp) of the rickettsial gltA gene, using primers CS-78 and CS-323 [2], and the other to amplify a 532-bp fragment of the rickettsial outer membrane protein A (ompA) gene [23]. Conventional PCR reactions consisted of a 25-µl mix (2 µl genomic DNA, 1 µl of each primer, 8.5 µl PCR water and 12.5 µl master mix; Thermo Fisher Scientific). PCR cycling parameters were as previously reported [2,22,24] (see Table 1). PCR products were electrophoresed in 1.5% agarose gels, and the results were captured with GelDoc (UVP BioDoc-It imaging system; Analytik Jena AG, Jena, Germany). Amplicons of the expected size were purified with ExoSAP-IT (Thermo Fisher Scientific) and sequenced using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Thermo Fisher Scientific) in an ABI 3500 Genetic Analyzer automated device (Applied Biosystems, Thermo Fisher Scientific).

Sequence and phylogenetic analyses
The obtained sequences were assembled and trimmed in SeqMan v 5.00 (DNASTAR, Madison, WI, USA) and subjected to BLASTn analyses (www. ncbi. nlm. nih. gov/ blast) to infer closest identities with the organisms available in GenBank [25]. Alignments with our sequences and GenBank-retrieved sequences were constructed with CLUSTAL W in BioEdit [26,27] manually edited with GeneDoc [28]. A phylogenetic tree for the tick mitochondrial 16S rRNA gene was implemented in MEGA X [29] and constructed using the neighbor-joining method with 1000 bootstrap replicons [30]. For alignment of the ompA rickettsial gene, we constructed a phylogenetic tree in MrBayes using Bayesian statistics [31], with four independent Markov chain runs for 1,000,000 metropolis-coupled Markov chain Monte Carlo (MCMC) generations, sampling a tree every 100th generation. The first 25% of the trees represented burn-in, and the remaining trees were used to calculate the Bayesian posterior probability.

Statistical analysis
The recorded observations were assembled and arranged in the spreadsheets of Microsoft Excel V 2016. The chisquare test (χ 2 ) was used for differences, and relative risk (RR) was analyzed using GraphPad Prism v. 5.00 (Graph-Pad Software Inc., San Diego, CA, USA), with the 95% confidence interval (CI). Significance was set at P < 0.05.

Tick infestation on different body parts of the equid host
Different parts of the body of the equid hosts were examined for ticks. The highest infestations were observed on the axial region of horses (63, 24%), the sternum of donkeys (44, 21%) and the belly of mules (19, 23%). The lowest tick-infested body parts were the dewlap in horses (37, 14%), the axial region in donkeys (25, 12%) and the scrotum of mules (10, 12%). Statistical analysis of infestation of the body regions of the equid hosts revealed significant differences (χ 2 = 18.66, 10; P < 0.045) ( Table 4).

Assessment of risk factors associated with tick infestations
Several independent variables (Table 5) potentially associated with tick infestations on equids were analyzed and the risk factors calculated. Analysis of infestation according to sex revealed that there were more tick infestations on male equids than on female equids, with the difference being statistically significant in the tick acquisition. When age was considered, the highest tick infestation was found on adult hosts aged 4-6 years compared to hosts aged 1-3 years. Most of the equids examined were in the plain areas, but tick infestations were predominantly found on equids in hilly areas, and the association between a high risk of tick infestation and equids in hilly areas was significant. When animal housing was considered, the tick infestation rate was found to be higher on equids in muddy systems than those maintained in concrete and semi-concrete housing. The risk of tick infestation was also higher in those managed in mud housing than in those in open grazing systems. Equids raised in a herd were found to have a higher risk of tick infestation than those reared solely. Providing fresh or stored food to the equids had no impact on the tick infestations ( Table 5).    H. anatolicum (MW172215). The phylogenetic analysis of the collected five species revealed an evolutionary relatedness with homologous species reported in Pakistan's neighboring countries. In a phylogenetic tree, the obtained sequences of prevalent tick species clustered with the same species reported in the bordering countries (Fig. 3). JN043507, CP000683, U59719). From these same ticks, a consensus partial sequence (592 bp) of the highly polymorphic Rickettsia ompA gene was 100% identical to a sequence of R. massiliae from GenBank (CP003319). The obtained sequences of gltA and ompA were uploaded to GenBank with accessions numbers MW922872 and MW928497 respectively. The obtained sequence of the ompA gene for Rickettsia spp. clustered with sequences of R. massiliae in a phylogenetic tree (Fig. 4).

Discussion
Ecological changes and an increasingly evident increase in global temperature have widely impacted the distribution of ticks and tick-borne diseases, ultimately increasing their threat to public and veterinary health [9,10]. Based on spatial distribution and molecular evidence, it has been suggested that global climate change has played a significant role in the spread of tick species [10]. Previously, we reported numerous tick species infesting diverse hosts, including equids (horse), in selected districts of KP, Pakistan [13,17]. However, a knowledge gap remained on the identity and effect of various factors affecting ticks infesting equid hosts and their associated Rickettsia spp. To fill this gap, we designed this study to investigate equid hosts for tick infestation, associated risk factors and the molecular phylogeny of these ticks in selected districts of KP, Pakistan. The collected ticks from five selected districts were categorized into two genera which included five medically important tick species: Rh. microplus, Rh. haemaphysaloides, Rh. turanicus, Rh. sanguineus (s.l.) and H. anatolicum. These ticks were found to be infesting equid hosts (horses, mules,  (27) 97.50 (16.26) donkeys) throughout the year in extreme cool to extreme hot weather conditions. The bacterium Rickettsia massiliae was detected in all Rh. microplus ticks collected in all nine study districts and in Rh. haemaphysaloides collected in Bajaur, Mardan and Charsadda districts. Ecological conditions and biotopes drive the geographic distribution and the risks associated with the possible harmful effects of tick species. The favorable environmental conditions in the region extend tick adaptation and also result in infestation of novel hosts [13,16,17]. Risk factors, such as age, altitude, housing type, grazing, living management and food conditions, of equid hosts with respect to tick infestation were statistically significant, with the exception of food conditions. A strong correlation was observed between the intensity of tick infestation and temperature and relative humidity, with the highest infestation recorded in the summer (June-August) and lowest infestation recorded in the winter (December-January); these results agree with previous reports in the region [16,17,32]. The variable climatic conditions of the KP favor tick survival and reproduction, resulting in high infestations on various hosts [17,33].
The younger age group (1-years) of equid hosts was found to have lower tick infestations, possibly correlated with a strong immune response, shorter exposure to questing ticks, shorter time in free grazing conditions and frequent grooming, all of which protect younger equids from tick infestation [34][35][36]. Higher tick infestations on adult hosts may be due to the larger surface area of the body, which enhances the chances  of tick attachment, frequent use and ultimate resistance to acaricide and low immunity [32,[36][37][38]. The muddykeeping system of hosts provided favorable conditions for the survival and reproduction of ticks, as compared to concrete and semi-concrete keeping systems. Host density in herds ultimately increases the chances of tick infestation, compared to animals kept alone, due to the availability of host-questing ticks detached in the herd area [39][40][41]. Pastures with long grasses in the hilly areas may provide favorable conditions for the questing ticks where they await and attach to open-grazing equids.
Ticks collected across the different study regions belonged to two genera and five species. The most prevalent tick was Rh. microplus, which has been reported to the dominant tick species in the region [13,16,17]. The highest numbers of ticks were collected in the months of June to August, while the lowest numbers were collected in the coldest months of November to January, which is in accordance with the results of a previous study [17]. It is important to mention that several Haemaphysalis and Hyalomma species other than H. anatolicum are also abundant tick species in Pakistan [13,16,17], but these latter species were not found to infest the equids inspected in this study throughout the study period, possibly due to the availability of other suitable and accessible hosts for blood feeding. Several markers have been utilized to delineate the accurate identity of a tick species at the molecular level. For example, 16S rRNA, COX1, 12S rRNA and ITS are well-known genetic markers that have been implied to separate closely related tick species [17,[42][43][44]51]. Since the objective of this work was only to identify tick species-and not population structurewe used 16S rRNA to identify the collected tick species and to establish their evolutionary relationship [43][44][45][46]. The phylogenetic analysis of the collected ticks showed the closest relation with the same species reported from other countries, including Afghanistan, China, South Africa and Taiwan.
Tick-borne Rickettsia spp. have been collected from several Rhipicephalus ticks infesting different hosts in various regions of the world [13,[47][48][49][50][51]. In Pakistan, reports are not available on ticks infesting equids associated with Rickettsia infection. In this study, Rickettsia massiliae was detected in Rh. microplus and Rh. haemaphysaloides using gltA and ompA rickettsial markers. There have been several reports of spotted fever due to R. massiliae in Europe [49]. Noteworthy, this is the first report of R. massiliae associated with equid ticks in Pakistan. Our findings highlight that human illness due to this spotted fever group of pathogens could be underreported in Pakistan, where to date tick-borne rickettsioses have not been reported.

Conclusions
To our knowledge this is the first report of equid ticks from selected regions of KP. The risk factor analyses showed that temperature, humidity, age, keeping and grazing systems, sex and host species considerably affected the intensity of tick infestations. The phylogenetic tree of the collected ticks based on the tick mitochondrial 16S rRNA gene (2 genera and 5 species) revealed a close resemblance to tick species reported from Afghanistan, China, South Africa and Taiwan. The human pathogen R. massiliae, a tick-borne pathogenic spotted fever group Rickettsia species, was detected in Rh. microplus and Rh. haemaphysaloides. These results demonstrate that it is important to investigate the vector potential of ticks for other pathogens infesting equids in different regions of Pakistan to avoid the emergence of zoonotic infections. They will also provide valuable data for use in developing integrated control management strategies against ticks and tick-borne diseases in Pakistan.