In Silico Survey and Characterization of Babesia microti Functional and Non-Functional Proteases

Human babesiosis caused by the intraerythrocytic apicomplexan Babesia microti is an expanding tick-borne zoonotic disease that may cause severe symptoms and death in elderly or immunocompromised individuals. In light of an increasing resistance of B. microti to drugs, there is a lack of therapeutic alternatives. Species-specific proteases are essential for parasite survival and possible chemotherapeutic targets. However, the repertoire of proteases in B. microti remains poorly investigated. Herein, we employed several combined bioinformatics tools and strategies to organize and identify genes encoding for the full repertoire of proteases in the B. microti genome. We identified 64 active proteases and 25 nonactive protease homologs. These proteases can be classified into cysteine (n = 28), serine (n = 21), threonine (n = 14), asparagine (n = 7), and metallopeptidases (n = 19), which, in turn, are assigned to a total of 38 peptidase families. Comparative studies between the repertoire of B. bovis and B. microti proteases revealed differences among sensu stricto and sensu lato Babesia parasites that reflect their distinct evolutionary history. Overall, this data may help direct future research towards our understanding of the biology and pathogenicity of Babesia parasites and to explore proteases as targets for developing novel therapeutic interventions.


Introduction
Human babesiosis caused by Babesia microti is a malaria-like tick-borne zoonotic disease, first described in the 1950s in the USA, with an increasing number of cases reported ever since in this and other countries around the world [1]. Infections proceed asymptomatic or are accompanied by mild or moderate signs in immunocompetent patients but often lead to severe disease and even death in neonates and the elderly or immunocompromised adults [2]. Some wild rodents act as natural reservoirs of B. microti, where the parasite is transmitted both by bites of Ixodes sp. ticks, as well as transplacentally [3,4]. Humans are dead-end hosts and suffer accidental infections mainly through tick bites. Transplacental and blood transfusion-related transmissions have also been documented [1,5,6].
Currently, there is no specific therapy for B. microti human babesiosis [7]. The recommended therapeutic drugs to treat B. microti infections are azithromycin plus atovaquone as the first choice or a combination of clindamycin and quinine as an alternative [7,8]. However, the reported appearance of B. microti parasites resistant to the first two drugs in chronically infected patients and the negative side effects of the latter two call for the development of alternative therapeutic strategies and increased investments in this field [2,[7][8][9][10]. Inactive T Protease types are color-coded: blue: aspartic, light yellow: cysteine, orange: threonine, gray: metallo, and yellow: serine proteases, with darker colors for peptidases predicted as active. The relevant amino acid positions in the active site needed for catalytic activity are included for each active protease. Inactive proteases have a predicted protease domain but lack one or more of the functional amino acids. Paralog groups within each family are underlined. When more than a paralog group is present in a family, different underlining styles are used for each group. (a) ROM4, (b) ROM7 (active protease, wrongly predicted as inactive in MEROPS), and (c) ROM8 [29]. (*) Paralogs of the M41 and the S59 families lacking peptidase domains are shown separately (Supplementary Tables S1 and S2). T: Transcribed genes in the intraerythrocytic stage [30]. H or M: High or medium levels of protein expression detected in Reference [31].
Twenty of the proteases presented in this study are not included in the MEROPS database and were identified by homology searches. In addition, five of the proteases listed in Table 1 presented duplicated MEROPS entries, likely because of the use of different sources of peptide sequences in this database (Table 1). On the other hand, a number of B. microti proteases have been annotated in GenBank as hypothetical proteins, uncharacterized proteins, or following the designation of another conserved domain also present in the sequences, and they are, thus, not identifiable by searches using keywords such as protease or peptidase. Importantly, despite the exhaustive search carried out to produce the list presented in Table 1, the presence of additional protease-coding genes in the B. microti genome that passed inadvertently in this study cannot be ruled out. In addition, it should be noted that the predicted catalytic activity for some threonine proteases (Table 1) could not be determined beyond doubt and needs to be confirmed experimentally.
All of the listed active and non active proteases are transcribed in B. microti merozoites, suggesting they likely fulfill a relevant functional role in this parasite stage [30]. In addition, 17 active and non active proteases were identified in the proteomic profile of B. microti during the acute infection of a hamster model (Table 1) [31]. Proteases that remained undetected might be expressed in low amounts at the intraerythrocytic stage or bear physicochemical characteristics that preclude detection by the experimental approach employed in this study [31].
Localization predictor algorithms located most of the identified proteases in the cytoplasm. Four proteases were predicted as extracellular and six as lysosomal. The latter might reach the extracellular milieu by the fusion of vacuoles with the plasma membrane, as has been shown for Tetrahymena thermophila, a free-living protozoon belonging together with Apicomplexa to Alveolata. However, this mechanism has not been demonstrated for Babesia spp. [32]. Other predicted locations include the nucleus, the mitochondria, a plastid (which would correspond to the apicoplast, in this case), the plasma membrane, the endoplasmic reticulum, the Golgi apparatus, and the peroxisomes (Supplementary Table S3). Importantly, these predictions are only tentative until they have been experimentally confirmed. Additionally, the used algorithm is not able to predict the location of proteases in Apicomplexa-specific secretory organelles, such as rhoptries and micronemes, where the trafficking of proteases has been shown to occur in Plasmodium and Toxoplasma [33].

Aspartic Proteases
Seven aspartic proteases were found in the B. microti-predicted proteome, all of which bear the aspartate and, in the case of the A1 family, also the phenylalanine or tyrosine residues in their active sites, needed to display catalytic activity (Table 1).
Interestingly, a recent transcriptomics study involving four of the five A1 aspartic protease genes of B. microti showed stage-associated expression for two of them. Thus, while BmR1_01G02485 (encoding cathepsin E-B or ASP2) displayed higher expression in mouse blood intraerythrocytic stages than in the stages present in I. ricinus gut or salivary glands, the opposite was true for BmR1_03g03850 (encoding ASP6). These results suggest a role for ASP2 in processes connected to the asexual reproduction of the parasite and/or gametocyte formation and, for ASP6, in zygote and/or kinete development, kinete dissemination in tick tissues, including salivary glands, and sporogony. For the other two studied A1 aspartic protease-encoding genes, BmR1_04g07350 and BmR1_04g05270, corresponding to ASP3 and ASP5, respectively, expression was similar in the three stages, suggesting a role in invasion both of erythrocytes and tick cells or in other cellular processes such as secretion or the trafficking of proteins [34].
Aspartic proteases have been proposed as chemotherapeutic targets against B. microti. Indeed, the aspartic protease inhibitors Lopinavir and Atazanavir, which are well-tolerated drugs used in HIV patients, were shown to be potent suppressors of B. microti infection in vitro, as well as in vivo, in a mouse model [24]. It is unknown which parasite aspartic protease is affected by these inhibitors, but one candidate is the signal peptide peptidase (SPP, XP_021338622), which has a critical role in the maintenance of the homeostasis of the endoplasmic reticulum. Consistent with this view, in the case of P. falciparum, these inhibitors were effective in blocking SPP activity and in vitro parasite growth [24,35]. Notably, B. microti and Plasmodium sp. SPP proteins are orthologous (results not shown) but do not have a counterpart in B. bovis (Supplementary Table S4) or any other Babesia sp.
Hemoglobin is certainly the main protein source available for the nutrition of intraerythrocytic parasites. The sequential steps for hemoglobin degradation by Plasmodium sp., as described by Guzman et al. 1994 [36], start with the unwinding of the molecule and partial digestion by aspartic proteases, followed by cysteine protease cleavage, which yields protein fragments that are finally degraded by exopeptidases, generating free amino acids useful for parasite nutrition. In Plasmodium sp., the first part of this process takes place in the food vacuole and involves the aspartic proteases Plasmepsins I-IV and is then followed by the action of papain-like cysteine proteases in the erythrocyte cytoplasm [37]. There is no evidence that a food vacuole is present in B. microti, and consistent with its absence, Plasmepsins I-IV homologs cannot be found in this parasite. The lack of these enzymes has been used as an argument to postulate that B. microti is not able to degrade hemoglobin [38]. However, it may be hypothesized that other B. microti aspartyl proteases of the A1 family (Table 1), likely secreted to the erythrocyte cytoplasm, are able to initiate hemoglobin degradation, such as pepsin A (XP_021337801), predicted to have a signal peptide and, thus, be exported to the erythrocyte cytoplasm (Supplementary Table S3). To find out whether this is the case or there is an alternative protein source available for the nutrition of the intraerythrocytic trophozoite and merozoite stages would need experimental exploration.

Cysteine Proteases
Cysteine proteinases are involved in the essential biological roles of Apicomplexan parasites [13,39,40]. They are present in B. microti with at least 27 members, of which 18 are predicted to be catalytically active (Table 1).
In P. falciparum, the papain-like falcipain-2 and falcipain-3 peptidases of the C1A family have attracted the most attention among cysteine proteases as potential therapeutic targets against malaria [41]. As mentioned above, these enzymes participate in hemoglobin degradation in the intraerythrocytic stage of the parasite, and, in addition, falcipain-2 has been shown to cleave erythrocyte cytoskeletal proteins during egress from the host cell [42,43]. Falcipain-2 orthologs have been characterized in B. bovis, B. bigemina, and B. ovis and have been named bovipain-2, babesipain, and ovipain-2, respectively. Similar to their P. falciparum counterpart, they are expressed inside merozoites and also released to the erythrocyte cytoplasm, consistent with the dual role described for falcipain-2 [44][45][46][47]. The significant impairment of the in vitro growth of B. ovis and B. bovis merozoites by antibodies against ovipain-2 and a papain-like C1A cysteine protease, respectively, indicate a relevant role of this type of enzymes in the propagation of the asexual stages of Babesia spp. [47,48].
The B. microti ortholog of falcipain-2 (XP_012650559) has four paralogs (three active proteases and one non-protease homolog), one of which is 100% identical (XP_012650562; Table 1). The corresponding genes for these two identical proteins are located on the same strand of chromosome 3, separated by a~5 kb intergenic region, where two unrelated genes are found in the opposite strand. Predictor algorithms localized XP_012650559 and XP_012650562, either within lysosomes or other vacuoles or secreted through a non classical pathway (Supplementary Table S3). This predicted localization agrees well with that described for their counterparts in B. bovis, B. bigemina, and B. ovis [44][45][46][47][48]. In a recent study, an enzymatically active recombinant form of B. microti XP_012650559 (rBmCYP) was expressed in E. coli. The activity of rBmCYP against a fluorescent peptide was significantly inhibited by recombinant forms of the cysteine protease inhibitors cystatins 1 and 2 of Riphicephalus haemaphysaloides ticks [49]. Although R. haemaphysaloides is not a typical B. microti-transmitting tick, it has been suggested as a potential vector for this parasite in China [50]. These results coincide with the inhibition exerted by R. microplus cystatins on a B. bovis C1A cysteine protease and suggest the involvement of these enzymes in tick host-pathogen interactions [51] Interestingly, the phylogenetic analysis of C1A cysteine protease paralog profiles of piroplasmids of the Babesia, Theileria and Cytauxzoon genera corroborates the assignment of analyzed species into Clades I-VI according to their 18S rRNA gene sequences [27,52].

Serine Proteases
At least thirteen functional serine proteases and eight non functional protease homologs belonging to nine families are encoded in the B. microti genome (Table 1).
A prominent group of serine proteases is constituted by the S54 family, which consists of rhomboid proteases (ROMs). ROMs were first described in Drosophila melanogaster and later shown to be present in all kingdoms of life, fulfilling various relevant roles, including cell signaling in animals, quorum sensing in bacteria, homeostasis regulation in mitochondria, and the dismantling of adhesion complexes in apicomplexan protozoa. They are characterized by having six to seven transmembrane domains and their active site embedded in the lipid bilayer [53,54].
ROMs have been thoroughly studied in the apicomplexans Toxoplasma gondii and Plasmodium spp. The former encodes ROM1-6, according to the nomenclature defined by Dowse and Soldati, 2005 [55], all of which have, with exception of ROM2, homologs in P. falciparum. The latter parasite has four additional ROMs that are not present in T. gondii, designated ROM7-10 [56,57]. T. gondii and Plasmodium sp. ROM4 proteases were shown to cleave parasite adhesins, thus dismantling the adhesive junctions formed between the membranes of the host and parasite, a process needed for parasite internalization into the host cell [53,57]. Due to their critical role in invasion, ROMs are regarded as potential targets for therapeutic interventions against apicomplexans [58]. Indeed, two ROM4 inhibitors were shown to specifically block the P. falciparum invasion of human erythrocytes [59]. Additionally, experimental vaccine formulations based on T. gondii and Emeria tenella ROM4 were able to partially protect mice and chickens, respectively, against challenges [20,60].
In a recent study, ROM-coding genes were identified in the genomes of several piroplasmids and shown to belong exclusively to the ROM4, ROM6, ROM7, and ROM8 types. While the latter three were always present in a single copy, two to five ROM4 paralogs could be found depending on the piroplasmid lineage analyzed [61]. B. microti has two ROM4 paralogs, one of which (XP_021338238) has been misannotated as "ROM3" in GenBank (Table 1). ROM4 proteinases are found exclusively throughout the phylum Apicomplexa, which is consistent with their predicted role in invasion of the host cell, a critical mechanism for these obligate parasites [56]. ROM6, on the other hand, is the only piroplasmid rhomboid not exclusive to apicomplexans and has been shown to participate in various processes, including mitochondrial homeostasis, apoptosis, and the electron transport chain [62]. Accordingly, a mitochondrial localization was predicted for B. microti ROM6 (XP_021338360 ; Supplementary Table S3). B. microti ROM7 (XP_012650510) and ROM8 (XP_021338098) were predicted to localize in the membranes of the endoplasmic reticulum and the Golgi apparatus, respectively (Supplementary Table S3). These two types of ROMs are present in Plasmodium sp. and piroplasmids but not in other apicomplexans. Their functions are unknown but could be related to processes shared by all Aconoidasida, such as those that take place during the intraerythrocytic stage [61]. Finally, three members of the "derlin" subfamily were found in B. microti (Table 1). Derlins are catalytically inactive members of the Rhomboid Superfamily and were first described in yeast and later found in mammals and other organisms. Their function is still unclear, but it has been suggested that they could be part of a channel through which misfolded proteins are retrotranslocated from the endoplasmic reticulum to the cytoplasm prior to their ubiquitination and degradation [63].
Notably, for B. bovis, one of the ROM-encoding genes (XP_001610128) was found to be significantly higher expressed in the parasite stages present in the hemolymph of Rhipicephalus microplus ticks as compared to the stages present in bovine blood, suggesting that the role of this protease is mostly associated with the development of the parasite in the tick [64]. It remains to be analyzed whether a similar scenario takes place for the corresponding orthologs in B. microti and other piroplasmids.
In an early study, the serine protease activity of B. bovis merozoite homogenates was found to be higher in two virulent than in two avirulent strains from Australia, and thus, these proteases were postulated as virulence determinants [65]. However, in a later study, all the genes encoding for active proteases (n = 66) were shown to be present and transcribed to similar levels in the asexual blood stages of a B. bovis virulent parental strain and an attenuated strain, obtained by successive blood passages in splenectomized bovines [26]. These data suggest that the virulent/attenuated phenotype in this parasite is not related to a different peptidase gene content or to changes in the transcriptional levels of any peptidase-coding gene. To establish whether or not parasite serine or other types of proteases are virulence determinants in Babesia spp. will need further experimental evidence, but in any case, their relevance for pathogenicity is based on the vital role they probably fulfill in the parasitic lifestyle.

Metalloproteases
Metalloproteases contain a metal ion at their active site, which acts as a catalyst in the hydrolysis of peptide bonds, and are represented by at least 17 active and two non active protease homologs in B. microti (Table 1) [29].
Among them, methionine aminopeptidases (MAPs), which are present with four members in B. microti (M24A family, Table 1), take care of the N-terminal methionine excision from polypeptides, general metabolism of amino acids and proteins, and regulation processes that imply the activation and inactivation of biologically active peptides [66]. Inhibitors of MAPs significantly reduced the in vitro growth of P. falciparum, B. bovis, B. bigemina, B. caballi, and T. equi, highlighting a relevant role of MAPs in the survival of these parasites [23,67]. Moreover, B. microti-infected mice treated with MAP inhibitors reached significantly lower parasitemia levels than untreated mice [23]. Additionally, one of the B. microti MAPs (XP_012649271) was tested in mice as a potential vaccine candidate for human babesiosis. Immunization with an E. coli-expressed recombinant form of this MAP induced a Th1 immune response characterized by IgG2a antibody titers and IFN-γ production, and provided partial protection against the challenge with B. microti [68]. Although the number of boosters and protein amounts needed to achieve this effect would not be practical to apply in humans, these results suggest a potential usefulness of MAPs in future vaccine formulations against B. microti.

Threonine Proteases and the Proteasome
The proteasome is a cylindrically shaped large complex of proteins in charge of degrading intracellular proteins destined for destruction that have been tagged with polyubiquitin chains, thereby controlling many cellular processes, such as cell cycle progression and cell signaling [69,70]. All B. microti threonine proteases are proteasome constituents (seven alpha and seven beta 20S proteasome subunits, Table 1). Additionally, a metalloprotease with a proteasome regulatory function is also listed among the B. microti proteases (XP_021338577 of the M67 family), while assignment of other proteases to this structure needs experimental confirmation. Due to their vital role in cell physiology, drugs targeting proteasome functions have been proposed as therapeutics against several parasitic diseases [71][72][73]. In the case of Babesia sp., the proteasome inhibitors epoxyketones and boronic acid were shown to reduce the chymotrypsin activity of the proteasome in lysates of B. divergens in vitro cultures, leading to the accumulation of poli-ubiquinated proteins and, also, impeding parasite growth in vitro [24]. One of the epoxyketones, carfilzomib, was also assayed in B. microti-infected mice. Carfilzomib is a covalent and irreversible peptide inhibitor of the β5 subunit of the human proteasome approved for the clinical treatment of multiple myeloma [74]. Blood lysates of B. microti-infected mice treated with carfilzomib also showed the accumulation of poli-ubiquinated proteins as compared to untreated mice. Moreover, carfilzomib treatment reduced the peak parasitemia levels without apparent toxic effects in the treated mice. Although the dose required to eliminate the parasite would be toxic when applied in humans, these studies indicate that specifically targeting the B. microti proteasome would be a possible chemotherapeutic approach against this parasite [24].

Comparison between B. microti and B. bovis Functional Proteases
The genome of B. microti is the smallest among Apicomplexans and encodes 7% less genes compared to that of B. bovis. This difference is mainly due to the large vesa and SmORF multigene families present in B. bovis, which are absent in B. microti [38]. These two gene families encode for highly variable proteins that are involved in escaping the immune system of the vertebrate host and cytoadhesion [75,76]. It remains unknown whether strategies to escape effectors of the immune system exist in B. microti. However, cytoadhesion, especially affecting brain capillaries, has not been described as a major pathogenic mechanism for this parasite [2]. Other unraveled differences include the lack of spherical body proteins in B. microti, consistent with a reduced apical complex [38]. Additionally, contrary to B. bovis, B. microti does not have an oligosaccharyl transferase in charge of transferring a (NAcGlc) 2 moiety from a lipid-linked oligosaccharide to a nascent protein destined for the secretory pathway in the endoplasmic reticulum. Thus, a significant difference among B. bovis and B. microti is the lack of ability of the latter to produce N-glycosylated proteins [77].
In the present study, we hypothesized that the differences between B. bovis and B. microti include the repertoire of active proteases encoded in their genomes. By orthology searches, we observed that most B. bovis active proteases have an ortholog in B. microti (Supplementary Table S4). The lack of orthology was connected in all but two cases to the expansion of a protease-coding ancestor gene into different numbers of paralogs, which most likely took place after the separation of the most recent common ancestor (MRCA) of B. bovis and B. microti, and thus, they differentiate both species. However, the S8 family of serine proteases is present with a single member, a subtilisin-like protein (XP_001610126), only in B. bovis but is absent in B. microti. The B. bovis subtilisin-like protein gene is syntenic with orthologous genes in B. divergens, B. ovata, and Babesia sp. Xinjang (data not shown). Importantly, characterization of the subtilisin-like protein of B. divergens showed that it localizes to dense granules and contains neutralization-sensitive B-cell epitopes, consistent with a relevant role in the invasion or establishment of the parasite in the infected erythrocyte, as observed for subtilisin-1 and subtilisin-2 in P. falciparum [78][79][80]. The other case is B. microti SPP aspartic protease (XP_021338622), which is absent in B.
bovis and other Babesia spp., as mentioned above. The identification of genes absent in B. microti and present in other Babesia spp. or vice versa can allow comprehending the minimum protein dotation needed to fulfill a basic Babesia sp. life cycle, as well as to identify which proteins are associated with species-specific peculiarities and can also be exploited for differential diagnosis, therapeutic, and vaccine developments. B. bovis and B. microti share important similarities in their life cycles, namely being tick-transmitted and having an asexual reproduction stage exclusively within the erythrocytes of their vertebrate hosts. However, they differ in tick and vertebrate host species, as well as by the presence or absence of transovarial transmission in the tick. Transovarial transmission is, indeed, a trademark of the "true" babesias or Babesia sensu stricto, such as B. bovis, while those members of the Babesia genus that do not have this trait, such as B. microti, are considered Babesia sensu lato [27,28,81]. These differences are undoubtedly connected with the evolutionary history of B. bovis and B. microti, which can be clearly visualized by their phylogenetic placement into two distant clades (Clades VI and I, respectively, according to Schnittger et al., 2012 and2021 [27,28].

Non-Peptidase Homologs
At least 25 non-peptidase homologs are encoded in the B. microti genome (Table 1). A conserved protease domain can be predicted in their sequences, but they lack one or more of the catalytically relevant amino acids. Non-peptidase homologs are commonly found among living organisms and believed to have evolved from catalytically active enzymes. They have lost their catalytic capacity but developed new functions, such as competitive inhibition regulating their active counterparts or even completely new nonprotease-related activities [82,83]. An extreme case of loss of function is observed with a group of paralogs of B. microti that include three metalloproteases of the M41 family and 14 other non-protease members. Different from other non-peptidase homologs, the latter do not have a recognizable protease active site region. According to their conserved domains, their functions include the hydrolysis of nucleoside triphosphates, fusion of vesicles, intracellular transport, and proteasome regulation (Supplementary Table S1).

Conclusions and Perspectives
Proteases are attractive targets against a large number of infectious agents, since many of them are druggable and participate in essential biological processes of pathogenic virus, bacteria, protozoa, and fungi [84]. Indeed, several protease inhibitors are commercially available, and some are successfully employed in the treatment of HIV and Hepatitis C [85,86]. The use of protease inhibitors against other relevant viruses, such as dengue and SARS-CoV-2, has also been postulated [87,88].
The present study was aimed at organizing the available information of B. microti proteases and extending the array of identified peptidases encoded in its genome. This information is expected to set the stage for future research directed to understand the biology and pathogenicity of this parasite and to explore proteases as targets for developing novel therapeutic interventions. Recent advances in B. microti gene editing will permit exploring the functional relevance of selected proteases [89,90]. In addition, the application of computer-based inhibitor screening and the use of optimized pipelines to test drug efficacies using in vitro cultures and animal models allows obtaining new therapeutics against human babesiosis in a relatively short period of time [34,91,92].

Materials and Methods
The proteases of B. microti, R1 strain, presented in this study were identified by three different search approaches: (i) extracting and organizing the data available for this parasite in the MEROPS database (www.ebi.ac.uk/merops/, accessed on 1 September 2021) [29], (ii) the identification of homologs of B. bovis proteases predicted as active, as reported by Mesplet et al. (2011) [26], and (iii) the search for paralogs of B. microti proteases identified in (i) and (ii). Orthology between B. bovis and B. microti proteases was defined using a BLASTp bidirectional best hit (BBH) approach [93]. Paralogs within the B. microti genome were determined by BLASTp (blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 1 September 2021), considering a threshold E value of 0.05. Peptidase domain names and locations were obtained from the Conserved Domains database of the NCBI.
For those proteases included in the MEROPS database and predicted as active, the relevant amino acids of the catalytic site were identified using the data available at this website. For the proteases not included in the MEROPS database, alignments of B. bovis and B. microti orthologs were carried out by Clustal omega [94] (https://www.ebi.ac.uk/Tools/ msa/clustalo/, accessed on 1 September 2021), and the relevant amino acids described for B. bovis in MEROPS were manually identified for the corresponding B. microti protease. The non-peptidase homologs included those described as such in the MEROPS database. In addition, the peptidases not present in MEROPS were listed as non-peptidase homologs whenever one or more of the catalytically relevant amino acid residues at the homologous positions were missing upon alignment with the sequence of an active proteinase homolog.

Conflicts of Interest:
The authors declare no competing interests.