Analysis of Bm86 conserved epitopes: is a global vaccine against Cattle Tick Rhipicephalus microplus possible?

The cattle tick Rhipicephalus microplus causes significant economic losses in agribusiness. Control of this tick is achieved mainly through the application of chemical acaricides, often resulting in contamination of animal food products and of the environment. Another major concern associated with acaricide use is the increasing reports of resistance of this tick vector against the active ingredients of many commercial products. An alternative control method is vaccination. However, the commercially available vaccine based on a protein homologous to Bm86 exhibits variations in efficacy relative to the different geographical locations. This study aimed to identify antigenic determinants of the sequences of proteins homologous to Bm86. Phylogenetic analyses were performed to determine the extent of divergence between different populations of R. microplus to identify the sequence that could be used as a universal vaccine against the multiple geographically distinct populations of R. microplus and related tick species. Considering the extensive sequence and functional polymorphism observed among strains of R. microplus from different geographical regions, we can conclude that it may be possible to achieve effective vaccination against these cattle ticks using a single universal Bm86-based antigen.


Introduction
The morpho taxonomy of Rhipicephalus microplus complex has been challenged in recent years, and until now, they are known as a complex composed of five taxa, namely R. australis, R. annulatus, R. microplus clade A sensu, R. microplus clade B sensu and R. microplus clade C sensu (BURGER et al., 2014;LOW et al., 2015;CSORDAS et al., 2016). This complex causes significant economic losses to livestock. In Brazil alone, the potential damage caused by this species amounts to US $3.24 billion per year (GRISI et al., 2014).
Control of these parasites is mainly achieved through the use of acaricides. However, the pressure of use has genetically selected tick populations resistant to these products, resulting in serious problems such as a shortage of new acaricidal chemical entities and the accumulation of residues in both the environment and products of animal origin (ANDREOTTI et al., 2011;RECK et al., 2014;SINGH et al. 2015). Therefore, the search for alternative control methods has become increasingly necessary.
Tick control through vaccines is an alternative to the use of acaricides. Some examples of vaccines that use the Bm86 protein as an immunogen against ticks include TickGard PLUS (formerly, Intervet Australia)  which is no longer commercially available and Gavac TM (Heber Biotec) (CANALES et al., 1997) which is the only one still available and, more recently, Bovimune Ixovac, Mexico (BOVIMUNE IXOVAC, 2018).
To understand why many vaccine strains should be tested, a better understanding of the history of the cattle tick R. microplus is necessary. Based on the history of the dissemination of R. microplus according to Barré & Uilenberg (2010), R. microplus originated in the southern and southeastern regions of Asia and was thus regarded as one of the most successful invasive tick species, colonizing wide areas in Central and South America, South East Asia, Australia, and islands in the Pacific Ocean.
Vaccines, such as TickGARD PLUS (formerly, Intervet Australia), against the cattle tick in Australia, which was previously named R. microplus until it was morphologically and genetically distinguished from Rhipicephalus australis and now is part of the R. microplus complex (LABRUNA et al., 2009;ESTRADA-PEÑA, et al. 2012;LOW et al., 2015;CSORDAS et al., 2016) are effective against susceptible strains of Australian origin. Numerous tests have been performed and might exhibit a different pattern of sensitivity and therefore induce poorer control when applied against R. microplus in America (SUTHERST & BOURNE, 2009).
Bm86-derived vaccines cause reduced weight of engorged adult female ticks, reduced egg mass weight, reduced numbers of ticks in the field over one generation, a decrease in the reproductive capacity of R. microplus females, decreased frequency of treatments with acaricides, and a positive impact on the implementation of integrated control programs (JONSSON et al., 2000;ODONGO et al., 2007;VARGAS et al., 2010). The TickGARD PLUS (formerly, Intervet Australia) and Gavac TM (Heber Biotec) vaccines were shown to be 46.9% and 49.2% effective, respectively, against a Brazilian R. microplus strain (ANDREOTTI, 2006). This effectiveness is lower than that reported globally for other strains (RAND et al., 1989;RICHARDSON et al., 1993;PATARROYO et al., 2002). Bm86 from the Campo Grande strain (GenBank: ACA57829) collected from cattle under field conditions in Campo Grande municipality, Brazil, exhibits differences in hydrophobic regions compared with the vaccine proteins. These differences may contribute to the binding of antibodies to the target protein, which may explain the variation in efficacy in other regions worldwide (ANDREOTTI et al., 2008).
Humoral and cellular immune responses are important to an organism to react against specific regions of pathogens via a process known as the adaptive response. Vaccines based on subunits or synthetic peptides, which is a strategy that focuses on epitopes recognized by B and T cells and the major histocompatibility complex (MHC), that could elicit a specific immune response have been studied in recent decades (BEN-YEDIDIA & ARNON, 1997). In silico epitope prediction tools are necessary to evaluate which regions of such proteins of interest could be candidates for a peptide-based vaccine. Initially, the predictions focused on epitopes that bind on MHC class I molecules for the development of cellular immune response against virus-infected cells and cancer (KAST et al., 1991;TOES et al., 1996). Examples of some classes of methods used in predictions are quadratic programming, linear programming and the use of sequence profiles obtained by clustering known epitopes of a given MHC allele to score candidate peptides (FLOREA et al., 2003).
The interaction of innate and acquired immune responses to ticks is what will determine the success of the host or the parasite. The innate response consists of an inflammatory reaction and a hemostatic process, such as vasoconstriction, clot activation and platelet aggregation (SZABÓ, 2008). Hypersensitivity related to histamine release by mast cells is noted in resistant animals (RIEK, 1962). Resistant animals exhibit increased cell infiltrate in the tick attachment site, which suggests that cellular immune response is important to confer natural resistance to the host (SZABÓ & BECHARA, 1999).
Studies on peptide-based vaccines against ticks have been performed in the last decade. In silico analysis of Bm86 gut glycoprotein revealed that it contains antigenic epitopes, one of which has elicited greater than 80% efficacy in immunized cattle, against R. microplus (PATARROYO et al., 2002). Monoclonal antibodies against a peptide designed based on B cell epitopes from Bd86 (a Bm86 ortholog from Rhipicephalus decoloratus) determinar a extensão da divergência entre diferentes populações de R. microplus com o objetivo de identificar a sequência que poderia ser usada como vacina universal contra as múltiplas populações geograficamente distintas de R. microplus e espécies de carrapatos relacionados. Considerando-se a extensa sequência e o polimorfismo observados entre linhagens de R. microplus de diferentes regiões geográficas, podemos concluir que pode ser possível obter uma vacinação efetiva contra esses carrapatos bovinos utilizando um único antígeno universal baseado em Bm86.
could bind to gut cells from several Rhipicephalus species (including R. microplus) (KOPP et al., 2009). Reverse vaccinology approaches to identify new tick antigens are currently under investigation. However, the number of tick genome databases available is limited given that these parasites are eukaryotes and have a large genome (LEW-TABOR & RODRIGUEZ-VALLE, 2016). The post-genomic era revolution has produced massive sequence data on genomes, transcriptomes and proteomes, including those of R. microplus (BARRERO et al., 2017). This study aimed to identify antigenic determinants of the sequences of proteins homologous to Bm86 that could be used as a universal vaccine against the multiple geographically distinct populations of R. microplus. Variations in the Bm86 locus sequence have been reported in the literature (RODRÍGUEZ et al., 1994;COBON et al., 1995;GARCÍA-GARCÍA et al., 1999). This study aimed to determine the effect of sequence variability within the Bm86 locus on the major epitopes of the Bm86 protein as a means of identifying possible isolates that confer protection against R. microplus strains of various geographic origins.

In silico analyses for epitope predictions
For the analyses, we used tools based on the major histocompatibility complexes (MHCs) of cattle (bovine leukocyte antigen (BoLA)) and mice (H-2 antigen). Mouse MHC class I and II predictions are exclusively based on mouse H-2 antigen because it is an animal model commonly used to screen candidate antigens before testing in cattle (AGUIRRE et al., 2016;CONTRERAS et al., 2016;LEW-TABOR & RODRIGUEZ-VALLE, 2016). A variety of algorithms were run to generate a combination of results capable of predicting qualities that altered the statistical probability of a peptide sequence having antigenic potential for the host immune system in question (bovine MHC I via BoLA loci and mouse MHC I and II via the H-2 locus). We used the following prediction parameters to predict three major epitopes based on the Brazilian Bm86-CG (ACA57829) amino acid sequences published in GenBank: cattle MHC class I binding epitopes for the BoLA-AW10, BoLA-N:00101 and BoLA-D18.4 alleles (IEDB, 2018a); binding to the mouse MHC I H-2-Qa1, H-2-Dd and H-2-Ld alleles (IEDB, 2018a); binding to the mouse MHC class II H-2-IAb and H-2-IAd alleles (IEDB, 2018b); linear epitopes for B lymphocytes (IEDB, 2018c); epitopes exposed on the protein surface in the tertiary structure (IEDB, 2018d); prediction of transmembrane helices (DTU bioinformatics, 2018a); prediction of intrinsically disordered protein regions (IUPRED2A, 2018); prediction of the signal peptide (DTU bioinformatics, 2018b); and prediction of glycophosphatidylinositol (GPI) anchors (PIERLEONI et al., 2008).
The function of each in silico tool mentioned above was to predict Bm86 regions with the potential to trigger cellular immune responses in cattle (the only available source for predictive tools for several animal species, netMHCpan, is based on MHC class I); predict regions with the potential to trigger cellular and humoral immune responses in mice (MHC class I was evaluated for comparison with cattle results) (VORDERMEIER et al., 2003;NENE et al., 2012;RODRIGUEZ-VALLE et al., 2013;AGUIRRE et al., 2016); predict B lymphocyte epitopes to ensure immunological memory and generate rapid antibody responses in immunized hosts exposed to the tick to avoid selecting epitopes in transmembrane helices, as these structures can mask epitopes because they are anchored to the cell membrane; evaluate intrinsically disordered regions, which are of some interest because they do not have a constant tertiary structure and the epitope is more likely to be exposed to the immune system when it is located in a flexible region of the protein, to avoid epitopes located in signal peptides, which are typically cleaved after the protein reaches its final destination; and predict glycophosphatidylinositol anchors (GPI) that may increase triggering of immune responses in the host.
To determine whether the epitopes predicted with the algorithms listed above are conserved among the different Bm86 sequences deposited in GenBank, all sequences were analyzed using the Mega 6.0 program (TAMURA et al., 2013), and an alignment was constructed using Gonnet modeling.
Protein BLAST (NCBI, 2018) was used to perform alignments of the epitopes predicted with Bos taurus protein sequences available in the SwissProt database using the BLASTP 2.5.1 program. This analysis aimed to ensure that the predicted sequences did not have homologies with host proteins, thereby minimizing autoimmune responses following vaccination.

Sequence alignment and phylogenetic tree construction
The Bm86-CG protein sequence (GenBank: ACA57829) and the three predicted epitope sequences were aligned with the sequences available from GenBank using the BLASTp program. In this way, a database was constructed that contained all similar sequences obtained from the analysis. The Mega 6.0 program (TAMURA et al., 2013) was applied to align the sequences taken from GenBank using the Gonnet protein weight matrix.
A Bayesian phylogenetic analysis was performed using the MrBayes 3.2.6 program (RONQUIST & HUELSENBECK, 2003). An amino acid analysis was performed using the Dayhoff model (DAYHOFF et al., 1978). For the data set used in this study, approximately 500,000 generations were found to be sufficient for topology and were plotted using the FigTree 1.4.2 program (TREE BIO, 2016). All analyses for Bm86 epitopes and partial protein of Bm86 were initiated with random starting trees and were run for 1 X 10 6 generations, with trees being sampled every 1,000 generations. To determine the stationarity of the Markov chain, the log-likelihood scores of sample points were plotted against generation time. The first 25% of samples was estimated as burn-in for each data set. The remaining samples were retained for generating consensus trees. Each sample included a tree topology that incorporates branch length and substitution model parameter values. These topologies were used to generate a 50% majority rule consensus tree, with the percentage of sample recovering any particular clade representing the posterior probability of a clade  (ODONGO et al., 2007, CANALES et al., 2008, were collapsed (tree's triangle). Species of the tick Hyalomma detritum (AEK31101), H. detritum (AEK31102), H. anatolicum (ACD14076) were used as outgroups in phylogenetic analyses.

Results
The bioinformatics analysis predicted three major relatively conserved regions (Figure 1) in the Bm86 sequences from different geographical regions. The Bm86-CG protein sequence (GenBank: ACA57829) was used as a reference because it originated from a strain of R. microplus from the study region (ANDREOTTI et al., 2008).
Epitope 1 (Figure 2) formed a distinct clade with high convergence for R. microplus populations from the geographical regions of Brazil, the United States and Mozambique, Israel and some populations from Thailand. A collapsed group was formed that included the variants from Thailand. Immediately below was formed with branch with variants from Thailand, Figure 1. Analysis of Bm86-CG integration using free algorithms available online, based on T and B cell epitope predictions. Each tool predicts a certain characteristic based on calculations that use the properties of each amino acid in the primary protein sequence to generate graphs indicating which regions of the primary structure have a statistically significant probability of having the desired characteristic. For this purpose, a threshold (horizontal line parallel to the x-axis) is plotted. All regions of the graph that exceed this threshold are considered in the present analysis. By superimposing the graphs from different algorithms, we can better visualize which regions share the largest combination of antigenic characteristics. USA, Mexico, Egypt and the species Rhipicephalus annulatus (GenBank: ACR19242 and ACL27210). The other clades were separated by species divergent for epitope 1.
For epitope 2 (Figure 3), the phylogenetic tree presented a clade with high convergence for R. microplus populations from the geographic regions of USA (FREEMAN et al., 2010), Mexico (CANALES et al., 2009), Israel (CANALES et al., 2009, and variants of Thailand (KAEWMONGKOL et al., 2015). High convergence with the R. annulatus species and R. decoloratus were observed in this same clade. The other clades were separated by species divergent for epitope 2.
For epitope 3 (Figure 4), a large clade with high convergence was formed for all of the Brazilian R. microplus populations reported by Sossai et al. (2005) together with the Brazilian strain reported by Andreotti et al. (2008). In the same clade, high convergence was observed for the geographic regions Other clade was formed for the variants from Thailand (KAEWMONGKOL et al., 2015), the United States (FREEMAN et al., 2010) and the orthologous R. annulatus species (GenBank: ADR19242, ADQ19693, ADQ19691 and ACL27210). The other clades were separated by species that showed divergence of the Bm86 protein.

Discussion
According to De Groot (2006), a vaccine must promote and trigger a strong B cell and/or T cell response to be effective. Because complex cellular interactions occur during the development of an immune response, mapping and characterization of epitopes is of fundamental importance for the selection of peptides with potent antigenicity as candidate vaccine components. The MHC genes of cattle were mapped from bovine autosome 23 (BTA 23) (FRIES et al., 1986;FRIES et al., 1993) and are referred to together as bovine leukocyte antigen (BoLA). This organization distinguishes BoLA from the MHC in humans and rodents, in which the MHC genes are located on chromosome 6 (HLA) and 17 (H-2), respectively (LEWIN, 1996). McShane et al. (2001) used FISH for mapping genes and revealed that tumor necrosis factor a (TNF-α), heat shock protein 70 (HSP70.1) and 21 steroid dehydrogenase (CYP21) are closely linked in the region of BTA23 band 22 along with BoLA class I genes. Therefore, this gene is important for immunological functions.
Few variations that might be involved in the process associated with the global variability in vaccine efficacy were observed in the three Bm86-CG epitopes that were predicted here to be the most antigenic. In fact that, these epitopes contained regions with major similarities to homologous domains of Bm86 from strains from different geographical regions (RAND et al., 1989;RICHARDSON et al., 1993;DE LA FUENTE et al., 1999;PATARROYO et al., 2002;ANDREOTTI, 2006;ANDREOTTI et al., 2008). Nevertheless, the variability observed in the Bm86 molecule may explain the differences in the efficacy observed when cattle vaccinated with rBm86-CG and rBm86 were challenged by infesting with the same strain, R. microplus Campo Grande (CUNHA et al., 2012;ANDREOTTI, 2006). Peconick et al. (2008) studied polymorphisms and observed that vaccine efficacy was not decreased despite changes in the amino acid sequences of polypeptide SBm7462 in two Brazilian isolates compared with the sequence of the polypeptide obtained from other isolates by Patarroyo et al. (2002). Despite the variability of proteins homologous to Bm86 and their impact on the immune response, it must be borne in mind that the genes of two MHC class II alleles of cattle have been described with a deletion of three base pairs in their nucleotide sequences, leading to the deletion of a lysine at the antigen recognition site. This deletion was correlated with increased immune response against the commercial vaccine TickGARD PLUS (formerly, Intervet Australia) (SITTE et al., 2002). Thus, the influence of the genetic variation of the MCH class II alleles on the cattle immune response against vaccine antigens is also evident.
The criterion established based on the prediction of B lymphocytebinding linear epitopes had the objective of selecting sequences against which the host could trigger a humoral response in the shortest time following re-exposure to the antigen and that could be used to stimulate immune memory with boosters using vaccine antigen formulations. Were predicted Bm86 regions exposed on the surface of the tertiary structure as a basis for selecting peptide sequences accessible for binding to antibodies produced by the host against each peptide.
Predictions of transmembrane helices were used to define the regions that should be avoided. These regions are usually located next to the lipid bilayer of the cell membrane, hindering their access by host antibodies. Bm86-CG lacks transmembrane helices according to the results of our TMHMM analysis (data not shown). The prediction of intrinsically disordered regions was performed to select peptide sequences that are interesting targets because these peptides do not have a constant tertiary structure. This phenomenon may facilitate exposure of the epitopes in the native protein, thereby facilitating access to the antibodies.
Epitopes 1 and 2 located near the C-terminus could have this characteristic to our results for the prediction of intrinsically disordered protein regions.
The signal peptide is usually cleaved during post-translational modification of the protein. Signal peptide prediction is used in reverse vaccinology to select proteins that are likely to be secreted by the target cell of the agent. However, these sequences should not be included in the selected peptide sequences because they typically do not form part of the final structure of the mature protein. The Bm86-CG protein does not appear to have a signal peptide according to the signalP algorithm (data not shown). This fact is due to the Bm86-CG sequence deposited in GenBank is a partial sequence of the whole protein, which lacks amino acids at the signal peptide site. On the other hand, full Bm86 sequences, such as from isolate N1 from a Thailand strain, present a signal peptide of 20 amino acids, as showed by signal algorithm (data not shown). The prediction of GPI anchors in combination with signal peptide prediction allows us to infer that the protein is anchored in a biological membrane. However, we did not observe a GPI anchor based on PderGPI results (data not shown). Additionally, we can infer its possible location, which is a factor that can enhance the triggering of an immune response in the host. Alignment of the sequences of Bm86 proteins from distinct isolates allowed us to identify highly conserved regions (Supplementary Material Figure S1) that could be used to create a vaccine that is based on subunits/peptides. However, these regions do not coincide with the antigenic regions of bovine Bm86 that are capable of binding to the bovine MHC. Thus, when the animal is exposed to Bm86, its antigen-presenting cells (APCs) will process the antigen but will present only peptides from variable regions to CD4+ T lymphocytes via MHC class II. This process may result in nonspecific immune responses against these different strains.
A review by Parizi et al. (2012) indicate that the decades-long search for a universal tick vaccine is making progress, with such a vaccine likely to be based on multiple antigens and consequently successful vaccinations against multiple tick species. Our study shows progress towards a universal tick vaccine, that presented three conserved epitopes for regions of the Bm86 protein. Investigation of Hue et al. (2017) an experimental vaccine was developed based on the orthologous R. australis Bm86 sequence identified from local R. australis strains and a recombinant protein expressed in Escherichia coli. The use of the vaccine reduced the population by 74.2% to each generation. These data once again demonstrate that it is possible to direct for a universal anti-tick vaccine.
De la  in the review work on vaccine strategies reports the caution should thus be placed on extrapolating the protective abilities of a specific antigen between tick species, as was seen in the case of Bm86 vaccines. It is worth mentioning, that new approaches to tick control depend on host-parasite interactions. Thus, the efficacy of a vaccine also stems from the clarification of the different immunological profiles presented by the hosts, where bovine tick resistance is a hereditary phenotype (KASHINO et al., 2005;CARVALHO et al., 2010).
The present study showed that based on Bm86, R. microplus populations were homogeneous and displayed high convergence among the isolates studied by García-García et al. (2000), Sossai et al. (2005)  microplus -THAILAND (KAEWMONGKOL et al., 2015); and Bm86 orthologs of R. microplus, such as R. annulatus (ABY58969) and R. decoloratus (ABY58970, ABG21130, ABG21131). Therefore, the results in the present study (Figures 4 and 5) suggest that the polymorphisms were not computed as significant changes and that all these peptides (Figures 2, 3 and 4) with polymorphic sites are more conserved when compared with the partial Bm86 protein sequence ( Figure 5).
Interestingly, the strain reported by Willadsen et al. (1989) (P20736) and from Hüe et al. (2017) (ATW75471 -76) remained in the clade along with the other Bm86 strains (Figure 4). It is worth remembering that the cattle tick, previously named R. microplus until its morphological and genetic distinction to R. australis (ESTRADA-PEÑA et al., 2012). These strains obtained high convergence with the other strains of R. microplus for epitope 3. As this strain could produce an ortholog of R. microplus Bm86, this peptide could become a vaccine candidate for both R. microplus and R. australis. Canales et al. (2009) showed that the Ba86 protein (an ortholog of Bm86 from R. microplus) is effective in controlling infestation with R. annulatus and R. microplus. This ortholog was shown to have low convergence in the tree with the partial Bm86 protein (ACR19242) ( Figure 5). However, the chosen epitopes have been preserved for this ortholog, and may continue to be efficient in cross-immunity. This shows even orthologs have conserved regions and could be effective antigenic and immunogenic regions.
Differences were observed in the phylogenetic trees of epitope 1 and for total protein Bm86 (Figures 2 and 5) in the Thailand variants, which formed a distinct clade from the other R. microplus sequences. These differences could hinder the development of a global vaccine for these variants. However, epitopes 2 and 3 (Figures 3 and 4) presented a larger and more homogeneous clade with convergence among the R. microplus Thailand populations (Figure 3), suggesting that this epitope was more conserved and could generate a potential immunogen among the Thailand variants described by Kaewmongkol et al. (2015). Kamau et al. (2016) demonstrated that intestinal cDNA from four Rhipicephalus appendiculatus field populations revealed genotypic polymorphisms and suggested that other factors such as exposure to innate immune components during blood feeding may direct the selection pressure that leads to the observed polymorphisms. Additional polymorphisms may have occurred during the evolution of the R. microplus populations (SOSSAI et al., 2005). Population studies performed by Csordas et al. (2016) with the cytochrome c oxidase subunit 1 (COI-1) gene indicated the presence of two haplotypes that differ from other Brazilian R. microplus populations. Thus, their antigenic epitopes may have become distinct from the epitopes of other populations. This possibility may indicate that most antigenic characteristics are not located in regions conserved in the Bm86-CG protein but instead are located in more variable regions.
Importantly, some convergence was observed for each phylogenetic tree of each selected epitope, resulting in their distribution among several geographic areas. Conversely, when the entire amino acid sequence of the Bm86-CG protein was analyzed (Figure 5), the polymorphisms formed several clusters associated with various geographic regions. This work presents peptides for a vaccine that has point polymorphisms between conserved regions of the epitopes, but this does not interfere in the immunogenic and antigenic regions and may be candidates for a global vaccine. The in silico analysis presented here can be an important tool for related tick populations that have not been challenged with the Bm86 antigen in order to minimize future economic problems arising from infestations.

Conclusion
The phylogenetic analysis of the isolated epitopes and the Bm86-CG protein revealed geographical patterns among R. microplus tick isolates collected from the database; these isolates were classified into a major clade based on their amino acid sequences. These sequences may contain epitopes with amino acid residues that have both conserved and antigenic regions. However, this characteristic is found in most variable regions of the protein.
Considering the extensive sequence and functional polymorphism observed among strains of R. microplus from different geographical regions, we can conclude that it may be possible to achieve effective vaccination against these cattle ticks using a single universal Bm86-based antigen. These point polymorphisms in some amino acids do not form in a less conserved peptide, not altering their antigenic or immunogenic capacity. As genomic technologies in vaccine development continue to improve and evolve, sequencing of various tick genomes can identify tick vaccine candidates globally applicable to various species of ticks.