Next Article in Journal
Knowledge and Perception of Rabies among School Children in Rabies Endemic Areas of South Bhutan
Next Article in Special Issue
Spatio-Temporal Patterns of Dengue Incidence in Medan City, North Sumatera, Indonesia
Previous Article in Journal
The COVID-19 Pandemic: Disproportionate Thrombotic Tendency and Management Recommendations
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
TropicalMed
Volume 6
Issue 1
10.3390/tropicalmed6010027
tropicalmed-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Identification, Distribution, and Habitat Suitability Models of Ixodid Tick Species in Cattle in Eastern Bhutan

by
Jamyang Namgyal
1,2,*,
Tim J. Lysyk
2,
Isabelle Couloigner
2,3,
Sylvia Checkley
2,
Ratna B. Gurung
4,
Tenzin Tenzin
4,
Sithar Dorjee
5 and
Susan C. Cork
2
1
District Veterinary Hospital, Department of Livestock, Ministry of Agriculture and Forests, Trashigang 42001, Bhutan
2
Department of Ecosystem and Public Health, Faculty of Veterinary Medicine, University of Calgary, Calgary, AB T2N 1N4, Canada
3
Department of Geography, University of Calgary, Calgary, AB T2N 1N4, Canada
4
National Centre for Animal Health, Department of Livestock, Ministry of Agriculture and Forests, Thimphu 11001, Bhutan
5
Khesar Gyalpo University of Medical Sciences of Bhutan, Thimphu 11001, Bhutan
*
Author to whom correspondence should be addressed.
Trop. Med. Infect. Dis. 2021, 6(1), 27; https://0-doi-org.brum.beds.ac.uk/10.3390/tropicalmed6010027
Submission received: 28 January 2021 / Revised: 13 February 2021 / Accepted: 15 February 2021 / Published: 19 February 2021
(This article belongs to the Special Issue Spatial Epidemiology of Vector-Borne Diseases)
Browse Figures
Versions Notes

Abstract

:
Tick infestation is the most reported parasitological problem in cattle in Bhutan. In May and June 2019, we collected ticks from 240 cattle in two districts of Eastern Bhutan. Tick presence, diversity, and infestation prevalence were examined by morphological identification of 3600 live adult ticks. The relationships between cattle, geographic factors, and infestation prevalence were assessed using logistic regression analyses. Habitat suitability for the tick species identified was determined using MaxEnt. Four genera and six species of ticks were found. These were Rhipicephalus microplus (Canestrini) (70.2% (95% confidence interval (CI): 68.7–71.7)), Rhipicephalus haemaphysaloides Supino (18.8% (95% CI: 17.5–20.1)), Haemaphysalis bispinosa Neumann (8.2% (95% CI: 7.3–9.1)), Haemaphysalis spinigera Neumann (2.5% (95% CI: 2–3)), Amblyomma testudinarium Koch (0.19% (95% CI: 0.07–0.4)), and a single unidentified Ixodes sp. Logistic regression indicated that the variables associated with infestation were: longitude and cattle age for R. microplus; latitude for R. haemaphysaloides; and altitude and cattle breed for H. bispinosa and H. spinigera. MaxEnt models showed land cover to be an important predictor for the occurrence of all tick species examined. These findings provide information that can be used to initiate and plan enhanced tick surveillance and subsequent prevention and control programs for ticks and tick-borne diseases in cattle in Bhutan.

1. Introduction

Bhutan is a small kingdom located in the Eastern Himalayas between latitudes 26°45′ N and 28°10′ N, and longitudes 88°45′ E and 92°10′ E. It shares borders with China to the north, the Indian states of Assam and West Bengal to the south, Arunachal Pradesh to the east, and Sikkim to the west. It has a total land area of 38,394 km2, out of which 70.8% is covered by forest mainly dominated by broadleaf and mixed conifers [1]. It is divided administratively into 20 districts (dzongkhags) and 205 subdistricts (gewogs). The population is 735,553 (as of 30 May 2017), of which 62.2% lives in rural areas, and their livelihoods depend on agriculture and livestock farming [2]. Geographically, the country is characterized by high mountains, dense forests, and fast-flowing rivers that form narrow valleys before flowing out onto the vast north Indian plains. There is a pronounced south–north elevation gradient (100–7500 m above sea level (masl)) and an inverse north–south precipitation gradient (500 to >2000 mm) [1,3]. This extreme variation in altitude and the impact of the North Indian monsoon have resulted in extremely diverse climatic conditions and ecosystems across the country, resulting in six agroecological zones (wet subtropical, humid subtropical, dry subtropical, warm temperate, cool temperate, and alpine) [4]. Vegetation coverage, cropping system, land use, and livestock farming are mainly determined by these agroecological zones.
Livestock farming plays an important role in the rural economy of Bhutan. It not only provides the major income for 49.1% of the population who are subsistence farmers [2] but also contributes toward meeting the increasing demands for dairy products in the country. Cattle, including yaks (Bos grunniens L.) and mithuns (Bos frontalis Lambert), are the predominant livestock farmed with the population of 0.35 million [5]. They are the most important livestock farmed with their diverse role in providing milk, draught power, transport, and organic manure [6]. The crossbreeds of local indigenous Siri cattle (Bos taurus indicus L.) and mithun comprise 56% of the total cattle population while European breeds (Bos taurus taurus L.) such as Jersey, Brown Swiss, and Holstein Friesian form 30%. The yaks and their crossbreeds comprise 14% of the country’s cattle population [5].
Tick infestation is a significant challenge for the livestock sector in Bhutan due to direct and indirect impacts on cattle health and production [7]. In 2019, 89% of the parasite related cases reported in cattle were due to tick infestation, and 42% of the cattle in the country were reported to have been treated for it, costing the government approximately 3.18 million Bhutanese Ngultrum (1 USD = Nu.70) for purchasing acaricides alone [8]. Tick-borne diseases like anaplasmosis, babesiosis, and theileriosis are also present in cattle in Bhutan, especially in the southern subtropical districts [7,8]. Despite ticks being prevalent throughout the country, there is currently limited data available on tick diversity, infestation prevalence, tick ecology, and geographic distribution. Current knowledge consists of one unpublished study conducted over a period of one year in cattle in Eastern Bhutan (Dr. Susan C. Cork, personal communication, 2018) and one government report [9] from Western Bhutan. The unpublished study may be the first to identify some tick species present on cattle in Bhutan. They reported Rhipicephalus microplus (Canestrini) as the most predominant tick species with Haemaphysalis spp. and Hyalomma spp. identified to the genus level. The other study by the Regional Livestock Development Center [9] in Western Bhutan identified the genera Rhipicephalus (Boophilus) spp., Rhipicephalus spp., Haemaphysalis spp., Ixodes spp., and Amblyomma spp. but no information on species is available.
In this study, we used targeted field sampling to: (1) determine the presence, diversity, and infestation prevalence of tick species in cattle in two districts of Eastern Bhutan; (2) examine the variation in infestation prevalence of tick species between agroecological zones in Trashigang and Pema Gatshel districts; (3) assess the relationships between cattle, geographic factors, and infestation prevalence using logistic regression analyses; and (4) develop habitat suitability models for tick species identified in Eastern Bhutan, under current environmental conditions, using MaxEnt [10]. The findings from this study are expected to provide information necessary to initiate and plan a targeted tick control program and to develop enhanced surveillance for ticks and tick-borne diseases in cattle in Bhutan.

2. Materials and Methods

2.1. Study Areas

The study was conducted in the districts of Trashigang and Pema Gatshel in Eastern Bhutan (Figure 1), covering an area spanning latitudes 26°75′–27°5′ N, and longitudes 91°–91° 85′ E. The study area covered the entire range of agroecological zones represented in Bhutan. These two districts also contain diverse breeds of cattle ranging from the indigenous breed Jaba (Bos taurus indicus) common in the lower subtropical plains to the yaks found at the higher alpine areas. The elevation of the study area ranged from 46 to 4571 masl. Trashigang, located toward the north, is predominantly a temperate district characterized by warm summers and cold winters, and shares a border with the Indian state of Arunachal Pradesh in the east. Pema Gatshel, located at the warmer south, is mainly a subtropical district characterized by hot-humid summers and cool winters, and shares a border with the Indian state of Assam in the south. The geographic and climatic features along with the cattle populations of the two districts are summarized in Table 1. Most of the land in both districts is associated with high mountainous terrains separated by narrow valleys. Human settlements and farming activities are generally limited to these narrow valleys and the gentle slopes on the mountains.

2.2. Sample Size

The number of sites sampled was determined by the administrative units, i.e., subdistricts/gewogs. Trashigang and Pema Gatshel have 15 and 11 subdistricts, respectively. All subdistricts in the two districts were included for sampling except the two subdistricts in the alpine zone of Trashigang district where we were unable to obtain tick specimens during the field survey conducted in May and June 2019. As the primary objective of the study was the detection of tick species presence, the sample size for each subdistrict was calculated using the formula n = [ 1 − ( 1 − p l ) 1 d ] ( N − d 2 ) + 1 [11], where n = required sample size, N = number of cattle in subdistricts, d = minimum number of tick infested cattle expected in the population, and pl = probability (0.95) of finding at least one infested cattle in the sample. The number of cattle (N) for each subdistrict was obtained from the annual livestock census data [5]. The expected proportion (d) was calculated as 0.3 for all subdistricts based on the tick infestation treatment records from the veterinary information system and the livestock population data of each subdistrict maintained by the Department of Livestock in Bhutan [5]. The sample size calculated for each subdistrict was 10 animals. Since the study area had 24 subdistricts to be sampled, the overall sample size was 240 animals.

2.3. Tick Sampling and Specimen Identification

For selecting the required number of animals, the list of households owning cattle in each subdistrict was obtained from subdistrict livestock offices and used as the sampling frame. Ten households owning cattle were selected from each subdistrict using simple random sampling in MS Excel 2016 (Microsoft Excel 2016, Redmond, WA, USA). Overall, this study covered 240 households and collected ticks from 240 cattle. Tick collection was conducted once per household in May–June 2019 in all 24 subdistricts. Prior to tick collection, the objectives of the study were explained to the owners, and verbal consent was obtained to collect ticks from their cattle. All selected households agreed to the sampling. In each farm linked to the selected household, one infested animal that could be properly restrained was selected for sampling (convenience sampling), and 15 ticks were randomly collected from different predilection sites. The ticks collected from each animal were placed separately in Eppendorf tubes containing 70% ethanol and labeled with a unique sample ID that included district and subdistrict codes and the animal number. Other associated information such as age, sex and breed of the animal, GPS coordinates of the location, site altitude, owners’ names, and dates of collection were recorded in an excel spreadsheet for data analyses and future reference. The samples were then transported to the veterinary laboratory at the National Centre for Animal Health (NCAH), Thimphu, Bhutan for morphological identification.
Ticks were identified under a stereo microscope using a two-step process of species identification. Ticks were first identified to genus using the keys described in Matthysee and Colbo [12] and Walker et al. [13]. Ticks were then identified to species using Anastos [14] and Robinson et al. [15] for the members of the genus Amblyomma; Matthysee and Colbo [12] and Estrada-Pena et al. [16] for the members of the genus Rhipicephalus (Boophilus); Walker et al. [13] and Anastos [14] for the remaining Rhipicephalus; Geevarghese and Mishra [17] and Nutall et al. [18] for the members of the genus Haemaphysalis; and Guo et al. [19] for the members of the genus Ixodes. The tick specimens, including the voucher specimens, are preserved at the veterinary laboratory, NCAH, Thimphu, Bhutan. Macro photography was done on some of the selected voucher specimens to get representative pictures of the identified ticks.

2.4. Environmental Data

Bioclimatic variables were downloaded from the WorldClim database (Version 2, Food and Agriculture Organization, Quebec City, Quebec, Canada) (https://www.worldclim.org/ (accessed on 28 January 2021)) at 30 s (1 km2) spatial resolution [20]. The digital elevation–shuttle radar topography information (DEM_SRTM) at 1 arc-second global was downloaded from the United States Geological Survey (USGS) database (https://www.usgs.gov/ (accessed on 28 January 2021)). The most recent land use and land cover data (LULC 2016) for Bhutan at 30 min spatial resolution was obtained from the National Land Commission (NLC) of Bhutan (https://www.nlcs.gov.bt/ (accessed on 28 January 2021)). LULC 2016 uses 17 classes: alpine scrubs (1); broadleaf (2); built-up (3); chuzhing (4); chirpine (5); fir (6); kamzhing (7); lake (8); landslides (9); meadows (10); mixed conifer (11); non built-up (12); orchards (13); rivers (14); rocky outcrops (15); shrubs (16); and snow and glaciers (17). Chuzhing refers to irrigated and bench terraced agriculture land for paddy cultivation, and kamzhing represents cultivated rain-fed dry land primarily used for cereal cropping [21]. DEM_SRTM and LULC 2016 were organized as raster (grid) type files to the spatial extent and resolution of the WorldClim layers. Data preparation was done using raster [22], rgeos [23], and rgdal [24] packages in R (R Core Team, Vienna, Austria) [25].
Seven bioclimatic variables (i.e., Bio 2, Bio 3, Bio 4, Bio 7, Bio 13, Bio 14, and Bio 15) were excluded, as they were not deemed ecologically relevant to the study area. Bio 2 (annual mean diurnal range) was excluded due to limited data range (2.3 °C) in the study area. Bio 3 (isothermality) was excluded, as there is little day to night temperature oscillation in the study area. Bio 4 (temperature seasonality) and Bio 7 (annual temperature range) was excluded, as their information is included in other temperature variables. Bio 13 (precipitation of the wettest month), Bio 14 (precipitation of the driest month), and Bio 15 (monthly precipitation variation) were excluded, as the precipitation in Bhutan is seasonal depending on the North Indian monsoon. The seasonal variations are more important than the monthly variations; therefore, the precipitation of the quarters was used. The environmental variables used in the MaxEnt are listed in Table 2.

2.5. Statistical Analyses

The raw data were collated in Excel spreadsheets and imported into R (R Core Team, Vienna, Austria) [25] for analyses. The tick species considered for statistical analyses were R. microplus, R. haemaphysaloides, H. bispinosa, and H. spinigera. The remaining two species, A. testudinarium, and Ixodes sp. were excluded from analyses, as they were infrequently collected. Descriptive analysis was performed with the entire dataset to calculate proportions and frequencies. Tick infestation prevalence was calculated as the number of cattle infested with ticks divided by the total number of cattle examined among the households sampled [26]. The 95% binomial confidence interval for the infestation prevalence was calculated using the Clopper–Pearson exact method using the “PropCI” package [27]. Pearson’s chi-squared test using the native “stats” package [25] was performed to compare the difference in proportions of tick infestation prevalence between the districts.
Logistic regression analyses were conducted using the cattle age (categorized as young or adult), sex, breed (categorized as European or indigenous), site altitude (in 100 m), latitude (in 1/10 decimal degrees), and longitude (in 1/10 decimal degrees) as the explanatory variables against the binary outcome variables of infestation prevalence of each tick species (categorized as infested or not). Correlation among the explanatory variables was assessed using the “Hmisc” package [28]. The explanatory variables with p-value ≤ 0.25 in univariable analyses were selected for multiple logistic regression analysis. A forward stepwise approach was used to manually build the final multiple logistic regression models. First, a variable with the smallest p-value in the univariable analysis was included into the model. Then, each of the remaining variables was individually added to the model to assess whether its addition improved the fit of the model significantly at p-value ≤ 0.05. A likelihood ratio test was used to select the variable that had the greatest improvement in the likelihood ratio statistic, and the process was repeated. Only the variables at the p-value ≤ 0.05 level of significance were retained in the final model. Multicollinearity of the explanatory variables in the models was assessed using the variance inflation factor (VIF) at the cut-off of 2.5 [29]. Interactions were assessed by adding a cross-product term (i.e., Latitude × Longitude). The odds ratio (OR) and its 95% confidence interval (CI) of the variables associated with the outcome variables were calculated from the final models. The goodness-of-fit for the final models were evaluated using the “LogisticDx” package [30]. The residual analysis of the final models was done using the “car” package [31].

2.6. Habitat Suitability Modeling

The habitat suitability modeling was conducted using MaxEnt v3.4.1 (New York, NY, USA) [10] in R (R Core Team, Vienna, Austria) [25] using the “dismo” package [32] and the “prepPara” function [33]. The presence points (i.e., latitude and longitude pair) for four tick species, R. microplus (204 points), R. haemaphysaloides (91 points), H. bispinosa (72 points), and H. spinigera (28 points) were used for the modeling. The presence data were transformed into the spatial data frame using the coordinates function in the “sp” package [34]. All environmental layers were stacked as raster layers and restricted to the study area using crop and mask functions in the “raster” package [22].
In the first step, MaxEnt was run with iterations and background points set to default using all the environmental variables that were transformed into the hinge feature. Variables contributing less than 1% to the increasing training gain or less than 1% permutation importance were considered non-significant and excluded from further analyses [35]. Variables selected from the first MaxEnt run for each tick species were examined for correlation using the corr.test function in the “psych” package [36].
In the second step, for each tick species, competing models were built using all selected variables and then subsequently reducing the model by removing the less important variables based on Akaike information criteria (AIC) [37]. AIC was calculated using the cal.aicc function from the “ENMeval” package [38], and model evaluations were done using the evaluate function from the “dismo” package [32]. Competing models for each tick species were then compared using threshold independent receiver operating characteristics (ROC) graphs [39], correlation, and the AIC [37,40]. The best model for each tick species were then selected based on AIC values, and whenever AIC values were indecisive, the most parsimonious models were selected.
Using “rasterVis” [41] and “RColorBrewer” [42] packages, the habitat suitability of the best model for each tick species was reclassified and mapped following Zuliani et al. [35] into 5 classes: class 1 (very low) for the probability level (p) between 0 and 0.2; class 2 (low) for p between 0.2 and 0.4; class 3 (moderate) for p between 0.4 and 0.6; class 4 (high) for p between 0.6 and 0.8; and class 5 (very high) for p between 0.8 and 1.0. For the best models, the response curves of the predictors were qualitatively assessed to understand the relationship between the variation in individual variables and the probability of tick species occurrence.

3. Results

3.1. Tick Diversity and Infestation Prevalence

A total of 3600 ticks were collected and identified to four genera and five species (Table 3 and Figure 2, Figure 3 and Figure 4). The majority were Rhipicephalus microplus (Canestrini) (70.3%) followed by Rhipicephalus haemaphysaloides Supino (18.8%), Haemaphysalis bispinosa Neumann (8.2%), and Haemaphysalis spinigera Neumann (2.5%). The remainder were seven Amblyomma testudinarium Koch and one unidentified species of Ixodes sp. Amblyomma testudinarium and Ixodes sp. were found only in Pema Gatshel and Trashigang district, respectively. The variation in proportions of tick infestation prevalence between the two districts is given in Table 4.

3.2. Factors Associated with Tick Infestation in Cattle

The characteristics of the cattle sampled and geographic variables are shown in Table 5. For R. microplus, simple logistic regression (Table S1) indicated that infestation varied between cattle age groups, but not between sexes or breeds. The infestation declined with each of increasing altitude, latitude, and longitude in simple logistic regression (Table S1). The odds of cattle being infested with R. microplus were 5.6 times more likely in the young animals than the adult animals (OR = 5.6 (95%CI: 1.5–35.5)), and there was a 40% (OR = 0.6 (95%CI: 0.42–0.75)) lower likelihood of R. microplus infestation with every 1/10th degree increase in longitude, when the other variable in the model was held constant, respectively (Table 6).
For R. haemaphysaloides, simple logistic regression indicated that infestation varied among cattle age groups, but neither sex nor breed (Table S1). It also indicated that infestation increased with each of altitude, latitude, and longitude (Table S1). However, the final multiple logistic regression showed that the odds of cattle being infested with R. haemaphysaloides increased 2.0 times with every 1/10th degree increase in latitude (OR = 2.0 (95%CI: 1.7–2.5)) (Table 6).
For both Haemaphysalis species (i.e., H. bispinosa and H. spinigera), simple logistic regression indicated no significant effect of cattle age and sex on infestation, but a significant effect of the breed (Table S1). Both species showed no relationship with latitude or longitude, but a significant effect of altitude (Table S1). However, H. bispinosa indicated a significant positive relationship with altitude while H. spinigera indicated a significant negative relationship with altitude (Table S1). In the multiple logistic regression analysis, infestation of both Haemaphysalis species varied with both altitude and cattle breed. The odds of cattle being infested with H. bispinosa were 1.8 times more likely in the indigenous breeds than that of the European breeds (OR = 1.8 (95%CI: 1.0–3.3)), and the infestation increased 1.1 times with every 100 m increase in altitude (OR = 1.1 (95%CI: 1.0–1.1)) when the other variable in the model was held constant, respectively (Table 6). For H. spinigera, the odds of cattle being infested were 2.7 times more in the indigenous breeds than that of the European breeds (OR = 2.7 (95%CI: 1.2–6.4)), and there was 20% (OR = 0.8 (95%CI: 0.7–0.9)) lower likelihood of H. spinigera infestation with every 100 m increase in altitude, when the other variable in the model was held constant, respectively (Table 6).

3.3. Habitat Suitability Modeling

For R. microplus, three variables (i.e., DEM_SRTM, LULC, and Bio 18) achieved more than 1% contribution and permutation importance in the first MaxEnt run (Table S2). There was a high correlation (rs = −0.93) between DEM_SRTM and Bio 18 (Table S3). Four competing models were built with the selected variables (Table S4). The best model showed that elevation was the most important predictor at 56.1%, followed by land cover at 43.9% (Table 7). The analysis of the response curve showed that the probability of R. microplus occurrence increased with elevation between 500 and 1000 masl, then declined (Figure S1). As with land cover, the probability of occurrence was most significant in land cover classes classified as kamzhing (0.84) and meadows (0.81) (Figure S1). Shrubs were associated with a 0.55 probability of occurrence. Overall, 49.7% of the study area was predicted as moderate to high suitable areas for R. microplus (Table S5). Based on the habitat suitability map (Figure 5), the northeastern part of the study area with elevation above 2000 masl and the southernmost part with elevation less than 500 masl were predicted as very low suitable areas for R. microplus occurrence.
For R. haemaphysaloides, six variables (i.e., LULC, Bio 18, Bio 16, Bio 10, DEM_SRTM, and Bio 8) achieved more than 1% contribution and permutation importance in the first MaxEnt run (Table S2). Correlation analyses indicated that all variables were highly correlated (rs > 0.90) (Table S3). Eighteen competing models were built using the selected variables (Table S4). The best model showed that land cover was the most important predictor at 46%, followed by Bio 16 (precipitation of the wettest quarter) and Bio 10 (temperature of the warmest quarter) at 43.1% and 10.9%, respectively (Table 7). The analysis of the response curve showed that the probability of occurrence was most significant in land cover classes classified as kamzhing (0.82), followed by shrubs (0.6) and chuzhing (0.59) (Figure S2). The probability of occurrence increased with the temperature of the warmest quarter (Bio 10) between 16 and 25 °C, beyond which there was a decline (Figure S2). The probability of occurrence also increased with precipitation of the wettest quarter (Bio 16) between 400 and 1100 mm, then declined when Bio 16 exceeded 1200 mm (Figure S2). Overall, 27.4% of the study area was predicted as moderate to high suitable areas for R. haemaphysaloides (Table S5). Based on the habitat suitability map (Figure 5), the northeastern and the southernmost part of the study area were predicted as low suitable areas for R. haemaphysaloides occurrence.
For H. bispinosa, five variables (i.e., LULC, Bio 18, Bio 16, Bio 12, and Bio 11) achieved more than 1% contribution and permutation importance in the first MaxEnt run (Table S2). Correlation analyses indicated that all variables were highly correlated (rs > 0.90) (Table S3). Thirteen competing models were built using the selected variables (Table S4). The best model showed that land cover was the most important predictor at 63.2%, followed by Bio 18 (precipitation of the warmest quarter) and Bio 16 (precipitation of the wettest quarter) at 30.6% and 6.2%, respectively (Table 7). The response curve showed that the probability of occurrence was most significant in Kamzhing (0.83) and meadows (0.82) (Figure S3). Shrubs and Chuzhing had 0.64 and 0.55 probabilities, respectively. The probability of occurrence increased with precipitation of the warmest quarter (Bio 18) between 600 and 2000 mm, beyond which there was a decrease (Figure S3). The probability of occurrence also increased with precipitation of the wettest quarter (Bio 16) between 400 and 900 mm, then declined when Bio 16 exceeded 1600 mm (Figure S3). Overall, 24.9% of the study area was predicted as moderate to high suitable areas for H. bispinosa (Table S5). Based on the habitat suitability map (Figure 6), the northeastern and the southernmost part of the study area were predicted as low suitable areas for H. bispinosa occurrence.
For H. spinigera, seven variables (i.e., LULC, Bio 8, Bio 11, Bio 12, Bio 16, Bio 17, and Bio 19) achieved more than 1% contribution and permutation importance in the first MaxEnt run (Table S2). All variables were highly correlated (rs > 0.70) (Table S3). Twenty-six competing models were built using the selected variables (Table S4). The best model showed that Bio 16 (precipitation of the wettest quarter) was the most important variable at 45.7%, followed by Bio 19 (precipitation of the coldest quarter) and land cover at 33.4% and 20.4%, respectively (Table 7). The response curve showed that the probability of occurrence increased with precipitation of the wettest quarter (Bio 16) at 400 mm, then declined when Bio 16 exceeded 1600 mm (Figure S4). With the precipitation of the coldest quarter (Bio 19), the probability of occurrence increased between 15 and 40 mm, beyond which there was a decrease (Figure S4). As with land cover classes, the probability of occurrence was 0.88 in kamzhing and 0.67 with shrubs (Figure S4). Overall, 33.2% of the study area was predicted as moderate to high suitable areas for H. spinigera (Table S5). Based on the habitat suitability map (Figure 6), the northeastern and the southernmost part of the study area were predicted as low suitable areas for H. spinigera occurrence.

4. Discussion

The level of tick infestation in animals is generally influenced by both host and environmental factors. The host (cattle in our study) factors such as age, sex, and breed can influence the susceptibility of animals to tick infestation [43,44]. Ticks are also generally dependent on the temperature and rainfall for their development and activity [16]. In Bhutan, climatic variables such as temperature and rainfall are primarily determined by altitude and latitude [45]. Thus, in the first part of our study, the relationships between cattle and geographic factors and infestation prevalence were assessed. In the second part, we modeled the relationship between tick species presence and the environmental variables, using the MaxEnt [10] species distribution modeling to identify environmental factors associated with the geographical distribution of the tick species found in the study area. MaxEnt modeling was selected because it can build a reliable model of species distribution using presence-only data and environmental variables without assuming species absence in locations not sampled or surveyed [46].
In this study, R. microplus was found to be the most predominant tick species infesting cattle in Eastern Bhutan (Table 3) and was collected from all subdistricts. Our study findings agreed with a previous unpublished study, which found R. microplus to be the primary tick species in cattle in the surveyed areas of eastern Bhutan. Further, a tick survey conducted in the western region of Bhutan also found R. microplus to be the most prevalent tick species [9]. This tick species is also one of the most predominant tick species infesting livestock in India [47]. Moreover, this tick species is distributed throughout the world, especially in tropical and subtropical regions, and considered to be the most important tick species of cattle in the world [48,49]. Recently, it was recognized that R. microplus is a complex species that is comprised of at least five taxa: R. australis, R. annulatus, R. microplus clade A, R. microplus clade B, and R. microplus clade C [50,51]. However, the R. microplus we have collected was identified based on morphological keys [12,16], and further studies should be conducted to determine which clade of R. microplus it belongs to.
Rhipicephalus microplus was found more frequently in animals sampled in Pema Gatshel (98.2%) compared with those sampled in Trashigang (78.3%). Generally, this one-host tick species prefers warm and humid conditions [52], and temperature and rainfall are the most important climatic factors driving its geographic distribution [53,54]. Pema Gatshel, located at the warmer south, has the majority of its area dominated by subtropical climate with higher temperatures, abundant rainfall, and high relative humidity, while Trashigang is mainly characterized by a temperate climate with moderate temperature, rainfall, and relative humidity. This difference in the climatic conditions could be one factor that could have led to the greater prevalence of R. microplus in Pema Gatshel.
The infestation prevalence of R. microplus was found to decrease with increasing longitude (Table 6). Geographically, the longitudes of our study area ranged from 91° to 92° E. The areas between the longitudes 91° and 91.5° E are at an elevation of below 2000 m. Between the longitudes 91.5° and 92° E, the areas are at an elevation of more than 2000 m. This variation in elevation toward the east with increasing longitude could be the reason for the decrease in R. microplus infestation as R. microplus is known to avoid higher altitudes such as mountains and plateaus, where low temperatures are prevalent [55,56]. Elevation was also found as the most important environmental variable in the best habitat suitability model generated for R. microplus. The highest probability of its occurrence corresponded to the low elevation range (500–1000 masl) and decreased after that threshold. This elevation range corresponds to the humid subtropical zone (600–1200 masl) in Bhutan, characterized by an annual mean temperature of 19.5 °C and an annual rainfall of 1200–2500 mm [1]. This agrees with a study from Zimbabwe, where elevation was found as the most influential factor for the geographic distribution of R. microplus and R. decoloratus [54]. In Bhutan, areas with higher elevation are characterized by colder temperatures and low rainfalls, and that might result in such areas to be less suitable habitats for R. microplus.
A greater prevalence of R. microplus was observed in younger cattle compared with that of adult cattle (Table 6), which is consistent with the similar studies conducted in the neighboring Indian state of West Bengal [57], Bangladesh [58], and Nigeria [59]. However, studies conducted in Pakistan [60], Ethiopia [61], Egypt [43], and Nigeria [62] indicated lower prevalence in young cattle due to factors such as frequent grooming by dams and innate immunity. In our study area, those calves that are reared in the stall-fed systems are deprived of grooming as the farmers tether them at a safe distance from their dam to avoid milk suckling. Furthermore, young calves, especially the males, are the least attended by farmers, unlike the milking adult females. In most cattle owning households in the study area, calves and young heifers are often let out for grazing in nearby pastures and forests while milking adults are kept tethered around the homestead. Lack of maternal grooming and free access to pastures and forests likely increase the risk of tick infestation, and this could be the reason for the greater prevalence of R. microplus in younger animals in this study.
Rhipicephalus haemaphysaloides was the second most predominant tick species in this study (Table 3) and was collected from 20 subdistricts. This typical three-host tick has a limited geographical distribution in the world [63]; however, it is widely distributed in the neighboring Indian states bordering Bhutan [47] and China [64]. Unlike R. microplus, the prevalence of R. haemaphysaloides was greater in animals sampled from Trashigang (56.9%) compared with those sampled in Pema Gatshel (15.4%). It also showed a positive association with latitude indicating northerly distribution in the study area (Table 6). Trashigang, located toward the north in the study area, is at a higher elevation compared with Pema Gatshel. Temperatures and rainfalls are also comparatively low. Generally, temperature is considered to be the key climatic factor influencing the biological performance of this tick [63]. However, our findings suggest that R. haemaphysaloides might be tolerating cooler temperatures than other tick species at some point in their life cycle. In the habitat suitability modeling, R. haemaphysaloides presence was related to Bio 10 (temperature of the warmest quarter) and Bio 16 (precipitation of the wettest quarter) with the highest probability of occurrence at 25 °C and 1100 mm precipitation, respectively. Generally, the relative humidity plays a significant role in regulating the questing activity of ticks, and, further, the mortality of ticks also depends on the loss of water [65]. This tick species is an exophilic tick that may lose water while questing for a host, and thus its survival may be dependent on its ability to retain or gain water. Therefore, precipitation might play a critical role in the survival of this tick in the environment as the relative humidity is largely influenced by precipitation.
Haemaphysalis bispinosa was the third most predominant tick species identified in this study (Table 3) and was collected from 22 subdistricts. In South Asia, this tick is distributed in India, Nepal, Sri Lanka, Pakistan, Bangladesh, and Myanmar [66]. In India, it is widely distributed throughout the country and has been reported from the Indian states of Assam, Arunachal Pradesh, West Bengal, and Sikkim [17,47,67] that border Bhutan. In China, this tick has been considered to exist in the southern parts of the country, but most of the H. bispinosa reported in the Chinese literature are considered to be, in fact, H. longicornis [68]. Haemaphysalis spinigera was collected from 14 subdistricts in the study area. This tick is found in the foothills of the Central and Eastern Himalayan region, through Nepal to West Bengal in India [69]. It is widely distributed in India and has been reported from the Indian states of Assam and West Bengal [17,47], which border Bhutan.
The prevalence of both H. bispinosa and H. spinigera was greater in indigenous breeds compared to European breeds. This could be attributed to the difference in preferred management systems used for European and indigenous breeds. The European breeds are mostly reared in a stall-fed system with limited access to forests, whereas the indigenous breeds are mostly reared in a free-grazing system with easy access to forests. The risk of exposure to ticks becomes greater when cattle graze in the forests. Studies have observed a higher prevalence of tick infestation in grazing cattle [58,70,71]. Furthermore, indigenous breeds of cattle in Bhutan are generally considered resistant to many pests and diseases, and the management practices like grooming, brushing, and acaricide application, are rarely practiced by the farmers [7].
Haemaphysalis bispinosa infestation increased with increasing altitude while H. spinigera infestation decreased with increasing altitude (Table 6). This might indicate that H. bispinosa can tolerate a cooler temperature at some point in its life cycle while H. spinigera prefers low altitude subtropical areas. However, altitude is not the only driving factor in the distribution of this tick species. Climatic factors such as temperature, rainfall, relative humidity, vegetation, and host availability also influence its distribution [72]. In Bhutan, the high mountains and narrow valleys cause an extreme variation in climate, and temperature is mainly affected by altitude and precipitation by latitude [45]. Further, the types of vegetation in the country are thought to be determined by the terrain and local climatic conditions; therefore, in the absence of detailed data on these various climatic and environmental factors, it is difficult to interpret what factors largely drive the distribution of ticks in Bhutan. Nevertheless, altitude could be one of the main driving factors determining the habitat suitability of tick species as it is believed to change the temperature by 0.5 °C with every 100 m change in altitude [45].
The best habitat suitability model generated for H. bispinosa showed that Bio 18 (precipitation of the warmest quarter) and Bio 16 (precipitation of the wettest quarter) were the two variables, besides land cover, predicting its distribution. For H. spinigera, Bio 16 (precipitation of the wettest quarter) and Bio 19 (precipitation of the coldest quarter) were the two important variables, besides land cover. This indicates that the distribution of both Haemaphysalis species was related to precipitation variables and land cover. The pattern of influence by Bio 16 (precipitation of the wettest quarter) was similar for both species—the probability of occurrence increased at 400 mm and then decreased when Bio 16 exceeded 1600 mm. Other studies from the neighboring countries like Bangladesh and India had similar findings for these two Haemaphysalis species. Haemaphysalis bispinosa seems to thrive well in areas with high summer rainfall in Bangladesh [73]. In India, the population of the engorged adults increased during the monsoon seasons when there is a lot of rain [17]. Similarly, H. spinigera is also known to prefer wet and moist habitats with heavy to moderate rainfall [17]. In South India, where Kysanuar forest disease (KFD) is endemic, the nymphal activity of H. spinigera (the principal vector of KFD) peaks in spring and summer seasons coinciding with the onset of monsoons [74]. Further, a habitat suitability modeling study (preprint) [75] from South India for H. spinigera has identified precipitation variables such as annual precipitation and wettest month precipitation as the important climatic variables influencing its geographical distribution. These findings indicate that precipitation is one of the important climatic factors that determine the distribution of H. bispinosa and H. spinigera.
The habitat suitability models developed indicated that the land cover classes such as kamzhing, meadows, chuzhing, and shrublands were suitable habitats for the tick species modeled. This might be associated with the presence of a host in these types of land. Kamzhing is the most important land type for dry land agricultural farming in Bhutan. Besides cultivating cereal crops, kamzhing is also used for growing fodder crops like oats and clover, especially in winters. Cattle are usually tethered and allowed to graze on the small patches of grasses and weeds in the kamzhing areas, and this might explain why kamzhing is a suitable land type for ticks, especially for a one-host tick like R. microplus. The high probability of tick(s) occurrence in the meadows could also be due to host availability. Most of the indigenous cattle breeds reared in a free grazing system are grazed in the meadows, and frequently, wild cervids like deer also graze in the meadows. Chuzhing, at least theoretically, should not be a suitable habitat as ticks would not be able to withstand flooding. However, the entire area categorized as chuzhing is not under rice cultivation. The insufficient irrigation supply, farm labor shortage, erratic rainfall, and human–wildlife conflicts, among many other factors, have pushed many farmers to leave their chuzhing fallow [76]. Consequently, these fallow lands have become a grazing area for domestic animals like cattle and wild cervids like deer. Occasionally, some farmers also grow fodder grasses on chuzhings that could not be utilized for rice cultivation. This changing pattern of land use might be the reason for chuzhing featuring as one of the suitable habitats for ticks. Shrublands were also identified as a suitable land type for some of the tick species modeled. Generally, shrublands in Bhutan are suitable for tick–host interactions as they are normal grazing areas for cattle and wild herbivores such as deer. Further, shrublands provide a conducive environment for small mammals, especially rodents, which can serve as reservoir hosts for the larvae and nymphs of multihost tick species.
For all tick species modeled, the northeastern and the southernmost parts of the study area were predicted to be unsuitable areas. The northeastern region is the area with higher elevations where temperatures and rainfall are low. In this region, the vegetation is predominantly mixed conifer forests. The low temperatures, scanty rainfall, and the occasional snow and frosts might be the limiting factors for tick survival. However, there might be some other species of ticks such as Hyalomma spp., which can tolerate high elevations. Hyalomma spp. were collected in Eastern Bhutan in areas with high elevation (above 2700 masl) (Dr. Susan C. Cork, personal communication, 2018). Therefore, additional studies using other sampling methods such as flagging and dragging may have to be conducted for at least a year to ascertain whether these areas are unsuitable for ticks. The southernmost part predicted as unsuitable areas are the subtropical floodplains at a very low elevation (less than 500 masl). The vegetation in these areas is predominantly broadleaf forests. Climatically, these areas will be suitable for ticks, but the host unavailability (especially the bovines) might be the limiting factor. There are no human activities such as settlement and farming in these areas. The frequent flooding during the monsoon will also be one of the limiting factors. However, there are some wild animals inhabiting these areas, which might be potential hosts for some species of ticks. Therefore, future efforts must be targeted toward tick collection from wild animals and the environment in these areas to understand the habitat suitability better.
Seven specimens of A. testudinarium were collected from seven subdistricts in Pema Gatshel, and one unidentified specimen of Ixodes was collected from a subdistrict in Trashigang. Amblyomma testudinarium is considered to be a rare tick species reported only from Asian countries like Malaysia, India, Japan, Korea, and China [47,77,78]. In India, it has been collected from cattle, mithuns, yaks, and wild animals such as the tiger, wild boar, barking deer, and elephant in Assam, Arunachal Pradesh, Mizoram, and West Bengal [47,79]. These reports suggest that this tick species is predominantly found on the animals that live in and around the forests of the Himalayan foothills. In this study, A. testudinarium was found only in the southeastern district of Pema Gatshel. Since this study lacked year-round sampling and was conducted only in cattle, it is difficult to interpret why A. testudinarium was not collected from Trashigang. As with Ixodes, there are 11 valid species reported from India [47], 24 valid species from China [19], and 15 valid species from Nepal [80]. This tick can infest a wide range of hosts, including domestic animals and humans, small mammals, wildlife such as deer, birds, and reptiles [80]. However, in this study, only one specimen of Ixodes was found on cattle. This could have been biased due to our sampling procedure as tick collection was done only from the cattle in May–June 2019.
This study has some limitations. Since the primary aim of our study was to understand the presence and diversity of tick species in cattle, our sampling was biased toward areas dominated by cattle rearing. The presence data used in this modeling study were based on tick collection from cattle, and it may not represent the exact location from where tick(s) originated. Therefore, future studies should include tick sampling from the environment and other hosts. The sampling for this study was conducted only once in May and June 2019 as this time period is thought to be the season for peak tick infestation in cattle in Bhutan. More species might have been found if the sampling was conducted throughout the year. To reduce these knowledge gaps, future studies should focus on active surveillance for some years targeting a wider range of hosts. Other sampling methods like flagging and dragging on vegetation should also be conducted to improve the knowledge of tick species diversity. The three-host ticks like R. haemaphysaloides, H. bispinosa, and H. spinigera depend on multiple hosts to complete their life-cycle; however, we collected ticks only from cattle, and the status of other potential hosts that these ticks might be feeding on is not known. In Bhutan, tick collection from a wider range of hosts could be done in collaboration with wildlife officials who often encounter wild animals and birds during rescue and relocation operations. The resampled bioclimatic variables may not represent the extreme climatic variations in the mountainous terrain of Bhutan, thereby increasing the uncertainty of the models. Therefore, more local climate or weather data should be used for future modeling studies. Further, the MaxEnt approach, as a correlative species distribution model, estimates the realized niche but not the fundamental niche of the species [81]; therefore, potential habitat is likely to be more than what is predicted by the habitat suitability maps. Nevertheless, MaxEnt is a powerful software tool for habitat suitability modeling with several advantages that include the need for few presence points and the ability to use both continuous and categorical variables [10].

5. Conclusions

The habitat suitability models developed in this study identified some potential climatic and environmental factors that can be used to predict the current distribution of common tick species in Eastern Bhutan. Land cover types such as kamzhing, meadows, chuzhing, and shrubland were identified as the preferred habitat of these tick species. Elevation was found to be the most important environmental variable influencing the distribution of R. microplus. Bio 10 (temperature of the warmest quarter) and Bio 16 (precipitation of the wettest quarter) were identified as the most important climatic variables driving the distribution of R. haemaphysaloides. Precipitation variables (i.e., precipitation of the warmest, the coldest, and the wettest quarters) were identified as the important climatic variables for two Haemaphysalis species (i.e., H. bispinosa and H. spinigera). Further, the habitat suitability maps for the tick species studied can be used for planning a targeted tick surveillance in Eastern Bhutan. Overall, the findings are expected to inform surveillance plans and help in the development and implementation of a targeted tick control program and subsequent prevention and control programs for ticks and tick-borne diseases in cattle in Bhutan.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2414-6366/6/1/27/s1: Figure S1. Response curve plotting the probability of R. microplus occurrence in Eastern Bhutan against the values of the top environmental variables: (A) elevation (DEM_SRTM) and (B) land use and land cover 2016 (LULC). The X-axis represents the variable value, and the Y-axis represents the probability of presence as predicted by the best MaxEnt model RM_B1. Figure S2. Response curves plotting the probability of R. haemaphysaloides occurrence in Eastern Bhutan against the values of the top environmental variables: (A) temperature of the warmest quarter (Bio 10), (B) land use and land cover 2016 (LULC), and (C) precipitation of the wettest quarter (Bio 16). The X-axis represents the variable value, and the Y-axis represents the probability of presence as predicted by the best MaxEnt model RH_D1. Figure S3. Response curves plotting the probability of H. bispinosa occurrence in Eastern Bhutan against the values of the top environmental variables: (A) precipitation of the warmest quarter (Bio 18), (B) land use and land cover 2016 (LULC), and (C) precipitation of the wettest quarter (Bio 16). The X-axis represents the variable value, and the Y-axis represents the probability of presence as predicted by the best MaxEnt model HB_C2. Figure S4. Response curves plotting the probability of H. spinigera occurrence in Eastern Bhutan against the values of the top environmental variables: (A) precipitation of the coldest quarter (Bio 19), (B) land use and land cover 2016 (LULC), and (C) precipitation of the wettest quarter (Bio 16). The X-axis represents the variable value, and the Y-axis represents the probability of presence as predicted by the best MaxEnt model HS_E1. Table S1. Results of simple logistic regression describing the relationships between infestation (infested or not) and cattle and geographic parameters. Table S2. Variable contribution and permutation importance values obtained during the first MaxEnt run performed with all the variables for predicting tick species presence. Only variables that achieved values of more than 1% in both metrics are shown. Table S3. Correlation matrix showing Spearman correlation coefficient (rs) for variables that achieved more than 1% contribution and permutation importance in the first MaxEnt run for all tick species. Table S4. MaxEnt models developed for predicting tick species occurrence in Eastern Bhutan. Table shows variables included into each different model. Variables, percent contribution (%), permutation importance (PI), AUC (training and testing), correlation (training and testing), corrected AIC, delta (Δ), Akaike weights (ω), and number of parameters (P) for each different model are given. The best model for each tick species is highlighted in bold. Table S5. Study area (in km2 and percentage) classified by the probability of tick species occurrence in Eastern Bhutan: RM_B1 for R. microplus; RH_D1 for R. haemaphysaloides; HB_C2 for H. bispinosa; and HS_E1 for H. spinigera. Dataset S1. Data (excel file) collected during the survey and used for analysis.

Author Contributions

Conceptualization, J.N., T.J.L., T.T., R.B.G., S.D., and S.C.C.; methodology, J.N., T.J.L., I.C., S.C., R.B.G., T.T., S.D., and S.C.C.; software, J.N., T.J.L., and I.C.; validation, T.J.L., T.T., R.B.G., S.D., S.C., and S.C.C.; formal analysis, J.N., T.J.L., and I.C.; investigation, J.N.; resources, S.C.C.; data curation, J.N., T.J.L., S.C., and I.C.; writing—original draft preparation, J.N.; writing—review and editing, J.N., T.T., T.J.L., I.C., S.C., and S.C.C.; visualization, I.C., T.T., T.J.L., S.C., S.D., and S.C.C.; supervision, S.C.C.; project administration, S.C.C.; funding acquisition, S.C.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by One Health Graduate Training Award, Faculty of Veterinary Medicine, University of Calgary, Canada and International Research Grant, University of Calgary, Canada. The Department of Livestock, Ministry of Agriculture and Forests, Royal Government of Bhutan provided in-kind support during the fieldwork.

Institutional Review Board Statement

The study protocol was approved by both the Animal Care Committee (ACC), University of Calgary, Canada (AC 19-0035) and the Research and Extension Division, Department of Livestock, Ministry of Agriculture and Forests, Royal Government of Bhutan (Animal Research Application Form-15/05/2019).

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in supplementary material here.

Acknowledgments

This work is an excellent example of the power of collaboration, and we are grateful to many institutions and individuals that helped with the field survey and laboratory work. The institutions include the National Centre for Animal Health, Thimphu, Bhutan; Regional Livestock Development Centre, Kanglung, Trashigang, Bhutan; Satellite Veterinary Laboratory, Nganglam, Pema Gatshel, Bhutan; and the Dzongkhag Livestock Sector, Trashigang and Pema Gatshel, Bhutan. We thank Carl S. Ribble and Shaun Dergousoff for their insights and guidance to this work. We thank Karma Phuntsho, Lungten, Tashi Norbu, Lhatru, Phurpa Wangdi, Kinley Tenzin and Rinzin Namgay for assisting us with the tick collection works, and Tenzinla and Tshewang Dema for assisting us with the countless hours of laboratory work. We also thank Karma Jigme for helping us with photography works. Finally, we thank our farmers who at every selected household helped us with tick collection work.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Ministry of Agriculture and Forests. MoAF BHUTAN RNR STATISTICS–2018. Available online: http://www.moaf.gov.bt/bhutan-rnr-statistics-2018/ (accessed on 9 March 2020).
  2. National Statistics Bureau. 2017 Population & Housing Census of Bhutan: National Report = ʼBrug Gi Mi Rlobs Dang Khyim Gyi Grangs Rtsis. 2017. Available online: http://www.nsb.gov.bt/publication/files/PHCB2017_national.pdf (accessed on 9 March 2020).
  3. NCHM Bhutan State of the Climate 2017. Available online: http://www.nchm.gov.bt/attachment/ckfinder/userfiles/files/Bhutan%20State%20of%20the%20Climate%202017.pdf (accessed on 9 March 2020).
  4. Banerjee, A.; Bandopadhyay, R. Biodiversity Hotspot of Bhutan and Its Sustainability. Curr. Sci. 2016, 110, 521–527. [Google Scholar] [CrossRef]
  5. DoL Livestock Statistics. Available online: http://www.dol.gov.bt/?page_id=394 (accessed on 9 March 2020).
  6. Hidano, A.; Dukpa, K.; Rinzin, K.; Sharma, B.; Dahal, N.; Stevenson, M.A. A Cross-Sectional Survey of Population Demographics, the Prevalence of Major Disease Conditions and Reason-Specific Proportional Mortality of Domestic Cattle in the Kingdom of Bhutan. Prev. Vet. Med. 2016, 130, 1–9. [Google Scholar] [CrossRef] [PubMed]
  7. Phanchung; Dorji, P.; Sonam, T.; Pelden, K. Smallholder Dairy Farming in Bhutan: Characteristics, Constraints, and Development Opportunities, Smallholder Dairy in Mixed Farming Systems of the HKH; International Centre for Integrated Mountain Development: Kathmandu, Nepal, 2002; pp. 19–34. [Google Scholar]
  8. NCAH Veterinary Information System. Available online: http://vis.ncah.gov.bt/php/index.php (accessed on 15 February 2020).
  9. RLDC Wangdue 9th Annual Progress Report of RLDC Wangdue (2018–2019). Available online: http://www.dol.gov.bt/?wpdmpro=rldc-wangdue-annual-progress-report-fy-2018-19. (accessed on 9 March 2020).
  10. Phillips, S.J.; Anderson, R.P.; Schapire, R.E. Maximum Entropy Modeling of Species Geographic Distributions. Ecol. Model. 2006, 190, 231–259. [Google Scholar] [CrossRef] [Green Version]
  11. Thrusfield, M. Veterinary Epidemiology, 3rd ed.; John Wiley & Sons: Hoboken, NJ, USA, 2018; ISBN 978-1-118-28028-7. [Google Scholar]
  12. Matthysee, J.G.; Colbo, M.H. The Ixodid Ticks of Uganda: Together With Species Pertinent to Unganda Because of Their Present Known Distribution; Entomological Society of America: College Park, MD, USA, 1987; ISBN 978-0-938522-31-7. [Google Scholar]
  13. Walker, A.R.; Bouattour, A.; Camicas, J.-L.; Estrada-Pena, A.; Horak, I.G.; Latif, A.A.; Pegram, R.G.; Preston, P.M. Ticks of Domestic Animals in Africa: A Guide to Identification of Species; Bioscience Reports; The University of Edinburgh: Edinburgh, Scotland, 2003; ISBN 0-9545173-0-X. [Google Scholar]
  14. Anastos, G. Scutate Ticks, or Ixodidae, of Indonesia; Brooklyn Entomological Society: Brooklyn, NY, USA, 1950. [Google Scholar]
  15. Robinson, L.E.; Nuttall, G.H.F.; Warburton, C. The Genus Amblyomma; Cambridge University Press: Cambridge, UK, 1926. [Google Scholar]
  16. Estrada-Pena, A.; Venzal, J.M.; Nava, S.; Mangold, A.; Guglielmone, A.A.; Labruna, M.B.; de la Fuente, J. Reinstatement of Rhipicephalus (Boophilus) australis (Acari: Ixodidae) With Redescription of the Adult and Larval Stages. J. Med. Entomol. 2012, 49, 794–802. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Geevarghese, G.; Mishra, A. Haemaphysalis Ticks of India, 1st ed.; Elsevier: Amsterdam, The Netherlands, 2011; ISBN 978-0-12-387812-0. [Google Scholar]
  18. Nutall, G.H.; Warburton, C.; Cooper, W.; Robinson, L.T. A Monograph of the Ixodoidea; Part III; Cambridge University Press: Cambridge, UK, 1915; pp. 349–550. [Google Scholar]
  19. Guo, Y.; Sun, Y.; Xu, R. The Genus Ixodes (Acari: Ixodidae) in China with Three New Record Species. Acta Parasitol. 2016, 61, 729–742. [Google Scholar] [CrossRef] [PubMed]
  20. Hijmans, R.J.; Cameron, S.E.; Parra, J.L.; Jones, P.G.; Jarvis, A. Very High Resolution Interpolated Climate Surfaces for Global Land Areas. Int. J. Climatol. 2005, 25, 1965–1978. [Google Scholar] [CrossRef]
  21. FRMD Land Use and Land Cover of Bhutan 2016, Maps and Statistics; Forest Resources Management Division Department of Forests & Park Services Ministry of Agriculture and Forests Thimphu: Thimphu, Bhutan, 2017; ISBN 978-99936-743-2-0.
  22. Hijmans, R.J. raster: Geographic Data Analysis and Modeling, R package version 3.0-12. 2020. Available online: https://CRAN.R-project.org/package=raster (accessed on 14 May 2020).
  23. Bivand, R.; Rundel, C. rgeos: Interface to Geometry Engine—Open Source (’GEOS’); R package version 0.5-2. 2019. Available online: https://CRAN.R-project.org/package=rgeos (accessed on 14 May 2020).
  24. Bivand, R.; Keitt, T.; Rowlingson, B. rgdal: Bindings for the “Geospatial” Data Abstraction Library; R package version 1.4-8. 2019. Available online: https://CRAN.R-project.org/package=rgdal (accessed on 14 May 2020).
  25. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018. [Google Scholar]
  26. Bush, A.O.; Lafferty, K.D.; Lotz, J.M.; Shostak, A.W. Parasitology Meets Ecology on Its Own Terms: Margolis. Revisited on JSTOR. Available online: https://www-jstor-org.ezproxy.lib.ucalgary.ca/stable/3284227?origin=crossref&seq=1#metadata_info_tab_contents (accessed on 18 March 2020).
  27. Scherer, R. PropCIs: Various Confidence Interval Methods for Proportions; R package version 0.3-0. 2018. Available online: https://CRAN.R-project.org/package=PropCIs (accessed on 13 March 2020).
  28. Harrell, F.E., Jr. with contributions from Charles Dupont and many others. Hmisc: Harrell Miscellaneous. R package version 4.4-0. 2020. Available online: https://CRAN.R-project.org/package=Hmisc (accessed on 19 April 2020).
  29. Allison, P. When Can You Safely Ignore Multicollinearity? Available online: https://statisticalhorizons.com/multicollinearity (accessed on 15 April 2020).
  30. Dardis, C. LogisticDx: Diagnostic Tests for Models with a Binomial Response; R package version 0.2. 2015. Available online: https://CRAN.R-project.org/package=LogisticDx (accessed on 25 January 2020).
  31. Fox, J.; Weisberg, S. An R Companion to Applied Regression, 3rd ed.; Sage publications: New York, NY, USA, 2019. [Google Scholar]
  32. Hijmans, R.J.; Phillips, S.; Leathwick, J.; Elith, J. dismo: Species Distribution Modeling; R package version 1.1-4. 2017. Available online: https://CRAN.R-project.org/package=dismo (accessed on 14 May 2020).
  33. Feng, X.; Gebresenbet, F.; Walker, C. Shifting from Closed-Source Graphical-Interface to Open-Source Programming Environment: A Brief Tutorial on Running Maxent in R.; No. e3346v1; PeerJ Inc.: Corte Madera, CA, USA, 2017; Available online: https://peerj.com/preprints/3346/ (accessed on 14 May 2020).
  34. Pebesma, E.J.; Bivand, R. S Classes and Methods for Spatial Data in R. R News 2005, 5, 9–13. [Google Scholar]
  35. Zuliani, A.; Massolo, A.; Lysyk, T.; Johnson, G.; Marshall, S.; Berger, K.; Cork, S.C. Modelling the Northward Expansion of Culicoides sonorensis (Diptera: Ceratopogonidae) under Future Climate Scenarios. PLoS ONE 2015, 10. [Google Scholar] [CrossRef] [PubMed]
  36. Revelle, W. psych: Procedures for Psychological, Psychometric, and Personality Research, R package version 1.9.12; Northwestern University: Evanston, IL, USA, 2019; Available online: https://CRAN.R-project.org/package=psych (accessed on 14 May 2020).
  37. Burnham, K.P.; Anderson, D.R. Multimodel Inference: Understanding AIC and BIC in Model Selection. Sociol. Methods Res. 2004, 33, 261–304. [Google Scholar] [CrossRef]
  38. Muscarella, R.; Galante, P.J.; Soley-Guardia, M.; Boria, R.A.; Kass, J.M.; Uriarte, M.; Anderson, R.P. ENMeval: An R Package for Conducting Spatially Independent Evaluations and Estimating Optimal Model Complexity for Maxent Ecological Niche Models. Methods Ecol. Evol. 2014, 5, 1198–1205. [Google Scholar] [CrossRef]
  39. Fielding, A.H.; Bell, J.F. A Review of Methods for the Assessment of Prediction Errors in Conservation Presence/Absence Models. Environ. Conserv. 1997, 24, 38–49. [Google Scholar] [CrossRef]
  40. Warren, D.L.; Seifert, S.N. Ecological Niche Modeling in Maxent: The Importance of Model Complexity and the Performance of Model Selection Criteria. Ecol. Appl. 2011, 21, 335–342. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Perpi, O.; Hijmans, R. rasterVis; R package version 0.47. 2019. Available online: http://oscarperpinan.github.io/rastervis/ (accessed on 14 May 2020).
  42. Neuwirth, E. RColorBrewer: ColorBrewer Palettes; R Package Version 1.1-2. 2014. Available online: {https://CRAN.R-project.org/package=RColorBrewer} (accessed on 22 May 2020).
  43. Asmaa, N.M.; ElBably, M.A.; Shokier, K.A. Studies on Prevalence, Risk Indicators and Control Options for Tick Infestation in Ruminants. Beni-Suef Univ. J. Basic Appl. Sci. 2014, 3, 68–73. [Google Scholar] [CrossRef] [Green Version]
  44. Bianchi, M.W.; Barré, N.; Messad, S. Factors Related to Cattle Infestation Level and Resistance to Acaricides in Boophilus microplus Tick Populations in New Caledonia. Vet. Parasitol. 2003, 112, 75–89. [Google Scholar] [CrossRef]
  45. Dorji, U.; Olesen, J.E.; Bocher, P.K.; Seidenkrantz, M.S. Spatial Variation of Temperature and Precipitation in Bhutan and Links to Vegetation and Land Cover. Mt. Res. Dev. 2016, 36, 66–79. [Google Scholar] [CrossRef] [Green Version]
  46. Elith, J.; Phillips, S.J.; Hastie, T.; Dudík, M.; Chee, Y.E.; Yates, C.J. A Statistical Explanation of MaxEnt for Ecologists. Divers. Distrib. 2019, 43–57. [Google Scholar] [CrossRef]
  47. Ghosh, S.; Bansal, G.C.; Gupta, S.C.; Ray, D.; Khan, M.Q.; Irshad, H.; Shahiduzzaman, M.; Seitzer, U.; Ahmed, J.S. Status of Tick Distribution in Bangladesh, India and Pakistan. Parasitol. Res. 2007, 101, 207–216. [Google Scholar] [CrossRef] [PubMed]
  48. Jongejan, F.; Uilenberg, G. The Global Importance of Ticks. Parasitology 2004, 129, S3–S14. [Google Scholar] [CrossRef] [PubMed]
  49. Spickler, A.R. Rhipicephalus (Boophilus) Microplus. 2007. Available online: http://www.cfsph.iastate.edu/DiseaseInfo/factsheets.php (accessed on 22 May 2020).
  50. Burger, T.D.; Shao, R.; Barker, S.C. Phylogenetic Analysis of Mitochondrial Genome Sequences Indicates That the Cattle Tick, Rhipicephalus (Boophilus) microplus, Contains a Cryptic Species. Mol. Phylogenet. Evol. 2014, 76, 241–253. [Google Scholar] [CrossRef]
  51. Low, V.L.; Tay, S.T.; Kho, K.L.; Koh, F.X.; Tan, T.K.; Lim, Y.A.L.; Ong, B.L.; Panchadcharam, C.; Norma-Rashid, Y.; Sofian-Azirun, M. Molecular Characterisation of the Tick Rhipicephalus microplus in Malaysia: New Insights into the Cryptic Diversity and Distinct Genetic Assemblages throughout the World. Parasit. Vectors 2015, 8, 341. [Google Scholar] [CrossRef] [Green Version]
  52. De Clercq, E.M.; Vanwambeke, S.O.; Sungirai, M.; Adehan, S.; Lokossou, R.; Madder, M. Geographic Distribution of the Invasive Cattle Tick Rhipicephalus microplus, a Country-Wide Survey in Benin. Exp. Appl. Acarol. 2012, 58, 441–452. [Google Scholar] [CrossRef] [Green Version]
  53. Estrada-Peña, A.; Acedo, C.S.; Quílez, J.; Del Cacho, E. A Retrospective Study of Climatic Suitability for the Tick Rhipicephalus (Boophilus) microplus in the Americas. Glob. Ecol. Biogeogr. 2005, 14, 565–573. [Google Scholar] [CrossRef]
  54. Sungirai, M.; Moyo, D.Z.; De Clercq, P.; Madder, M.; Vanwambeke, S.O.; De Clercq, E.M. Modelling the Distribution of Rhipicephalus microplus and R. decoloratus in Zimbabwe. Vet. Parasitol. Reg. Stud. Rep. 2018, 14, 41–49. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Estrada-Peña, A.; García, Z.; Sánchez, H.F. The Distribution and Ecological Preferences of Boophilus microplus (Acari: Ixodidae) in Mexico. Exp. Appl. Acarol. 2006, 38, 307–316. [Google Scholar] [CrossRef] [PubMed]
  56. Lynen, G.; Zeman, P.; Bakuname, C.; Di Giulio, G.; Mtui, P.; Sanka, P.; Jongejan, F. Shifts in the Distributional Ranges of Boophilus Ticks in Tanzania: Evidence That a Parapatric Boundary between Boophilus microplus and B. decoloratus Follows Climate Gradients. Exp. Appl. Acarol. 2008, 44, 147–164. [Google Scholar] [CrossRef] [PubMed]
  57. Debbarma, A.; Pandit, S.; Jas, R.; Baidya, S.; Mandal, S.C.; Jana, P.S. Prevalence of Hard Tick Infestations in Cattle of West Bengal, India. Biol. Rhythm Res. 2018, 49, 655–662. [Google Scholar] [CrossRef]
  58. Kabir, M.H.B.; Mondal, M.M.H.; Eliyas, M.; Mannan, M.A.; Hashem, M.A.; Debnath, N.C.; Miazi, O.F.; Mohiuddin, C.; Kashem, M.A.; Islam, M.R.; et al. An Epidemiological Survey on Investigation of Tick Infestation in Cattle at Chittagong District, Bangladesh. Afr. J. Microbiol. Res. 2011, 5, 346–352. [Google Scholar] [CrossRef]
  59. Eyo, J.E.; Ekeh, F.N.; Ivoke, N.; Atama, C.I.; Onah, I.E.; Ezenwaji, N.E.; Ikele, C.B. Survey of Tick Infestation of Cattle at Four Selected Grazing Sites in the Tropics. Glob. Vet. 2014, 12, 479–486. [Google Scholar] [CrossRef]
  60. Rehman, A.; Nijhof, A.M.; Sauter-Louis, C.; Schauer, B.; Staubach, C.; Conraths, F.J. Distribution of Ticks Infesting Ruminants and Risk Factors Associated with High Tick Prevalence in Livestock Farms in the Semi-Arid and Arid Agro-Ecological Zones of Pakistan. Parasit. Vectors 2017, 10, 190. [Google Scholar] [CrossRef]
  61. Kemal, J.; Tamerat, N.; Tuluka, T. Infestation and Identification of Ixodid Tick in Cattle: The Case of Arbegona District, Southern Ethiopia. J. Vet. Med. 2016, 2016. [Google Scholar] [CrossRef] [PubMed]
  62. Lorusso, V.; Picozzi, K.; de Bronsvoort, B.M.; Majekodunmi, A.; Dongkum, C.; Balak, G.; Igweh, A.; Welburn, S.C. Ixodid Ticks of Traditionally Managed Cattle in Central Nigeria: Where Rhipicephalus (Boophilus) microplus Does Not Dare (Yet?). Parasit. Vectors 2013, 6, 171. [Google Scholar] [CrossRef] [Green Version]
  63. Diyes, G.C.P.; Apanaskevich, D.A.; Rajakaruna, R.S. Lifecycle of Rhipicephalus haemaphysaloides under Laboratory Conditions. Med. Vet. Entomol. 2017, 31, 327–332. [Google Scholar] [CrossRef] [PubMed]
  64. Chen, Z.; Yang, X.; Bu, F.; Yang, X.; Yang, X.; Liu, J. Ticks (Acari: Ixodoidea: Argasidae, Ixodidae) of China. Exp. Appl. Acarol. 2010, 51, 393–404. [Google Scholar] [CrossRef] [PubMed]
  65. Estrada-Peña, A. Ticks as Vectors: Taxonomy, Biology and Ecology: -EN- -FR- Les Tiques En Tant Que Vecteurs: Taxonomie, Biologie et Écologie -ES- Las Garrapatas Como Vectores: Taxonomía, Biología y Ecología. Rev. Sci. Tech. OIE 2015, 34, 53–65. [Google Scholar] [CrossRef] [Green Version]
  66. Hoogstraal, H.; Lim, B.-L.; Anastos, G. Haemaphysalis (Kaiseriana) bispinosa Neumann (Ixodoidea: Ixodidae): Evidence for Consideration as an Introduced Species in the Malay Peninsula and Borneo. J. Parasitol. 1969, 55, 1075–1077. [Google Scholar] [CrossRef] [PubMed]
  67. Brahma, R.K.; Dixit, V.; Sangwan, A.K.; Doley, R. Identification and Characterization of Rhipicephalus (Boophilus) microplus and Haemaphysalis bispinosa Ticks (Acari: Ixodidae) of North East India by ITS2 and 16S RDNA Sequences and Morphological Analysis. Exp. Appl. Acarol. 2014, 62, 253–265. [Google Scholar] [CrossRef]
  68. Chen, Z.; Li, Y.; Ren, Q.; Liu, Z.; Luo, J.; Li, K.; Guan, G.; Yang, J.; Han, X.; Liu, G.; et al. Does Haemaphysalis bispinosa (Acari: Ixodidae) Really Occur in China? Exp. Appl. Acarol. 2015, 65, 249–257. [Google Scholar] [CrossRef]
  69. Ghalsasi, G.R.; Dhanda, V. Taxonomy and Biology of Haemaphysalis (Kaiseriana) spinigera (Acarina: Ixodidae). Orient. Insects 1974, 8, 505–520. [Google Scholar] [CrossRef]
  70. Ogden, N.H.; Swai, E.; Beauchamp, G.; Karimuribo, E.; Fitzpatrick, J.L.; Bryant, M.J.; Kambarage, D.; French, N.P. Risk Factors for Tick Attachment to Smallholder Dairy Cattle in Tanzania. Prev. Vet. Med. 2005, 67, 157–170. [Google Scholar] [CrossRef] [PubMed]
  71. Tiki, B.; Addis, M. Distribution of Ixodid Ticks on Cattle in and Around Holeta Town, Ethiopia. Glob. Vet. 2011, 7, 527–531. [Google Scholar]
  72. Estrada-Peña, A.; de la Fuente, J. The Ecology of Ticks and Epidemiology of Tick-Borne Viral Diseases. Antivir. Res. 2014, 108, 104–128. [Google Scholar] [CrossRef]
  73. Islam, M.K.; Alim, M.A.; Tsuji, N.; Mondal, M.M.H. An Investigation into the Distribution, Host-Preference and Populationdensity of Ixodid Ticks Affecting Domestic Animals in Bangladesh. Trop. Anim. Health Prod. 2006, 38, 485–490. [Google Scholar] [CrossRef] [PubMed]
  74. Ajesh, K.; Nagaraja, B.K.; Sreejith, K. Kyasanur Forest Disease Virus Breaking the Endemic Barrier: An Investigation into Ecological Effects on Disease Emergence and Future Outlook. Zoonoses Public Health 2017, 64, e73–e80. [Google Scholar] [CrossRef]
  75. Pramanik, M.; Singh, P.; Dhiman, R. Identification of Bio-Climatic Determinants and Potential Risk Areas for Kyasanur Forest Disease in Southern India Using MaxEnt Modelling Approach. BMC Infect. Dis. 2020. [Google Scholar] [CrossRef] [Green Version]
  76. Kuensel Rice Cultivation Area Shrinks–Kuensel Online. Available online: https://kuenselonline.com/rice-cultivation-area-shrinks/ (accessed on 24 May 2020).
  77. Chao, L.-L.; Lu, C.-W.; Lin, Y.-F.; Shih, C.-M. Molecular and Morphological Identification of a Human Biting Tick, Amblyomma testudinarium (Acari: Ixodidae), in Taiwan. Exp. Appl. Acarol. 2017, 71, 401–414. [Google Scholar] [CrossRef]
  78. Zheng, W.; Xuan, X.; Fu, R.; Tao, H.; Xu, R.; Liu, Y.; Liu, X.; Jiang, J.; Wu, H.; Ma, H.; et al. Preliminary Investigation of Ixodid Ticks in Jiangxi Province of Eastern China. Exp. Appl. Acarol. 2019, 77, 93–104. [Google Scholar] [CrossRef]
  79. Chamuah, J.K.; Bhattacharjee, K.; Sarmah, P.C.; Raina, O.K.; Mukherjee, S.; Rajkhowa, C. Report of Amblyomma testudinarium in Mithuns (Bos frontalis) from Eastern Mizoram (India). J. Parasit. Dis. Off. Organ Indian Soc. Parasitol. 2016, 40, 1217–1220. [Google Scholar] [CrossRef] [Green Version]
  80. Clifford, C.M.; Hoogstraal, H.; Keirans, J.E. The Ixodes Ticks (Agarina: Ixodidae) of Nepal. J. Med. Entomol. 1975, 12, 115–137. [Google Scholar] [CrossRef] [PubMed]
  81. Kearney, M.R.; Wintle, B.A.; Porter, W.P. Correlative and Mechanistic Models of Species Distribution Provide Congruent Forecasts under Climate Change. Conserv. Lett. 2010, 3, 203–213. [Google Scholar] [CrossRef]
Tropicalmed 06 00027 g001 550
Figure 1. (a) Map of Bhutan showing the study area (pink shade) and; (b) Elevation map of the study area showing tick sampling sites (black dots). The map was prepared using Quantum GIS, QGIS Development Team (2019), QGIS Geographic Information System, Open-Source Geospatial Foundation Project (http://qgis.osgeo.org (accessed on 28 January 2021)) and was not taken from another source.
Figure 1. (a) Map of Bhutan showing the study area (pink shade) and; (b) Elevation map of the study area showing tick sampling sites (black dots). The map was prepared using Quantum GIS, QGIS Development Team (2019), QGIS Geographic Information System, Open-Source Geospatial Foundation Project (http://qgis.osgeo.org (accessed on 28 January 2021)) and was not taken from another source.
Tropicalmed 06 00027 g001
Tropicalmed 06 00027 g002 550
Figure 2. Dorsal and ventral views of male (A) and female (B) Rhipicephalus microplus; and male (C) and female (D) Rhipicephalus haemaphysaloides.
Figure 2. Dorsal and ventral views of male (A) and female (B) Rhipicephalus microplus; and male (C) and female (D) Rhipicephalus haemaphysaloides.
Tropicalmed 06 00027 g002
Tropicalmed 06 00027 g003 550
Figure 3. Dorsal and ventral views of male (A) and female (B) Haemaphysalis bispinosa; and male (C) and female (D) Haemaphysalis spinigera.
Figure 3. Dorsal and ventral views of male (A) and female (B) Haemaphysalis bispinosa; and male (C) and female (D) Haemaphysalis spinigera.
Tropicalmed 06 00027 g003
Tropicalmed 06 00027 g004 550
Figure 4. Dorsal and ventral views of male (A) and female (B) of Amblyomma testudinarium; and female (C) Ixodes sp.
Figure 4. Dorsal and ventral views of male (A) and female (B) of Amblyomma testudinarium; and female (C) Ixodes sp.
Tropicalmed 06 00027 g004
Tropicalmed 06 00027 g005 550
Figure 5. Habitat suitability maps developed by the best models (a) for Rhipicephalus microplus using land use and land cover (LULC) and digital elevation–shuttle radar topography information (DEM_SRTM) and (b) for Rhipicephalus haemaphysaloides using LULC, Bio 16, and Bio 10.
Figure 5. Habitat suitability maps developed by the best models (a) for Rhipicephalus microplus using land use and land cover (LULC) and digital elevation–shuttle radar topography information (DEM_SRTM) and (b) for Rhipicephalus haemaphysaloides using LULC, Bio 16, and Bio 10.
Tropicalmed 06 00027 g005
Tropicalmed 06 00027 g006 550
Figure 6. Habitat suitability maps developed by the best models (a) for Haemaphysalis bispinosa using LULC, Bio 18, and Bio 16 and (b) for Haemaphysalis spinigera using LULC, Bio 19, and Bio 16.
Figure 6. Habitat suitability maps developed by the best models (a) for Haemaphysalis bispinosa using LULC, Bio 18, and Bio 16 and (b) for Haemaphysalis spinigera using LULC, Bio 19, and Bio 16.
Tropicalmed 06 00027 g006
Table
Table 1. Geographic and climatic characteristics of the study area.
Table 1. Geographic and climatic characteristics of the study area.
CharacteristicTrashigangPema GatshelTotal
Total land (km2)220410233227
Agro-ecological zones (in % of land area)
Wet Subtropical (100–600) *0.120.96.7
Humid Subtropical (600–1200)4.441.716.2
Dry Subtropical (1200–1800)14.828.119
Warm Temperate (1800–2600)29.29.322.9
Cool Temperate (2600–3600)32.7022.3
Alpine (>3600–7500)18.8012.8
Forest coverage (in %)87.681.6
Vegetation typeTemperate - broadleafBroadleaf
ClimateWarm and temperateHot and humid
Annual rainfall (in mm)1115.61649.6
Annual Temperature (min-max in °C)10.3–24.213.1–22.1
Cattle population40,6858252
* Numbers in bracket are the elevation range in masl.
Table
Table 2. Environmental variables used in building MaxEnt habitat suitability models.
Table 2. Environmental variables used in building MaxEnt habitat suitability models.
Variable, UnitCodeSource
Annual mean temperature, °CBio 1WorldClim
Maximum temperature of warmest month, °CBio 5WorldClim
Minimum temperature of coldest month, °CBio 6WorldClim
Mean temperature of wettest quarter, °CBio 8WorldClim
Mean temperature of driest quarter, °CBio 9WorldClim
Mean temperature of warmest quarter, °CBio 10WorldClim
Mean temperature of coldest quarter, °CBio 11WorldClim
Annual precipitation, mmBio 12WorldClim
Precipitation of wettest quarter, mmBio 16WorldClim
Precipitation of driest quarter, mmBio 17WorldClim
Precipitation of warmest quarter, mmBio 18WorldClim
Precipitation of coldest quarter, mmBio 19WorldClim
Elevation (in masl)DEM_SRTMUSGS
Land Use and Land Cover 2016 *LULC 2016NLC, Bhutan
* Categorical variable.
Table
Table 3. Tick species sampled from farm animals in Trashigang and Pema Gatshel districts, Eastern Bhutan (2019).
Table 3. Tick species sampled from farm animals in Trashigang and Pema Gatshel districts, Eastern Bhutan (2019).
Tick SpeciesInfested Cattle (n = 240)Total Ticks (n = 3600)
Rhipicephalus microplus204 (85%, 79.8–89.3) *2530 (70.3%, 68.7–71.7)
Rhipicephalus haemaphysaloides91 (37.9%, 31.7–44.4)677 (18.8%, 17.5–20.1)
Haemaphysalis bispinosa72 (30%, 24.3–36.2)295 (8.2%, 7.3–9.1)
Haemaphysalis spinigera28 (11.7%, 7.9–16.4)90 (2.5%, 2.0–3.0)
Amblyomma testudinarium7 (2.9%, 1.2–5.9)7 (0.19%, 0.07–0.4)
Ixodes sp.11
* Numbers in bracket represents percentage and 95% confidence interval, and n represents the total number.
Table
Table 4. Tick infestation prevalence in cattle in Trashigang and Pema Gatshel districts.
Table 4. Tick infestation prevalence in cattle in Trashigang and Pema Gatshel districts.
Tick SpeciesTrashigang (n = 130)Pema Gatshel (n = 110)Infested Cattle (n = 240)χ2p-Value
Rhipicephalus microplus96 (73.8%, 65.4–81.1) *108 (98.2%, 93.6–99.8)20425.8<0.001
Rhipicephalus haemaphysaloides74 (56.9%, 47.9–65.6) 17 (15.4%, 9.3–23.6)9141.78<0.001
Haemaphysalis bispinosa39 (30%, 22.3–38.6)33 (30%, 21.6–39.5)7201
Haemaphysalis spinigera11 (8.5%, 4.3–14.6)17 (15.4%, 9.3–23.6)282.190.139
* Numbers in bracket represents percentage and 95% confidence interval.
Table
Table 5. Characteristics of the explanatory variables (n = 240) used for the logistic regression analyses.
Table 5. Characteristics of the explanatory variables (n = 240) used for the logistic regression analyses.
VariablesCategoriesTrashigang (n = 130)Pema Gatshel (n = 110)Total (n = 240)%95%CI
Cattle age * Adult1037918275.869.9–81.1
Young27315824.218.9–30.1
Cattle sex Female1118519681.776.2–86.4
Male19254418.313.6–23.8
Cattle breedEuropean738015363.857.3–69.8
Indigenous57308736.230.2–42.7
Altitude ΨMean 16871177
Range 897–2187720–2060
Latitude ϑMean 27.2826.99
Range27.11–27.4126.84–27.12
Longitude ϑMean 91.6191.33
Range91.44–91.7691.09–91.46
* The cattle under calf and heifer were categorized as “young” while the rest were categorized as adults (3 years and above), Ψ altitude in masl, ϑ latitude and longitude in decimal degrees. CI: confidence interval.
Table
Table 6. The final multiple logistic regression models to understand factors associated with infestation prevalence of each tick species (categorized as infested or not) in cattle.
Table 6. The final multiple logistic regression models to understand factors associated with infestation prevalence of each tick species (categorized as infested or not) in cattle.
VariableEstimate ± SEOR (95% CI)p-Value
Rhipicephalus microplus
Intercept506.59 ± 133.11 <0.001
Age (Young) ϖ1.72 ± 0.755.57 (1.5–35.5)
Longitude × 10−0.55 ± 0.150.57 (0.42–0.75)
Rhipicephalus haemaphysaloides
Intercept−191.81 ± 27.31 0.003
Latitude × 100.70 ± 0.102.02 (1.67–2.48)
Haemaphysalis bispinosa
Intercept−2.12 ± 0.55 0.009
Breed (Indigenous) μ0.62 ± 0.291.85 (1.04–3.29)
Altitude/1000.07 ± 0.041.07 (1.002–1.14)
Haemaphysalis spinigera
Intercept0.35 ± 0.74
Breed (Indigenous)1.00 ± 0.432.72 (1.18–6.38)<0.001
Altitude/100−0.21 ± 0.060.81 (0.72–0.90)
ϖ Adult as the referent category, and μ European breed as the referent category. SE: standard error, OR: odds ratio, CI: confidence interval.
Table
Table 7. Variable contribution, permutation importance, corrected Akaike information criteria (AICc), and number of parameters of the best MaxEnt models developed for predicting tick species occurrence in Eastern Bhutan.
Table 7. Variable contribution, permutation importance, corrected Akaike information criteria (AICc), and number of parameters of the best MaxEnt models developed for predicting tick species occurrence in Eastern Bhutan.
ModelVariables% ContributionPermutation ImportanceAICcParameters
Rhipicephalus microplusDEM (Elevation)56.176.43084.3515
LULC (Land use and land cover)43.923.6
Rhipicephalus haemaphysaloidesLULC (Land use and land cover)469.41360.4713
Bio 16 (Precipitation of wettest quarter)43.129.6
Bio 10 (Mean temperature of warmest quarter)10.961
Haemaphysalis bispinosaLULC (Land use and land cover)63.239.31095.3213
Bio 18 (Precipitation of warmest quarter) 30.644.7
Bio 16 (Precipitation of wettest quarter)6.216.1
Haemaphysalis spinigeraBio 16 (Precipitation of wettest quarter)45.726.3433.918
Bio 19 (Precipitation of coldest quarter)33.466.9
LULC (Land use and land cover)20.86.8
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Namgyal, J.; Lysyk, T.J.; Couloigner, I.; Checkley, S.; Gurung, R.B.; Tenzin, T.; Dorjee, S.; Cork, S.C. Identification, Distribution, and Habitat Suitability Models of Ixodid Tick Species in Cattle in Eastern Bhutan. Trop. Med. Infect. Dis. 2021, 6, 27. https://0-doi-org.brum.beds.ac.uk/10.3390/tropicalmed6010027

AMA Style

Namgyal J, Lysyk TJ, Couloigner I, Checkley S, Gurung RB, Tenzin T, Dorjee S, Cork SC. Identification, Distribution, and Habitat Suitability Models of Ixodid Tick Species in Cattle in Eastern Bhutan. Tropical Medicine and Infectious Disease. 2021; 6(1):27. https://0-doi-org.brum.beds.ac.uk/10.3390/tropicalmed6010027

Chicago/Turabian Style

Namgyal, Jamyang, Tim J. Lysyk, Isabelle Couloigner, Sylvia Checkley, Ratna B. Gurung, Tenzin Tenzin, Sithar Dorjee, and Susan C. Cork. 2021. "Identification, Distribution, and Habitat Suitability Models of Ixodid Tick Species in Cattle in Eastern Bhutan" Tropical Medicine and Infectious Disease 6, no. 1: 27. https://0-doi-org.brum.beds.ac.uk/10.3390/tropicalmed6010027

Article Metrics

Back to TopTop