Climate change implications for the distribution of the babesiosis and anaplasmosis tick vector, Rhipicephalus (Boophilus) microplus

Climate change ranks among the most important issues globally, affecting geographic distributions of vectors and pathogens, and inducing losses in livestock production among many other damaging effects. We characterized the potential geographic distribution of the ticks Rhipicephalus (Boophilus) microplus, an important vector of babesiosis and anaplasmosis globally. We evaluated potential geographic shifts in suitability patterns for this species in two periods (2050 and 2070) and under two emissions scenarios (RCPs 4.5 and 8.5). Our results anticipate increases in suitability worldwide, particularly in the highest production areas for cattle. The Indo-Malayan region resulted in the highest cattle exposure under both climate change projections (2050), with increases in suitability of > 30%. This study illustrates how ecological niche modeling can be used to explore probable effects of climate change on disease vectors, and the possible consequences on economic dimensions.


Introduction
Rhipicephalus (Boophilus) microplus is the most important tick in transmission of bovine parasitic diseases around the world [1]. The principal hosts for this species are cattle, but interactions have been shown with buffalo, horses, donkeys, dogs, deer, sheep, and goats [2][3][4]. High incidence of this tick is associated with economic losses, particularly in cattle [5]. This tick is responsible for transmission of the protozoa Babesia bovis and B. bigemina, and the bacterium Anaplasma marginale, which are the main pathogens of bovine babesiosis and anaplasmosis, respectively [6][7][8]. These diseases induce extreme emaciation in livestock, culminating in death.
Losses of US$13-18 billion are caused by these pathogens globally each year [9,10]. Developing countries have seen the strongest consequences caused by R. microplus (e.g. in South America; [5,11]). Such countries generally lack effective control mechanisms for the tick, so that economic losses are exacerbated and livestock production is reduced markedly [12]. Estimates regarding milk losses caused by this tick species are 90.2 L per cow yearly; this reduction and losses in milk products generate production drops of around US$922 million yearly [13]. Cattle at high risk regarding this tick are valued at US$3 billion annually in Brazil alone [5], one of the most important developing countries in terms of livestock, responsible for exporting ~ 43 million tons of milk and meat [14].
Population growth and establishment are related to the following: (1) historical contingencies and geographic barriers [15]; (2) biological factors such as host availability, competition; and (3) environmental conditions such as temperature and humidity [16,17]. Future changes in climate include modifications of temperature and precipitation regimes; these environmental factors are crucial in delimiting species distributions and determining the success of population establishment [18]. Changes in these factors may modify the life cycle, abundance, and distribution potential for R. microplus [19,20] and other disease vectors or species with importance in human and animal health [21][22][23] However, many based-ENM studies of related disease arthropods have not used environmental variables related directly to the physiology of these species instead relying on general climate datasets that are easily available. In this paper, we assess possible potential areas for R. microplus under present-day and future climate conditions for two greenhouse gas emissions scenarios (RCP 4.5 and RCP 8.5) in two time periods (2050 and 2070), including a variable known to be crucial to the physiology of this species. We evaluate potential hotspots for this species and their coincidence with livestock concentrations to determine the most important risk areas under climate change for bovine parasitic diseases transmitted by this tick species.

Occurrence data
Occurrences were obtained from two different sources: (1) published data based on searches of different databases (Web of Science, Scopus, and Google Scholar); we used the keywords "Rhipicephalus", "Rhipicephalus microplus", "Boophilus microplus" and "Rhipicephalus (B.). microplus". We also obtained (2)  Data were collected for the period 1970-2018. Data lacking georeferencing (obtained mainly from published papers) were assigned coordinates via searches in Google Earth. We reduced biasing effects of spatial autocorrelation in occurrence data using a distance filter of 22 km in the spThin R package [38]. We chose a random 50% of the occurrence data for calibrating models, and used the remaining 50% to evaluate the models. Our initial 1487 occurrences for R. microplus in America reduced to 531 with spatial filtering (Figure 1). We also considered occurrence data for this species from Africa [39, 40]; these 145 African occurrences were used as independent evaluation data, with the same spatial filter, for an additional model evaluation.

Environmental data
We used 15 bioclimatic variables from WorldClim version 1.4 (Table 1) [41], excluding four variables known to include spatial artefacts [42]. WorldClim variables were derived from climatic data for 1950-2000. To add a variable known to be important to the physiology of this species [43], we obtained relative humidity (RH) from the Coupled Model Intercomparison Project [44], which we downscaled by the Delta Method, commonly used in climate data [45]. To summarize future conditions, we used outputs from 20 general circulation models (GCM) available from Climate Change, Agriculture and Food Security [46] (Additional file 1). We used two greenhouse gas emissions scenarios (RCP 4.5 and RCP 8.5) for two time periods (2050 and 2070) to explore model-to-model variation. The environmental data were used at a spatial resolution of 0.2° (~ 22 km) under both present-day and future conditions. Dimensionality was reduced by calculating Pearson correlations over the entire calibration area, removing one from each pair of variables with correlations ≥ 0.80. Uncorrelated variables were according to different variable sets: we used all possible variable combinations (120 sets) for model calibration and evaluation [47]. Variables used as candidates for inclusion in models therefore included annual mean temperature, temperature seasonality, minimum temperature during the coldest month, annual precipitation, precipitation during the driest month, precipitation seasonality, and relative humidity.

Model calibration and evaluation
Models were calibrated with the kuenm R package, running [48] Maxent 3.4.1 [49] and using model selection approaches [50]. We used significance, performance (omission rate), and model complexity to choose optimal parameter settings from among candidate models. All possible combinations of linear (l), quadratic (q), product (p), threshold (t), and hinge (h) feature types were tested, as were different regularization multiplier values (0.1, 0.3, 0.5, 0.7, 1, 2, 3, 5, 7, and 10). Models are built using all possible combinations of the seven environmental variables. Hence, we explored a total of 5400 candidate models. Significance testing was via partial ROC [51], with acceptable omission error of E = 0.05. Finally, we evaluated the model complexity using the Akaike information criterion with a correction for small sample size (AICc), via a code derived from Warren et al. [50].
We used the accessible area concept as a means of choosing the area over which to calibrate our models [15,52], using a buffer of 200 km around occurrence data as a proxy to this area ( Figure 1). Final models were summarized as the median of 10 bootstrap replicates of the model corresponding to the best model parameter set, and transferred to future conditions worldwide. We used the kuenm package [48] both to evaluate final models and to transfer models to future conditions. Areas presenting extrapolative conditions were identified with a MOP analysis [53], comparing conditions between calibration and transfer areas across the 20 GCM × 2 RCPs × 2 time periods that made up our future scenarios.
Model transfers were summarized and simplified via a binarization process. Median model transfers were binarized using an acceptable omission rate of E = 0.05. Binary maps were used to determine climate uncertainty in the models, which we summarized as agreement among multiple scenarios between present and future to determine areas of stability. We used a threshold of ≥ 60% (12) of agreement among GCM as a relatively clear signal of presence or absence of suitable conditions. Finally, to provide an additional model evaluation, we used occurrence data from the African range as an evaluation of the model's predictive ability. We obtained cattle abundance data from the Food and Agriculture Organization of the United Nations (FAO) [54] and Robinson et al. [55] to evaluate implications of future changes in the range of the tick for cattle production. We used different categories (cattle abundances: Additional file 2) to evaluate possible impacts of the tick under future projections. We evaluated the highest potential range expansion of the tick with regards to the different world cattle production areas.

Results
We evaluated 5400 candidate models; 3687 of these models were statistically significant (P < 0.05), of which 1348 showed good performance (i.e., omission error ≤ 0.05); however, a single model was selected on the basis of low complexity (AICc = 10,398.67), in that the difference in AICc between it and the next best model was large (26.4786). The best model included all feature types (lqpth: Linear, Quadratic, Product, Threshold and Hinge) with a regularization parameter of 1. Variables selected were annual mean temperature, temperature seasonality, annual precipitation, precipitation of driest month, and relative humidity ( Table 1). The biggest relative contribution of environmental variables to the model was from annual mean temperature (26.6%), whereas the smallest contribution was from annual precipitation (14.6%; Table 1). Our independent evaluation (i.e., predicting the African distribution) was significantly different from random predictions according to the pROC evaluation (P < 0.001). Omission rate was 7.5%, with only 11 failures out of 145 evaluation points.
In the Americas, the present-day model shows high suitability for R. microplus in North and South America: in Brazil (central, western, southern), Uruguay (northern), Argentina (northern, eastern), and across Central America, Mexico, and the southern USA ( Figure 2). Models also indicate high suitability across much of sub-Saharan Africa, except for the interior of South Africa and Botswana. Western Europe, Southeast Asia, and coastal parts of Australia also had high suitability (Figure 2 [44] were used for modeling. Five variables a were selected for modelling via our model selection processing [47]. b Variables excluded because they have unrealistic spatial artefacts [42]. were in much the same regions. The areas that presented extrapolation risk from the calibration area are shown in Figure 3.

Acronym Description % variable contribution
Our models show high suitability in the present day, and increases in suitability in the future, in places with the highest abundances of cattle around the world (Tables 2 and 3, Additional file 2). Low cattle abundances (0 to 1 individuals/10 km 2 ) were the areas most likely to see increases in suitability for the tick, with a possible increase of ~ 21% in the world for this category (Tables 2 and 3); the Palearctic, Neotropical, and Indo-Malayan regions were those most likely to see increases in suitability. Highest cattle abundances (> 100 individuals/10 km 2 ) in the Indo-Malayan region were projected to see suitability increases of ~ 34% in 2050 and ~ 16% in 2070 (Tables 2 and 3). All study regions show increases in suitability (~ 1% to ~ 134%); major changes were noted in the Nearctic Region (Tables 2 and 3).

Discussion
Babesiosis and anaplasmosis may be related to the same environmental factors as R. microplus because they depend on this tick as a main vector [6][7][8]. Different wild animals serve as hosts for this tick; the best-studied is the white-tailed deer Odocoileus virginianus [3]: livestock reservoirs of these pathogens include goats [4]. Host species are important for dispersal, although this tick is not a specialized ectoparasite; this generalist habit facilitates dispersal by wild and livestock animals [6][7][8].
The potential distribution of R. microplus is related to a diverse suite of ecological and environmental factors around the world [1,17,19,20,[56][57][58][59][60]. Moreover, in order for the tick populations to spread in a region, individuals must first be introduced, then go through a phase of adaptation to the local hosts and then reach a population density that allows mating and reproductive success of adults [1,31,58]. Complex relationships with other species, especially wild hosts, are particularly relevant for individual dispersal and population occurrence in suitable environments [3,61]. Our models show suitable areas in Europe and Asia, places without species present; however, the environmental conditions correspond to species establishment. Climate factors affect the ticks life cycle and geographic distribution [62,63]. The most important factors identified in our model construction were annual mean temperature, precipitation seasonality, and relative humidity [17].
Climate seasonality is an important factor in the R. microplus life cycle; variation in this factor influence the number of generations (three to four per year), increasing the population size and potentially facilitating dispersal [62,63]. We found important environmental variables similar to those identified in field analyses of this species [17]. Annual mean temperature, seasonality in precipitation, and variables derived from humidity were crucial to these models [17]; however, no previous study has used relative humidity in model construction, which is a fundamental variable affecting this species' development [64,65].
Currently, R. microplus does not occur in some regions that our models signaled as suitable, which can be explained by several factors; for example: (1) some trophic activities or other biological interactions that reduce the possibility of occurrence in Europe, Australia, and parts of northern Asia. (2) Historical conditions in particular, major dispersal barriers, decrease the access of this species to Europe, Australia, and northern Asia. In addition, for the species to be introduced into these areas, it is important to consider the number of individuals required to overcome demographic limitations and produce permanent populations, sufficient host density to support tick populations, and immune responses of the host [17].
Economic losses in cattle are generated by tick infestations in cattle herds around the world. Countries with high cattle populations have experienced significant losses in meat production (e.g., Brazil, with US$3.4 billion per year [5]; Tanzania, with US$364 million per year [66]; and Mexico, with US$573 million per year [67]). Thirtyeight percent of the world's cattle population is located in India; babesiosis is a common disease there, driving important economic losses [68]; indeed, ~ 4% of cattle in northeastern India died from this disease [69]. Our results anticipate increases in suitability for the tick, particularly in northern India, potentially increasing losses under future climates (Tables 2 and 3 (Tables 2 and 3, Additional file 2). Increases in suitability for the tick in West Africa (Côte d'Ivoire, Benin), increases cattle exposure to babesiosis and anaplasmosis [31]; of particular note is that importation of Brazilian cattle in these countries creates a situation optimal for tick establishment a situation confirmed by independent models from De Clercq et al. [31].
Changes in the potential distribution of this tick in relation to climate change have been discussed and documented in previous studies [17,19,56]. However, our work used an uncertainty evaluation that explored 20 GCM, two greenhouse gas emission scenarios, and two time periods (Additional files 3 and 4): we also used independent data for model testing (from Africa) and incorporated novel data on relative humidity. Our model is strongly consistent and accurate (Additional files 3, 4, 5), but has two strong limitations. (1) We did not include biotic factors (interactions) in the model, especially because this species requires a host; however, given the broad host range of this tick, we did not have any logical means of assessing all the different possible hosts in the world. (2) We did not consider relationships between numbers of ticks and numbers of cattle in different regions of the world. However, our models provide a view of suitability for the tick under future conditions, which we translate into metrics of cattle exposure for different cattle populations worldwide. The results were consistent with future potential increases of this invasive species around the world, despite the biotic factors not evaluated in our models (e.g. vegetation, host availability and biological interactions).