New insights from molecular characterization of the tick Rhipicephalus (Boophilus) microplus in Brazil

The Rhipicephalus (Boophilus) microplus complex currently consists of five taxa, namely R. australis , R. annulatus , R. (B.) microplus clade A sensu , R. microplus clade B sensu , and R. (B.) microplus clade C sensu . Mitochondrial DNA-based methods help taxonomists when they are facing the morpho-taxonomic problem of distinguishing members of the R. (B.) microplus complex. The purpose of this study was to perform molecular characterization of ticks in all five regions of Brazil and infer their phylogenetic relationships. Molecular analysis characterized 10 haplotypes of the COX-1 gene. Molecular network analysis revealed that haplotype H-2 was the most dispersed of the studied populations (n = 11). Haplotype H-3 (n = 2) had the greatest genetic differentiation when compared to other Brazilian populations. A Bayesian phylogenetic tree of the COX-1 gene obtained strong support. In addition, it was observed that the population of R. (B.) microplus haplotype H-3 exhibited diverging branches among the other Brazilian populations in the study. The study concludes that the different regions of Brazil have R. (B.) microplus tick populations with distinct haplotypes.

In most cattle-producing countries, control methods against this tick are costly and require the use of acaricides. However, acaricide-resistant R. (B.) microplus populations have become a worldwide problem, and molecular ecology studies and new tick control technologies are now required to preserve cattle production (RODRÍGUEZ-VIVAS et al., 2007;ANDREOTTI et al., 2011;GUERRERO et al., 2014).
Although the history of the dissemination of R. (B.) microplus is not well documented, the species is known to have originated in India (HOOGSTRAL, 1986). According to Barré & Uilenberg (2010), R. (B.) microplus originated in the southern and southeastern regions of Asia and was spread throughout the tropical and subtropical belts via cattle, arriving in Brazil between the 16th and 17th centuries.
Mitochondrial and nuclear genome markers are increasingly being used to better understand phylogenetic relationships in order to elucidate whether evolution has occurred among populations of a species (KANDUMA et al., 2012).
The DNA barcode has been proposed as a universal tool for identifying biological diversity, and can be used as a molecular marker (HEBERT et al., 2003). DNA barcoding is based on information gathered from a fragment of approximately 688 base pairs of mitochondrial DNA base sequences (mtDNA) from the cytochrome oxidase I gene (COX-I) of different species. A number of studies have shown that the DNA barcode is a universal code highly effective for the identification of species (HEBERT et al., 2004;BARRETT & HEBERT, 2005).
In addition to the DNA barcode, ribosomal DNA (rDNA) may be used as a nuclear molecular marker. The internal transcribed spacer 2 (ITS-2) is located between the 5.8S and 28S ribosomal subunits of rDNA (CRUICKSHANK, 2002). These genes are arranged in repeated units known as ribosomal cistrons, which have repeated copies and have been used to study phylogenetic relationships (CAMPBELL et al., 1993;SONG et al., 2011;BURGER et al., 2014).
The objective of this study was to infer the phylogenetic and phylogeographic relationships among R. (B.) microplus tick populations in Brazil based on COX-I mitochondrial DNA gene sequences and ITS-2 nuclear DNA, and compare them with R. (B.) microplus tick populations from other countries.

Tick collection
Ticks were obtained from 22 locations in different geographical regions of Brazil ( Figure 1, Table 1). A pool of larvae and 22 engorged females were collected from natural environments. The morphology and identification key was based on Barros-Battesti et al. (2006). Samples were stored in an ultra-freezer at -80 °C for subsequent analysis.

DNA extraction and PCR amplification
DNA extraction was performed on a pool of larvae and individual samples of engorged females using DNAzol  Reagent (Invitrogen) according to the manufacturer's recommendations. The DNA was quantified using a Nano Drop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA), and an A260nm/A280nm ratio above 1.6 was established as acceptable. The sample concentration was corrected to 50 ng µL -1 . Next, the DNA was used in a polymerase chain reaction (PCR) for COX-I and ITS-2 region amplification with new primers designed by our group. PCR generated products of 643bp and 580bp for COX-1 and ITS-2, respectively.
COX-I gene fragments were amplified using the COX-I.2F (5'-CTTCAGCCATTTTACCGCGA-3') and COX-I.2R (5'-CTCCGCCTGAAGGGTCAAA-3') starter oligonucleotides. ITS-2 fragments were amplified with the ITS-2F (5'-CGGATCACATATCAAGAGAG-3') and ITS-2R (5'-CCCAACTGGAGTGGCCCAGTTT-3') primers. The PCR was standardized for a final volume of 25 µL with 1X buffer, 1.5 mM MgCl, 25 mM dNTPs, 10 µM of each primer, 0.5 unit of Platinum  Taq DNA Polymerase High Fidelity (Invitrogen by Life Technologies TM , Massachusetts, USA), and 1µL of DNA at 50 ngµL -1 . PCR conditions were optimized for each reaction, and the annealing temperature was adjusted to suit the primers used. General PCR conditions were 94 °C for 2 min, followed by 35 cycles of 94 °C for 1 min, 56 °C/60 °C (COX-I/ITS-2, respectively) for 30s, 68 °C for 1 min, and a final extension of 68 °C for 3 min. PCR products were visualized in a 2% agarose gel stained with ethidium bromide. The PCR-amplified products were purified using a Purelink  kit (Invitrogen by Life Technologies TM , Massachusetts, USA). Fragments were then cloned with pGEM  -T Easy Vector Systems cloning vector according to the manufacturer's manual (Promega, Madison, WI-USA) and sequenced as described below.

Genetic analysis
All analyses generated in this study were performed in silico. Sanger sequencing was performed in a 48 well plate with robotic instrumentation using T7 Transcription Start (5'-CTAATACGACTCACTATAGGG-3') universal sequencing primer (Promega, Madison, WI-USA). Samples were sequenced using an Applied Biosystems TM ABI 3730 DNA Analyzer and the conditions described by Song et al. (2011). The sequencing reactions were done using the Big Dye  Terminator v3.1 cycle sequencing kit. The runs were performed in 36-cm capillaries using the POP7 polymer, and the sequences were generated by the Sequencing Analysis Software v5.3.1 through the Caller KB Phred program .
Plasmid sequences were identified and removed using an NCBI (VecScreen) tool (ALTSCHUL et al., 1997), and a hash search algorithm was used to remove the contaminant sequences.
Consensus sequences were aligned with the GenBank data using the BLASTN program (ALTSCHUL et al., 1997).
A median-joining analysis implemented in the program Network Version 5.0 (BANDELT et al., 1999) was used for the intraspecific analysis of the evolutionary relationships among haplotypes. Uncorrected (p) pairwise genetic distances were calculated using ARLEQUIN program, version 3.5.1.2 The ARLEQUIN software revealed population structure by means of AMOVA, which uncovers the existence of population differentiation at both intra-and inter-population levels (EXCOFFIER & LISCHER, 2010).

Nucleotide analyses and Haplotype
The phylogenetic and phylogeographic relations of the Brazilian R. (B.) microplus were inferred from the COX-I sequences (22 samples), R. (B.) microplus Brazil (GenBank: KC503261) and the ITS2 sequences (20 samples) obtained from all five Brazilian geographic regions (Figure 1). This is the only study that has compared a representative sample of R. (B.) microplus, which together represent the whole country and the authors were careful to represent each region by just having differentiated biogeography. Assays performed by Dantas-Torres (2015), can demonstrate that certain factors such as climate change and biodiversity boost the expansion of tick populations.
The COX-I dataset showed haplotype diversity (0.767) greater than that revealed by the ITS-2 (0.567) gene. However, the ITS-2 showed greater nucleotide diversity (0.00277) than the COX-I (0.00432) dataset. Haplotype network analysis (COX-I) revealed ten distinct haplotype clusters among Brazilian populations, with no clear separation by tick populations/geographical areas, indicating that an overlap of the ten genetically divergent haplotypes of R. (B.) microplus (Figure 2A) exists.

Genetic distance, genetic differentiation and gene flow
In the uncorrected "p" distance matrix, the COX-I gene indicated stronger resolving power than ITS-2 (1.36-11.00) ( Table 2). The genetic distance of the COX-I H-2 haplotype of R. (B.) microplus indicated 0.73% of divergence within this population (red value, Table 2). Overall, a high level of significant genetic differentiation was observed among all populations (P < 0.001).
Great genetic differentiation (FST = 1.00), revealed by the COX-I gene, was found with in all haplotypes (Table 3). Additionally, genetic differentiation between H-2 haplotypes and others varied between 0.47-0.89 (Table 3). High levels (30.00) of gene flow (Nm) were observed among tick populations as shown by the COX-I (data no shown).
Among the ten haplotypes studied, R. (B.) microplus haplotype H-3 from Rondônia (RO) and Paraíba (PB) were genetically diverse as supported by the highest haplotype and nucleotide diversity based on COX-I (Table 3). Higher inter-population levels (82.69%) and (FST = 0.8269) revealed the existence of population structure in all haplotypes analyzed by means of AMOVA (Table 4).
Livestock transition in Brazil is very intense. Whether these regions actually have different populations not found in other regions can only be answered by analyzing a larger number of samples from various municipalities in the States of Rondônia and Paraíba.
All haplotypes are derived from haplotype H-2. However, it is not possible to assert where haplotype H-2 has originated. Additionally, there is no haplotype H-2 in northeastern Brazil.
The fact that the country has cattle production for imported and local breeds and that animals are transported between different regions do not allow us to infer the exact origin of haplotypes. The natural dispersal ability of ticks might occur in parallel with human-mediated dispersal, causing genetic admixture and high gene flow in Brazilian populations.

Phylogenetic analyses
Bayesian analyses produced phylogenetic trees with high valued probabilities (Figure 2 and Figure 3). The 10 haplotypes (H1-H10) generated from the COX-I sequences of Brazilian R. (B.) microplus were subjected to phylogenetic analyses. By the topology of the phylogenetic trees generated, the R. (B.) microplus Brazilian population comprises two clades: clade A (H3) and clade B (H1, H2 and H4 to H10) (Figure 2A). Different than results revealed by the haplotype network analysis ( Figure 2B), according to Bayesian phylogenetic tree, haplotypes H5, H6 and H8 were associated with ticks from the North region of the country and are derived from haplotype H-4. The tree also indicates that all haplotypes derive from H-2.  (Figure 3).
The intraspecific genetic distance or haplotype frequency, like that of haplotype H2 reported here, are notably higher than previously described for R. (B.) microplus sensu Low et al. (2015). Furthermore, the two, and perhaps three, genetic assemblages inferred from Brazilian R. (B.) microplus are more genetically diverse than those reported from other regions (Figure 3).
This study sought to obtain high support for the topology of tree and the COX-I molecular marker was found to have the advantage of presenting greater mutation rates than the nuclear marker, corroborating the data presented by Lv et al. (2014). However, the ITS-2 molecular marker failed to resolve discrepancies between the populations of R. (B.) microplus and this study corroborates the findings of Burger et al. (2014) in which they showed that the ITS-2 molecular marker has insufficient power to distinguish between very close species.
The comparative analysis in this study shows that R. (B.) microplus ticks represent at least two different populations. Moreover, this study also attempted to investigate whether haplotypes are ) microplus clade B is more closely related to R. annulatus, with which it forms a sister group relationship, corroborating the data by Burger et al. (2014) (Figure 3). distributed according to the phylogeographical areas. It is unknown whether this indicates a cause-effect association of haplotypes and geographical areas.
However, the prevalence of ticks has been typically associated with cattle breeds, the highly invasive and widespread movement of R. (B.) microplus may facilitate its occurrence (IBELLI et al., 2012;BUSCH et al., 2014).

Conclusion
This study provides new insights into the distinct genetic assemblages of the tick R. (B.) microplus in Brazil and reveals that the Brazilian R. (B.) microplus population consists of at least two different populations. In addition, future studies using a large number of samples from the Americas using COX-1 may be able to show how these populations differ in their inter-and intra-relationships.