Molecular survey and genetic diversity of piroplasmids in equids from Midwestern Brazil

We evaluated the distribution of piroplasmids in equids from the Mato Grosso state in Midwestern Brazil using molecular methods and the interspecific genetic diversity. For this, 1,624 blood samples of equids from 973 farms were examined by PCR, using primer pairs that amplify a fragment of the genes rap-1 and ema-1 of Babesia caballi and Theileria equi, respectively. For molecular characterization and phylogenetic studies, 13 and 60 sequences of the rap-1 and ema-1 genes, respectively, were used to build a dendogram using maximum parsimony. B. caballi and T. equi were detected in 4.11% and 28.16% of the farms, respectively, and molecular prevalence was 2.74% for B. caballi and 25.91% for T. equi. The location of the farms and animals raised in the Pantanal ecoregion influence the probability of equids testing positive for B. caballi and T. equi. Moreover, age and herd purpose were variables significantly associated with T. equi infection. The sequences of B. caballi presented 1.95% intraspecific variability, contrasting with 2.99% in T. equi. Dendrograms for both species demonstrated the presence of subgroups with high values of support of branches. However, it is not possible to associate these groups with geographic origin and/or ecoregion.


Introduction
Equine piroplasmosis is a tick-borne disease of horses and other Equidae (donkeys, mules, and zebras) caused by the hemoprotozoan parasites Theileria equi and Babesia caballi (WISE et al., 2013;WISE et al., 2014).These two parasites are classified within the phylum Apicomplexa, and despite their biological differences, they cause a similar pathology and have similar life cycles and vector relationships.Acute disease is characterized by fever, malaise and reduced appetite, increased pulse rate and respiration, anorexia, constipation followed by diarrhea, tachycardia, petechiae, splenomegaly, thrombocytopenia, and hemolytic anemia leading to hemoglobinuria and icterus (MEHLHORN & SCHEIN, 1998;UILENBERG, 2006;ROTHSCHILD, 2013;WISE et al., 2013;HABIBI et al., 2016).
Horses that recover from acute disease remain chronically-infected carriers without any overt signs of disease and can be reservoirs for the transmission of these protozoan pathogens with ticks as a major vector.Parasitemia is often too low to detect on blood smears, but infected animals can be identified by serology or polymerase chain reaction (PCR) (WISE et al., 2013;WISE et al., 2014;SCOLES & UETI, 2015).Ticks are the definitive hosts for both T. equi and B. caballi, and in Brazil, these two parasites are mainly transmitted by Rhipicephalus (Boophilus) microplus and Dermacentor nitens, respectively (SCOLES & UETI 2015).Furthermore, equine piroplasmosis has an economic importance because of the potential effect of the international movement of these horses for trade and for participation in international sporting events (WISE et al., 2013;WISE et al., 2014;SCOLES & UETI, 2015).
Piroplasmosis occurs in most countries around the world and infection is maintained within equine populations as long as competent vectors are present.According to the World Organization for Animal Health (OIE), Central and South America, Cuba, Africa, Asia, the Middle East, and Southern Europe are considered to be endemic regions.Within South America, the disease is readily identified in all regions, except for the southernmost areas of Chile and Argentina (WISE et al., 2013).In Brazil, previous studies using serological methods have shown that equine piroplasmosis is endemic, with prevalence rates varying from 19.4% to 100% (KERBER et al., 1999;HEIM et al., 2007;KERBER et al., 2009;BALDANI et al., 2010;SANTOS et al., 2011;MACHADO et al., 2012;TORRES et al., 2012;VIEIRA et al., 2013;PROCHNO et al., 2014;FERREIRA et al., 2016;GUIMARÃES et al., 2016;BRAGA et al., 2017).
Despite the endemicity of piroplasmosis among equids in Brazil, studies involving the molecular detection of the parasite in carrier animals (HEIM et al., 2007;BALDANI et al., 2010;PECKLE et al., 2013;MACHADO et al., 2012;BARROS et al., 2015;FERREIRA et al., 2016;BRAGA et al., 2017) that have also examined genetic diversity (BRAGA et al., 2017;PECKLE et al., 2018) are scarce.The aim of this study is to survey the distribution of piroplasmid (T.equi and B. caballi) parasites using molecular methods in equids from the Mato Grosso state and evaluate the interspecific genetic diversity.

Study area
With nearly one thousand square kilometers, the state of Mato Grosso is located in Midwestern Brazil, encompassing biotic elements from Amazonia, Cerrado (Brazilian savannah), and Pantanal ecoregions.Southern Amazon forests correspond to nearly 53% of the state and are located in its northern portion.The Cerrado -a huge tropical savannah that extends across the Brazilian Central Plateau -occupies nearly 40% of the state, in its central portion; and the Pantanal -the largest tropical wetland in the world -occupies another 7% of the state area, in its southwestern portion (BRASIL, 2014).
Blood samples were collected during the course of a previous study on the epidemiology of equine infectious anemia (EIA) within the three ecoregions in the state of Mato Grosso.This study was conducted from September to December 2014 (BARROS, 2017).The stored samples were made available for this study after the conclusion of the previous study.
Concomitant with sampling, each farmer was asked to fill out a questionnaire for subsequent epidemiological data analysis.

Sampling procedures
According to the Mato Grosso Institute for Agricultural Defense (INDEA/MT, 2014), there are currently 306,931 equids on 53,710 farms distributed through three ecoregions (Cerrado, Pantanal and Amazonia) in the state, as observed in table 1.
The sample size was determined by the statistical formula (SCHAEFFER et al., 2011).This formula considered population size (306,931 animals); 50% estimated prevalence (due the absence of prior studies and the maximum number of occurrences stay around 50% on the normal distribution); 3% maximum error (d); 95% confidence interval (CI); and a design effect of 1.45 (deff).The sample size was adjusted from 1,542 to 1,624 animals (1,389 equines, 227 donkeys and eight asinines) for the possibility of loss, and then distributed among the three ecoregions (Table 1).Considering that 77% of the farms had no more than five equids, the following criteria were used: only one animal was sampled on farms with five or fewer equids, while five animals were randomly sampled on farms that had more than five equids.Sampling was performed in two stages: selection of farms and selection of Equidae.Selection procedures were as follows: a pre-established number of random farms (primary sampling units) was selected, which resulted in the sampling of 973 farms (Table 1).Next, a random selection of equids on each farm (secondary sampling units) was performed.Figure 1 shows the distribution of the farms sampled for the screening of B. caballi and T. equi infection in the three ecoregions (Amazonia, Pantanal and Cerrado) of the state of Mato Grosso.

DNA extraction and PCR
DNA was extracted from 60 μg of each blood clot sample by the phenol-chloroform-isoamyl method as described by Sambrook & Russel (2001).In order to verify the presence of amplifiable DNA, 163 randomly selected samples were submitted to internal control PCR assays targeting fragments of mammalian glyceraldehyde-3-phosphate dehydrogenase (GAPDH), as previously described (BIRKENHEUER et al., 2003).Each sample of extracted DNA was submitted for PCR analysis, using species-specific primer pairs based on the genes rap-1 for B. caballi and ema-1 for T. equi, which amplified DNA fragments of 570 and 750 base pairs (bp), respectively.This PCR analysis procedure was adapted from Alhassan et al. (2005).The reaction mixture consisted of 1.5 μl template DNA (B.caballi) or 2 μl (T.equi), 2.5 μl 10× PCR buffer, 4.0 μl MdNTP mix, 0.15 μl Taq DNA polymerase (Sigma-Aldrich), 0.5 μl each of 10 pmol B. caballi and T. equi specific primers, and 15.35 μl of distilled water (B.caballi) or 15.85 μl (T.equi), for a total volume of 25 μl.The amplification conditions for B. caballi and T. equi were as follows: 40 cycles with an initial enzyme activation at 96 °C for 10 min, denaturation at 96 °C for 1 min, primer annealing at 60.5 °C for 1 min, extension at 72 °C for 1 min and final extension at 72 °C for 10 min.For all PCR assay, B. caballi and T. equi DNA samples obtained from naturally infected horse (BARROS et al., 2015) and ultra-pure sterile water were used as positive and negative controls, respectively.PCR product sizes were analyzed on 1.5% agarose gel stained with GelRed™ Nucleic Acid Gel Stain, and visualized in a ChemiDoc XRS system.

Sequencing and phylogenetic analysis
For molecular characterization and phylogenetic studies of B. caballi and T. equi isolates, fragments of 569 and 675 bp from rap-1 and ema-1 genes, respectively, were used.DNA from positive samples was purified and sequenced in an automated sequencer (ABI-PRISM 3500 Genetic Analyzer, Foster City, CA).
All the sequences obtained were aligned with sequences previously determined for B. caballi and T. equi, available in GenBank, using ClustalX (THOMPSON et al., 1997).The sequences were adjusted manually using GeneDoc (NICHOLAS et al., 1997).
The rap-1 and ema-1 sequences were used to build a dendogram using maximum parsimony, as implemented in PAUP version 4.0b10 (SWOFFORD, 2002), with 500 bootstrap replicates, random stepwise addition starting trees (with random addition sequences), and TBR branch swapping.

Prevalence and statistical analysis
Piroplasmid prevalence on each farm was analyzed based on the logic that a farm would be considered positive if at least one animal tested positive.A multiple logistic regression model was created in several steps as follows: 1.For continuous variables, a linearity test was performed using graphs for log odds.When the visual evaluation indicated a nonlinear trend, the continuous variable was categorized.The optimization methodology of the software package CatPredi R was used (R DEVELOPMENT CORE TEAM, 2015).2. For the categorical variables, a chi-square or exact Fisher test was performed on the results and various predictor variables.All variables with a value of p <0.20 were incorporated into the model.3. The selected variables were tested for collinearity including a multicollinearity analysis according to Vatcheva et al. (2016), to ensure a mean variance inflation factor (VIF) <10 prior to applying them in the models.4. The model was constructed by inserting all the variables selected in the previous steps.The less-significant variables (according to the Wald test) were then removed and the logistic regression analysis was then repeated.This process was repeated, and the adjusted model was compared to the previous one by a likelihood ratio test to check for confounding factors.When there was a change of more than 30% in the parameter estimates, the removed variable was considered to be a confounding factor and again included in the model.
The variables used to construct the herd model were ecoregion (biome), type of holding (farm, breeding facility or horseback riding facility), herd purpose (service, sport/recreation or reproduction), type of management (extensive, intensive or semi-extensive), needle/syringe sharing among animals, flooded area at the farm, veterinary assistance, and animal density (greater than 5 hectares per animal or less than 5 hectares per animal).
The animals were analyzed following the same steps used to analyze the farms.However, a complex analysis, according to Korn & Graubard (1999), was applied to the tests to correct the effect of the weighted values used for the selection of animals.The following formula was used to determine the weighted value of each animal: (# total animals/# total animals sampled) x (# total animals on the farm/# total animals sampled on the farm).
The model was built using the following variables: ecoregion (biome), animal species, sex, age (classified as greater than six months and less than 3 years, greater than 3 years and less than 10 years, or greater ten than years), and animal purpose (work or sports/ recreation).Statistical analysis was performed using the R statistical software package (R DEVELOPMENT CORE TEAM, 2015).
This study was approved by the Animal Research Ethics Committee of the Federal University of Mato Grosso, under reference no.23108.016503/14-9.
In the multivariate analysis, after the model was applied, location (ecoregion) was the only variable found to be significantly associated with the prevalence of babesiosis on farms.Farms in the Pantanal ecoregion were 73% less likely to have equine infected with B. caballi.Table 2 describes the relationship between the positive farms and the variable of the logistic model.Furthermore, animals raised in the Pantanal ecoregion were 3% less likely to be infected with B. caballi (Table 3).
Two variables were found to be significantly associated with the prevalence of theileriosis on farms: location (ecoregion) in the state, and herd purpose.Table 4 describes the relationship between the positive farms and the variables of the logistic model.Farms in the Pantanal ecoregion were 57% less likely to have equine infected with T. equi, and farms that used equids to perform service were 84% less likely to have T. equi infected animals.Age and location of equids were variables significantly associated with T. equi infection.Equids between six months and three years old were 1.2 times (20%), more likely to be infected with T. equi, and equines older than 10 years were 10% less likely to have T. equi infection.Furthermore, animals raised in the Pantanal ecoregion were 20% less likely to be infected with T. equi.Table 5 describes the relationship between the prevalence of infection by T. equi in the Equidae population and the variables of the logistic model of the 1,624 animals sampled.
Thirteen sequences of B. caballi were obtained from 13 farms and 60 sequences of T. equi were obtained from 52 properties.The sequences of B. caballi presented 1.95% intraspecific variability, contrasting with 2.99% in T. equi.Dendrograms based on maximum parsimony for B. caballi (Figure 2) and T. equi (Figure 3) demonstrated the presence of subgroups with high support values of branches (Figure 3).The GenBank nucleotide sequence accession numbers generated in the present study are: MG906574-MG906586 for partial sequences of the rap-1 gene of B. caballi and MG906587-MG906646 for partial sequences of the ema-1 gene of T. equi.

Discussion
In this study, we investigated the occurrence of piroplasmids (T.equi and B. caballi) using molecular methods in 1,624 equids (1,389 equine, 227 donkeys, eight asinines) from 973 farms from an extensive region of Midwestern Brazil, which covers three different ecoregions (Amazonia, Cerrado and Pantanal) of the Mato Grosso state.Furthermore, we evaluated variables (risk factor) thought to be associated with the occurrence of infection in animals and Equidae herds, as well as the interspecific genetic diversity of B. caballi and T. equi.
Despite the fact that B. caballi and T. equi have previously been detected in the Pantanal ecoregion (BARROS et al., 2015), this is the first study, to our knowledge, to show the prevalence of piroplasmids in Equidae in the entire state.Herein, we confirm that infection is widespread in the studied region, with higher prevalence of T. equi (25.91%) than B. caballi (2.74%) among equids.These findings are similar to those found in other studies conducted in Brazil (HEIM et al., 2007;VIEIRA et al., 2018) and other countries (POSADA-GUZMÁN et al., 2015;GUVEN et al., 2017).This may reflect differences in the biology of the vectors, as well as differences in persistence of the parasite in the equine    Unexpectedly, when the purpose of the herd was for service, the likelihood of a farm to be positive for T. equi decreased by 83.02%, once horse's activities related to work have already been described as factors that affect the likelihood of infection by T. equi (PECKLE et al., 2013;ABUTARBUSH et al., 2012).Furthermore, the activities that the horse exerts on the farm, mainly, cattle handling, allows move into large pasture perimeters and may increase their chances of infestation by potential tick vector (PECKLE et al., 2013).
In the present study, farms located in the Pantanal are less likely to have animals infected with B. caballi and T. equi.Moreover, animals raised in the Pantanal ecoregion were less likely to be infected with both protozoans (B.caballi and T. equi).In fact, the lowest prevalence among farms and the Equidae population, both B. caballi and T. equi, were observed in the Pantanal ecoregion.
The Pantanal ecosystem is characterized by flood and drought cycles, where over 80% of the plains are submerged during the rainy period (ALHO, 2008).Therefore, the total or partial flooding of pastures may eliminate ticks from the vegetation, and could also affect tick survival.Ramos et al. (2016) have studied ecological relationships between Nellore cattle (Bos indicus) and ticks within the Brazilian Pantanal, and although R. (B.) microplus (the known vector of T. equi) was the most abundant tick species on cattle (primary tick-host), the overall highest tick prevalence on bovines in the dry season was that of Amblyomma sculptum nymphs.Additionally, during the entire study, R. (B.) microplus was not found in the environment.
According to Campbell & Glines (1979), conditions of greater or lesser flooding of land have a great influence on the biology of Ixodidae, in addition to other environmental factors such as ambient temperature and humidity.Louzada & Daemon (2003) and Santos et al. (2012) have verified the effects of the immersion of engorged females of R. (B.) microplus, and observed a negative effect on reproduction parameters, through a negative effect over the weight of the eggs, alterations in the egg-laying and in the hatching.Moreover, Paula et al. (2004) have observed that immersion in water is deleterious to the free-living phase of D. nitens larvae (the known vector of B. caballi).Thus, the differences between ecoregions and their edaphoclimatic conditions should influence the epidemiology of equine piroplasmosis.
T. equi seroprevalence rates increase with age, with the older group related to a chronic infectious status (DE WAAL, 1992).Santos et al. (2011) have observed in 714 equids (702 horses and 12 mules) from the Rio de Janeiro state, that animals aged six to 12 months showed 87.5% of T. equi seropositivity.Furthermore, considering the absence of maternal antibodies in this age group, they have suggested that the first T. equi infection occurred between six and 12 months of age or possibly even younger since there is no passive immunity in animals over six months old.We also observed by molecular methods that equines older than six months and less than three years have a 1.2-fold increased risk of T. equi infection, however animals older than 10 years have a 10.08% decrease in the chance of being infected with T. equi.Our results are in agreement with Bartolomé Del Pino et al. (2016), who also showed prevalence for T. equi decreasing with age as detected by PCR.However, they are contrary to some reports in literature that have highlighted a cumulative age-dependent increase in the PCR-prevalence of T. equi (RÜEGG et al., 2008) and to others that have shown no association with age and infection rates (PECKLE et al., 2013).Hence, the circumstances that influence positivity rates of T. equi during the lifetime of a horse, require further study due to the sensitivity of the various methods used in studies involving host-parasite interactions (BARTOLOMÉ DEL PINO et al., 2016).Indeed, the associations of age and infection that are reported in the literature may be due to the small sample size of animals in some categories, which could possibly interfere with the analysis of the data (SANTOS et al., 2011).
Characterizations of piroplasm sequences from Brazilian horses are scarce.Intraspecific analysis of the sequences obtained from equids of different ecoregions in the state of Mato Grosso showed 1.95% sequence divergence for B. caballi and 2.99% for T. equi.Studies based on ribosomal (18S) genes demonstrated variability between 0.5% and 2.1% for T. equi, allowing the segregation into different genotypes (BRAGA et al., 2017;PECKLE et al., 2018).
The intraspecific divergence of T. equi found within the ema-1 gene was greater than that presented by the 18S gene.However, it was not possible to distinguish between subgroups/genotypes related to the different ecogeographic regions studied.The absence of clear patterns related to the geographic region and/or ecoregion may be due to the internal migrations of the animals, making it impossible to form the barriers or geographic isolation needed to identify clusters.
Further studies with the generation of new sequences of B. caballi and T. equi from different regions and Brazilian biomes could minimize the effects of animal traffic in the analysis of intraspecific variability.

Conclusions
Equine piroplasmosis occurs throughout the three ecoregions of the Mato Grosso state in Midwestern Brazil, with the Amazonia ecoregion containing the largest number of animals and farms infected by both agents (B.caballi and T. equi), while farms and animals raised in the Pantanal ecoregion were less likely to be infected with both protozoans.Animals between six months and two years and eleven months old are at an increased risk of T. equi infection.Interestingly, animals over the age of ten years were shown to be at reduced risk of T. equi infection.Intraspecific analysis of the sequences obtained for B. caballi and T. equi did not allow us to distinguish between subgroups/genotypes related to the different ecogeographic regions under area studied.

Figure 1 .
Figure 1.Distribution of farms sampled for the detection of Babesia caballi and Theileria equi in equids from the three ecoregions (Pantanal, Cerrado and Amazonia) of the Mato Grosso state, Brazil, between September and December 2014.

*
Animal age, classified as greater than six months and less than 3 years, greater than 3 years and less than 10 years and greater than ten years.

Figure 2 .
Figure 2. Maximum parsimony tree inferred from rap-1 gene sequences of Babesia caballi.Numbers at nodes are the support values for the major branches (bootstrap over 500 replicates).

Figure 3 .
Figure 3. Maximum parsimony tree inferred from ema-1 gene sequences of Theileria equi.Numbers at nodes are the support values for the major branches (bootstrap over 500 replicates).

Table 1 .
Columns left to right: (1) Ecoregion in Mato Grosso state, Brazil.(2) Total number of farms in each ecoregion.(3) Number of farms sampled in the current survey.(4) Total number of equids in each ecoregion.(5) Number of equids sampled in the current survey.All the data were collected between September and December 2014.

Table 2 .
Results of the logistic regression model for risk factors associated with Babesia caballi herd prevalence on the 973 farms sampled in the state of Mato Grosso, Brazil, between September and December 2014.

Table 5 .
Results of the logistical regression model for risk factors associated with the prevalence of infection by Theileria equi in individual equids from the state of Mato Grosso, Brazil, between September and December 2014.

Table 3 .
Results of the logistical regression model for risk factors associated with the prevalence of infection by Babesia caballi in individual equids from the state of Mato Grosso, Brazil, between September and December 2014.

Table 4 .
Results of the logistic regression model for risk factors associated with Theileria equi herd prevalence on the 973 farms sampled in the state of Mato Grosso, Brazil, between September and December 2014.