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

**Monica Florin-Christensen 1,2,\* , Sarah N. Wieser 1,2, Carlos E. Suarez 3,4 and Leonhard Schnittger 1,2**


**Abstract:** 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.

**Keywords:** human babesiosis; *Babesia microti*; therapeutic drugs; peptidases

## **1. 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–10].

**Citation:** Florin-Christensen, M.; Wieser, S.N.; Suarez, C.E.; Schnittger, L. In Silico Survey and Characterization of *Babesia microti* Functional and Non-Functional Proteases. *Pathogens* **2021**, *10*, 1457. https://doi.org/10.3390/ pathogens10111457

Academic Editors: Estrella Montero, Jeremy Gray, Cheryl Ann Lobo and Luis Miguel González

Received: 18 September 2021 Accepted: 6 November 2021 Published: 10 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

*B. microti* belongs to the Apicomplexa phylum and, as such, has a mandatory parasitic lifestyle that alternates between its definitive tick host and its intermediate mammalian hosts. Complex physiological processes and molecular interactions between the pathogen and host are needed for invasion, egress, parasite development in the tick stages, and migration processes that lead to the completion of the parasite life cycle and its efficient perpetuation and dissemination. Among the molecules involved in these events, parasite proteases, i.e., enzymes that catalyze proteolytic cleavages, are bound to be of paramount importance [2,11,12].

Indeed, proteases of the model Apicomplexan protozoans *Toxoplasma gondii* and *Plasmodium falciparum* have been shown to participate in several essential physiological processes, including nutrient acquisition and processing, invasion and egress from host cells, protein recycling, posttranslational processing, and signal transduction, among others [13–15]. Due to their vital roles and the fact that they show low or no identity with host-encoded peptidases, parasite proteases have been proposed as potential drug targets and/or vaccine candidates [16–20]. In the case of *Babesia* spp., the importance of peptidases for parasite survival and their potential as therapeutic targets are highlighted by several studies showing that different protease inhibitors significantly impede parasite growth in vitro and/or in vivo [21–25].

The present study aims to shortlist the proteases encoded in the *B. microti* genome by organizing the information available in the MEROPS protease database, as well as identifying additional peptidases by homology searches for paralogs within the *B. microti* genome and orthologs of previously described active proteases of *B. bovis* [26]. We also tested the hypothesis that the repertoires of functional proteases encoded in the genomes of *B. bovis* and *B. microti* differ, possibly due to the peculiarities displayed in their life cycles and their different phylogenetic placements [27,28]. The information recorded in this study can be applied to future research aimed at understanding the biology of this emergent pathogen and designing new therapeutic interventions.

#### **2. Results and Discussion**

#### *2.1. Survey of B. microti Proteases*

The present study shows that the *B. microti* genome encodes for at least 64 active proteases and 25 non active protease homologs. These proteases belong to the cysteine (*n* = 28), serine (*n* = 21), threonine (*n* = 14), aspartic (*n* = 7), and metallopeptidase (*n* = 19) types, which, in turn, are assigned to a total of 38 peptidase families (Table 1).

The classification into peptidase types refers to the nature of the nucleophile in the hydrolytic reaction, which can be the thiol of a cysteine in cysteine peptidases, the hydroxyl of a serine, or a threonine residue in serine and threonine peptidases, respectively, or water bound to aspartic acid or to a metal ion in aspartic and metallopeptidases, respectively. An additional protease group has been described, the glutamic peptidases, in which the nucleophile is water bound to a glutamic acid residue, but these enzymes are absent in Apicomplexan protozoa. Peptidases of each type are assigned into families according to sequence similarities. Non active protease homologs are characterized by bearing a conserved protease domain but lacking in the active site one or more of the critical amino acids needed for catalysis [29].

**Table 1.** Proteases belonging to the aspartic, cysteine, threonine, serine, and metallopeptidase types encoded by *B. microti*.




**Table 1.** *Cont.*







**Table 1.** *Cont.*



**Table 1.** *Cont.*



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].

#### *2.2. 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. (not shown).

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.

#### *2.3. 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–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–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].

#### *2.4. 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.

#### *2.5. 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*.

#### *2.6. 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–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].

#### *2.7. 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)<sup>2</sup> 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–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 and 2021 [27,28].

#### *2.8. 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).

#### *2.9. 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].

#### **3. 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.

The presence of transcripts and translated proteins in the blood parasite stages was evaluated in PiroplasmaDB [95] (piroplasmadb.org/piro/app, accessed on 10.09.2021) and in the proteomic database provided in Reference [31], respectively. The subcellular location of each protease was evaluated by the presence of a signal peptide (SignalP 5.0 server, www.cbs.dtu.dk/services/SignalP/, accessed on 10 September 2021) [96] and transmembrane domains [97] (TMHMM server, www.cbs.dtu.dk/services/TMHMM/, accessed on 10 September 2021) and using the localization predictor DeepLoc-1.0 [98] (www.cbs.dtu.dk/services/DeepLoc/, accessed on 10 September 2021) with the settings for eukaryotic sequences.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/pathogens10111457/s1, Table S1: M41 metalloprotease paralogue family of *B. microti*, Table S2: S59 paralogue family of *B. microti*, Table S3: Predicted subcellular localization of *B. microti* proteases, Table S4: Comparison between the repertoire of *B. microti* and *B. bovis* proteases predicted as active at least in one of either species.

**Author Contributions:** Conceptualization, M.F.-C. and L.S.; investigation: M.F.-C. and S.N.W.; original draft preparation: M.F.-C.; and writing—reviewing and editing: C.E.S., S.N.W., L.S., and M.F.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** We acknowledge the support of projects 2019-PD-E5-I102, 2019-PE-E5-I109, and 2019- PE-E5-I105 from the National Institute of Agricultural Technology (INTA, Argentina) and CRIS 2090-32000-039-000-D from ARS-USDA (USA). SNW received a doctoral fellowship from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no competing interests.

## **References**


## *Article* **Plasmepsin-like Aspartyl Proteases in** *Babesia*

**Pavla Šnebergerová 1,2 , Pavla Bartošová-Sojková 1 , Marie Jalovecká <sup>2</sup> and Daniel Sojka 1,\***


**\*** Correspondence: sojka@paru.cas.cz

**Abstract:** Apicomplexan genomes encode multiple pepsin-family aspartyl proteases (APs) that phylogenetically cluster to six independent clades (A to F). Such diversification has been powered by the function-driven evolution of the ancestral apicomplexan AP gene and is associated with the adaptation of various apicomplexan species to different strategies of host infection and transmission through various invertebrate vectors. To estimate the potential roles of *Babesia* APs, we performed qRT-PCR-based expressional profiling of *Babesia microti* APs (BmASP2, 3, 5, 6), which revealed the dynamically changing mRNA levels and indicated the specific roles of individual BmASP isoenzymes throughout the life cycle of this parasite. To expand on the current knowledge on piroplasmid APs, we searched the EuPathDB and NCBI GenBank databases to identify and phylogenetically analyse the complete sets of APs encoded by the genomes of selected *Babesia* and *Theileria* species. Our results clearly determine the potential roles of identified APs by their phylogenetic relation to their homologues of known function—*Plasmodium falciparum* plasmepsins (PfPM I–X) and *Toxoplasma gondii* aspartyl proteases (TgASP1–7). Due to the analogies with plasmodial plasmepsins, piroplasmid APs represent valuable enzymatic targets that are druggable by small molecule inhibitors—candidate molecules for the yet-missing specific therapy for babesiosis.

**Keywords:** aspartyl protease; plasmepsin; apicomplexa; piroplasmida; *Babesia*

## **1. Introduction**

Babesiosis (also known as piroplasmosis) is a malaria-like disease caused by the parasites from the genus *Babesia* of the apicomplexan order Piroplasmida. This order was initially represented by three separated lineages—*Babesia*, *Theileria*, and *Cytauxzoon*—but more recent phylogenetic analyses have indicated approximately six lineages of Piroplasmida, out of which the approximately 100 *Babesia* species are represented in at least three distinct clades [1]. Analogously to their malaria-causing relatives (genus *Plasmodium*, order Haemosporidia), *Babesia* parasites are transmitted to their vertebrate hosts via a bloodfeeding arthropod (tick), where they follow an asexual erythrocytic growth cycle. With the global distribution of ticks and their dynamically changing distribution in recent decades, babesiosis represents an important worldwide veterinary threat, and an emerging risk to humans [2].

Despite routine epidemiological surveillance, babesiosis has long been recognised as an economically important disease affecting livestock, with growing incidence in both domesticated and wildlife animals [3]. Bovine babesiosis, commonly called red water fever, is the economically most important arthropod-transmitted disease affecting cattle, causing mortalities, miscarriages, and decreased meat production [4]. Recently, increased attention has been devoted to the alarming increase in severe "dog babesiosis", caused mainly by *Babesia canis* transmitted by the ornate dog tick *Dermacentor reticulatus*, which has mosaic distribution throughout Europe [5]. Humans are accidental hosts of *Babesia*,

**Citation:** Šnebergerová, P.;

Bartošová-Sojková, P.; Jalovecká, M.; Sojka, D. Plasmepsin-like Aspartyl Proteases in *Babesia*. *Pathogens* **2021**, *10*, 1241. https://doi.org/10.3390/ pathogens10101241

Academic Editors: Estrella Montero, Jeremy Gray, Cheryl Ann Lobo and Luis Miguel González

Received: 3 August 2021 Accepted: 22 September 2021 Published: 26 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

but numerous factors indicate human babesiosis as an emerging zoonosis. The clinical features of human babesiosis are similar to malaria, and can be fatal, particularly in the elderly and in immunocompromised individuals [2]. The majority of human infections are reported from the USA, where the principal agent *Babesia microti* is the most commonly transmitted pathogen during blood transfusions [6,7]. Babesiosis is also an emerging problem in Europe, where *Ixodes ricinus*-transmitted *Babesia divergens* is considered the major causative agent of human babesiosis [8]. Incidence rates of human babesiosis are probably underestimated due to misdiagnosis of malaria in the overlapping distribution areas of these parasitic diseases [9].

Protection against bovine babesiosis is mostly based on the debatable vaccination of young cattle with attenuated parasites [10]. Treatment of human babesiosis is non-specific and relies on the combination of antimalarial drugs and antibiotics, such as atovaquone and azithromycin [6]. However, a remarkable increase in parasite resistance, together with significant numbers of relapsed immunocompromised and asplenic individuals [11], has made this widely used treatment regime less effective [12]. Alternative therapy with clindamycin and quinine is toxic and has never been tested in clinical trials. Reports of babesiosis within new geographical regions, as well as identifications of new *Babesia* species as agents of severe human disease, suggest rapid changes in *Babesia* spp. epidemiology, and make human babesiosis a serious public health concern [2]. Hence, novel approved drugs and veterinary control strategies based on *Babesia*-specific molecular targets are highly desirable [10,12].

Parasite-derived proteolytic enzymes (proteases) have been adapted for various functions connected with the parasitic lifestyle [13]. Aspartyl (aspartic, aspartate) proteases use an activated water molecule bound to one or more aspartate residues for catalysis of their peptide substrate. The peptidase database MEROPS [14] identifies five clans of aspartic proteases (AA, AC, AD, AE, and AF), each representing an independent evolution of the same active site and mechanisms. Proteases in the clan AA are either bilobed (family A1, the pepsin family) or homodimeric (all other families in the clan, including the A2 family of retropepsins) [15]. Each lobe consists of a single domain with a closed beta-barrel, and each lobe contributes one aspartate to form the active site. The monomers (retroviruses, retrotransposons, and badnaviruses) are structurally related to one lobe of the pepsin molecule. Due to their ability of specific cleavage within protein bonds, pepsin family (A1) aspartyl proteases (APs) have been evolutionary selected and adapted for unique cellular and physiological roles [16]. Apicomplexa, in contrast to metazoan parasites, have used rapid evolution of the single ancestral A1 protease, resulting in multiple AP-encoding genes being found in their genomes, as best represented by the 7 different APs of *Toxoplasma gondii* (TgASP1–7) [17] and 10 plasmepsin isoenzymes encoded by the genome of *P. falciparum* (PfPM I–X) [18]. The multiple Apicomplexan APs are phylogenetically clustered into six clearly distinguishable clades, tagged A–F. These clades reflect the function-driven evolution and various biological roles of these enzymes within the life cycle of these parasitic organisms [16,17]. Due to their essential roles and affordable selective targetability by small molecule inhibitors, APs have been considered and validated as valuable therapeutic targets [16,19].

In this work, we use the available EuPathDB and NCBI GenBank database records to identify and phylogenetically analyse the yet unrevealed plasmepsin and other apicomplexan AP homologues from selected *Babesia* and *Theileria* species. We estimate the roles of newly identified entries by their phylogenetic clustering to the six clades of apicomplexan APs, as well as their evolutionary relation to *P. falciparum* plasmepsins and *T. gondii* ASP homologues of known functions. Our estimation is supported by the dynamic expression profiling throughout selected life stages of *Babesia microti* using our previously developed *B. microti–Ixodes ricinus*–mouse acquisition model [20].

#### **2. Results and Discussion**

#### *2.1. Expression Profiling of BmASPs Indicate Their Various Roles throughout the Life Cycle of Babesia microti*

Based on the previously published phylogenetic analysis of parasite cathepsin Dlike proteases [16], we initially assessed the expression dynamics of four identified *B. microti* APs (BmASPs) to estimate their potential roles within the life cycle of *B. microti*. This was enabled by the use of our established *B. microti* acquisition model, involving BALB/c mice and *I. ricinus* tick nymphs [20]. Our results clearly show *B. microti* ASP3 (BmASP3, BmR1\_04g07350) to be the most abundantly expressed AP isoenzyme of the parasite intraerythrocytic stage (Figure 1A). This finding is in line with the phylogenetic clustering of BmASP3 to clade C of apicomplexan APs, together with the *P. falciparum* plasmepsins PfPM IX and PfPM X [21], and *T. gondii* ASP3 [22], which are the apical complex-associated proteases playing essential roles in parasite invasion and egress from the host cell (further discussed in Section 2.2.). Analogously to the clade C plasmepsins, BmASP3 is also expressed in *B. microti* stages developing in the midgut and salivary glands of *I. ricinus* nymphs (Figure 1B,C), indicating the potential involvement of BmASP3 in tick tissue invasion.

*B. microti* ASP6 (BmASP6; BMR1\_03g03850) mRNA represents the relatively most abundant AP encoding mRNA (54%) in the gut of fully fed tick nymphs, and its abundance in this tissue remains strong even 6 days post the detachment of tick nymphs from their host (24%) (Figure 1B,C). This indicates the involvement of BmASP6 in *Babesia* zygote development and resembles the mRNA expression profile of the *P. falciparum* orthologue PM VII. This enzyme is present in all sexual stages of *P. falciparum* [23], although the protein first appears in diploid stages upon fertilization, presumably due to regulated protein translation [23]. Since BmASP6 is also expressed in *B. microti* stages developing in the tick midgut 6 days post tick detachment, we hypothesise its potential role during the *B. microti* zygote–kinete stage transition, and the subsequent development of kinetes occurring in the tick midgut wall. Supported by the predominant expression of BmASP6 in the salivary glands of fully fed and detached *I. ricinus* nymphs (Figure 1B,C), BmASP6 presumably plays a role in the dissemination of kinetes to other tick tissues, including salivary glands, and subsequently in the early phase of *B. microti* sporogony. The BmASP6 *T. gondii* orthologue TgASP6 is expressed in sporulated oocyst-containing haploid sporozoites [17]. Since piroplasms do not create oocysts in tick midguts, and sporogony occurs in tick salivary glands, we speculate that BmASP6 might be synthesised as a precursor during early sporogony, and that the enzyme could be catalytically activated upon tick feeding. Mature BmASP6 might thus be involved in the release of sporozoites to tick saliva and/or sporozoite invasion of host red blood cells.

Interestingly, the clade B member *B. microti* ASP2 (BmASP2, BMR1\_01G02485) is relatively more expressed in *B. microti* blood stages. This contrasts with its orthologue PfPM VI, which is specifically expressed in the vector stages, and is involved in the correct transition of sporozoites from the oocyst [24]. BmASP2 is relatively lowly expressed in monitored tick stages (Figure 1B,C), while it appears abundant (up to 27%) in host blood (Figure 1A). Thus, BmASP2 might be already produced in gametocytes in host blood, and the protein is first produced upon fertilization in tick midguts due to regulated translation in the same manner as the clade E PfPM VII [23]. Alternatively, it might play a different role during the *B. microti* life cycle than its orthologue PfPM VI [25].

The clade D member *B. microti* ASP5 (BmASP5, BmR1\_04g05270), the PEXEL (TEXEL) cleaving *P. falciparum* PM V [26,27], and *T. gondii* ASP5 [27,28] orthologue, are expressed throughout all three selected timepoints of the *B. microti* life cycle (Figure 1). Since the exported proteins of *B. microti* were proposed to lack the PEXEL-like motifs [29] and their trafficking is most likely mediated via vesicles [30], BmASP5 supposedly plays other roles throughout the parasite life cycle, e.g., in the intracellular trafficking of proteins to specific organelles [31], in the development of gametocytes analogous to PfPM V [32], or in the secretion of proteins interacting with tick tissues during transmission. However, the role of

*B. microti* PEXEL processing should not be fully excluded, as only more complex studies on the *B. microti* secretome can fully address this point in the future.

**Figure 1.** Dynamic expression profile of aspartyl proteases (BmAPs) during the life cycle of *B. microti*. The qPCR results were obtained with cDNA templates prepared from the total RNA isolated: *A*—from mouse blood on the 6th (red asterisk) and 10th (black asterisk) days post parasite injection; *B*—from infected midgut and salivary glands of *I. ricinus* nymphs at fully fed tick stage; *C*—6 days post detachment of *I. ricinus* nymphs from the host. DPI: days post-injection; FF: fully fed stage; 6DPD: 6 days post detachment; ASP2: BmASP2 (BMR1\_01G02485, clade B); ASP3: BmASP3a (BmR1\_04g07350, clade C); ASP5: BmASP5 (BmR1\_04g05270, clade D); ASP6: BmASP6 (BMR1\_03g03850, clade E). Relative expression of individual BmASP-encoding mRNAs was counted as percentage ratios of the BmASP mRNA with the highest Ct value (100%) for each cDNA template. The obtained ratios were used to create the pie charts (100% total).

## *2.2. Data-Mining and Phylogenetic Analysis of Piroplasmid APs Reveals the Presence of Multiple AP Isoenzymes Clustering to Several Apicomplexan AP Clades* 2.2.1. Clade A

Our data search throughout available Piroplasmida sequence databases did not reveal a single AP isoenzyme that would cluster into clade A (Figure 2) represented by the digestive vacuole-residing hemoglobinolytic *Plasmodium falciparum* plasmepsins I, II, and IV (PfPM I, PfPM II, and PfPM IV), or the histo-aspartyl protease (HAP/PfPM III) [33]. This clade has apparently formed under the evolutionary demand to digest host haemoglobin. Although both apicomplexan sister orders Haemosporida and Piroplasmida are obligate intracellular parasites whose propagation strictly depends on nutrients provided by the host cell [34], their feeding mechanisms likely differ, which is also reflected in their classification: *Plasmodium* spp. intraerythrocytic stages digest a substantial amount (up to 80%) of host cell haemoglobin (Hb) [35]. Hb proteolysis releases heme, and globin serves as a source of free amino acids. During Hb digestion within the acidic environment of the food vacuole (FV), clade F plasmepsins tightly cooperate with cysteine proteases—falcipains, metalloprotease falcilysin, and other aminopeptidases that complete the degradation process [36]. The digestive plasmepsins are the most closely related isoenzymes among the 10 *P. falciparum* plasmepsins; they share 50–70% amino acid identity, and their encoding genes are located together on one chromosome [37]. Interestingly, the plasmodial species outside the primate-infecting group of *P. falciparum* encode for a single digestive plasmepsin, analogous to PfPM IV [37]. The heme moiety originating from host Hb does not

appear to be metabolised or recycled; instead, it is aggregated and polymerised to the dark pigmented hemozoin (Hz) (the name of the order Haemosporida) [38–40].

In contrast to *Plasmodium* spp., piroplasms degrade little if any Hb during host erythrocyte infection, and no pigmented Hz is formed [41–43]. Knowledge of the exact routes of nutrient uptake and processing remains rather unclear and is mainly based on electron microscopy observations of trophozoites. The FV formation via cytostome has been observed in *Theileria* spp., where the optical density of this vacuole indicated potential Hb presence [44]. In addition, *T. equi* possibly uses another tubular feeding structure that is formed via the invagination of an erythrocyte plasma membrane enclosing only a minor amount of ingested cytoplasmic material [45]. *Babesia* does not possess a cytostome; thus, the true FV that emerges from its constriction cannot be formed. The ovoid bodies in *Babesia* trophozoites, initially considered to be true FVs, have since been recognized as invaginations of host cell cytoplasm [46]. In *B. microti*—the basal piroplasmid species (*Babesia microti*-like group)—the obscure coiled structure that protrudes from the parasite into the host cytoplasm has been observed and speculated to contain digestive enzymes. Similarly, possible haemoglobin-containing vesicles were also observed from *B. divergens* trophozoites, indicating—yet not confirming—the endocytic uptake of hemoglobin from the cytoplasm of host erythrocytes [47]. The role of hemoglobin invaginations in trophozoites of different species of *Babesia* remains unclear, but could be relevant to the biology of piroplasms, e.g., by limited digestion of hemoglobin providing the essential source of heme to these heme auxotrophic organisms. Moreover, the special organelle of *B. microti* allegedly contains ferritin, which may be used as a nutrient source for *Babesia* [46]. Although piroplasms do not encode direct homologues of Hb-digesting plasmepsins, they do encode several papain-like cysteine proteases [36]. However, these orthologues of hemoglobinolytic *Plasmodium* spp. falcipains have not yet been subjected to biochemical and functional characterization that would reliably validate their contribution to Hb digestion in erythrocyte-residing piroplasms. Thus far, the only evidence of their essential role in the survival of *Babesia* spp. parasites has come from the observation of a hampering effect on *B. bovis* erythrocyte invasion and in vitro replication by cysteine protease inhibitors [48], and from the lowered parasitemia observed in *B. ovis* erythrocyte cultures exposed to antibodies against ovipain-2, the *B. ovis* orthologue of falcipain-2 from *P. falciparum* [49]. The same applies to the falcilysin protein family represented by one or two isoenzymes in *Babesia* and *Theileria*, respectively [36]. A homologue of the *P. falciparum* heme detoxification protein (HDP) has already been identified in *Babesia* and *Theileria*. However, this homologue supposedly plays a different role than in *Plasmodium* spp., since piroplasms do not form the hemozoin pigment [36]. Overall, the molecular basis of feeding and catabolic metabolism in piroplasms, and the involvement of individual aspartyl and cysteine proteases in protein digestion, remain unclear, and should be addressed by future experimental studies.

**Figure 2.** Phylogenetic tree of apicomplexan pepsin family aspartyl proteases (AP). Data-mined APs from selected piroplasmid species cluster into previously determined clades A–F together with their apicomplexan homologues. The image displays the unrooted maximum likelihood phylogenetic tree of 106 selected apicomplexan AP sequences, including those from the basal free-living species *Vitrella brassicaformis* and *Chromera velia;* clades A–F are tagged and highlighted with colours. Sequences were retrieved from EuPathDB and GenBank, and the source organisms are indicated (upper left). Nodal supports were calculated from 1000 bootstrap replicates; those lower than 50 are not depicted. For better orientation within the tree, the 7 *Toxoplasma gondii* ASPs and the 10 *Plasmodium falciparum* plasmepsins are highlighted (arrows: yellow-green and magenta dots, respectively). Data-mined piroplasmid APs are tagged with dark blue dots. Note: The branch leading to clade D was shortened to 25% of its original length for optimal image display. Multiple alignment data are accessible as an online supplementary file at Mendeley Data, link: http://dx.doi.org/10.17632/ds3f2j32ny.1, accessed on 23 September 2021.

## 2.2.2. Clade B

Clade B is an independent cluster originally recognized as the ancestral group of apicomplexan APs by Jean et al. [50]. Some members of this group—for example, the proteases of *T. gondii* (TgASP2 and TgASP4), *Eimeria tenella* (eimepsin), *Cryptosporidium parvum* (EAK89992), and *Theileria annulata* (TA02510)—are predicted to be GPI-anchored, while plasmepsins VI and VIII lack a sufficiently long hydrophobic region at the C-terminus for the GPI anchor prediction [17]. Our analyses confirm the phylogenetic relation of clades B and E (discussed below) and highlight the existence of three subclusters within the clade B. In addition to the basal subgroup represented by the two *Chromera velia* isoenzymes, the two other independent subgroups are represented by plasmepsin VI (PM VI)/TgASP2 and plasmepsin VIII (PM VIII)/TgASP4, respectively (Figure 2). Interestingly, all analysed piroplasmid species encode for a single clade B AP orthologue clustering with the first listed group of PM VI/TgASP2. Analogously to PM VI, no GPI anchor is predicted in this enzyme (data not shown). Clade B member plasmepsins, studied mostly in the rodent malaria parasite *Plasmodium berghei*, are primarily expressed in oocysts and sporozoites the parasite transmission stages [24,51]. They play a role in midgut sporozoite development. Functional genomic analyses involving gene disruptions did not confirm the essential roles of these enzymes for the blood stages of malaria [24,52], but the *pm vi* gene knockout did affect sporozoite development from oocysts, resulting in an unsuccessful transmission of the parasite through the mosquito vector [24]. The role of the plasmepsin VI piroplasmid orthologues might thus be connected to the penetration and migration of piroplasmid stages through tick tissues, offering a suitable target to develop transmission-blocking therapy. However, our dynamic expression profiling (Figure 1) analysis controversially revealed a higher abundance of BmASP2-encoding mRNA in the blood stages of *B. microti* in comparison to the tick tissue isolates. This might be explained by the expression of clade B piroplasmid APs already encoding mRNA in gametocytes in host blood, and by regulated protein translation of the clade B AP enzyme later in tick tissues (upon fertilization).

The *P. berghei pm viii* gene knockout lineage developed in the mosquito midgut, but it showed a limited ability of the parasites to egress from oocysts. This was accompanied by a drastic decrease in the number of salivary gland and haemolymph sporozoites with a defect in gliding motility, leading to the block of transmission to hosts [51]. Interestingly, the egress phenotype mirrors that seen with PfPM X knockout lineages in intraerythrocytic parasites (see below). We propose that the reason piroplasms do not possess a direct homologue of the PfPM VIII/TgASP4 subclade is because they do not create—and do not need to egress from—oocysts in the tick gut tissue.

#### 2.2.3. Clade C

In recent years, Clade C of apicomplexan APs has gained a lot of attention, as this group of enzymes have been demonstrated as regulators of invasion and egress of host cells. This clade comprises *P. falciparum* plasmepsins IX and X (PfPM IX and PfPM X), and their *T. gondii* homologue TgASP3 [21,22]. Our analysis clearly identified two different isoenzymes clustering to two groups of clade C piroplasmid APs. While the first piroplasmid group (here termed *Babesia* ASP3b) firmly clusters together with PfPM X (bootstrap value of 97), and its members might thus be considered as the PfPM X direct orthologues, the sister relationship of the second piroplasmid cluster (here termed *Babesia* ASP3a) and the PM IX + TgASP3 group is not well supported (bootstrap value <50). Clade C APs are believed to undergo self-catalytic activation, upon which they proteolytically process/activate a vast number of protein precursors associated with secretory organelles of the apical complex (AC) [19,22,53]. Thus, PfPM IX and PfPM X are the major self-activating maturases standing at the top of the proteolytic machinery regulating *P. falciparum* cell invasion and egress [54]. Although the mechanism of invasion and egress employing the secretory and non-secretory parts of the AC is relatively conserved among apicomplexan parasites, species-specific alterations of the generic mechanisms have evolved [55,56]. *Babesia* and *Plasmodium* replicate solely inside host erythrocytes, and the invasion process involves

initial attachment of the parasite to the host cell via a receptor-ligand-mediated interaction [57]. Both parasites interact with glycophorin receptors on the RBC surface, although their surface adhesins differ [58,59]. *Plasmodium* merozoites attach via merozoite surface proteins (MSPs), erythrocyte-binding ligands (EBLs), and reticulocyte binding-like homologues (RHs)—protein families [60] that undergo proteolytic shedding by PfPM X [53]. Proteins of the merozoite surface antigen (MSA) family that are found on the *Babesia* merozoite surface also undergo processing by a yet unidentified sheddase [61]. This role supposedly might be played by the newly identified *Babesia* ASP3b isoenzyme. Similarly, *Theileria*'s surface coat is shed by an undescribed protease, suggesting an identical role of the *Theileria* ASP3b orthologues [62]. Initial attachment of the parasite to the host erythrocyte is followed by the reorientation of the apical tip towards the host cell, and the establishment of the moving (tight) junction (MJ) [63]. The components of MJ—such as the parasite membrane-associated apical membrane antigen 1 (AMA1) and the host cell membrane-anchored rhoptry neck protein complex serving as a ligand structure—are also conserved among piroplasms [63,64]. In *P. falciparum*, AMA1 is directly cleaved by PfPM X, and subsequently shed via PfPM X-activated integral membrane serine protease subtilisin 2 (SUB2) [21]. Since AMA1 of *Babesia* parasites is also proteolytically processed at several positions [65–67], we propose that ASP3b might be analogously involved in AMA1 maturation. However, this concept needs to be experimentally validated because other proteases, such as the intramembrane-cleaving rhomboids that have been identified from *Babesia* [68], cannot be excluded as *Babesia* AMA-1 sheddases working in the same manner as shown in *T. gondii* tachyzoites [69]. Moving junction protein orthologues are most likely dispensable for the non-motile merozoites of *Theileria* that enter the host cell in any orientation via a passive process known as zippering [64,70]. Thus, processing of *Theileria* AMA1 via the *Theileria* ASP3b homologues might be crucial in different motile invasive stages [71]. Upon its internalization within the host cell, *P. falciparum* remains surrounded by parasitophorous vacuole membrane (PVM) and undergoes schizogony. Later, exoneme-specific serine protease subtilisin 1 (PfSUB1), which is cleaved both autocatalytically and by PfPM X, is in charge of PVM lysis, erythrocyte plasma membrane poration, and red blood cell rupture [72]. Only minutes prior to the merozoite's egress from the erythrocyte, PfSUB1 is secreted into the parasitophorous vacuole (PV) space, where it cleaves the pseudoprotease serine repeat antigen 5 (SERA5)—the negative regulator of *P. falciparum* egress [72,73]—as well as other PV-resident proteins important for egress [74]. Additionally, PfSUB1 initiates the primary processing of merozoite surface protein complex (MSP1/6/7) [75,76]. Piroplasmid genomes encode for a single subtilisin-like protease that is believed to represent the direct orthologue of PfSUB1. This protease presumably undergoes analogous post-translation processing during secretory transport [77]. However, the analogy in identical protein substrate cleavage by *Plasmodium* and piroplasmid SUB1 proteases [78], as well as the involvement of the activated SUB1 enzyme in *Babesia* and *Theileria* parasite egress from erythrocytes, remains to be experimentally confirmed. Notably, PfPM X mRNA is also transcribed in gametocytes, while the protein can be found in gametes, zygotes, and ookinetes. In the late stage of ookinete development, PfPM X cleaves the cell-traversal protein of ookinetes and sporozoites (CelTOS) [21]. This processing is necessary for ookinetes to pass through arthropod cells to the site of oocyst formation. Moreover, the disruption of CelTOS abolishes liver infection by *P. berghei* sporozoites [79]. As CelTOS is conserved among apicomplexan parasites [80], we propose that its processing by clade C APs plays an important role during the transmission of piroplasmid species through their tick vectors.

PfPM IX, the second *P. falciparum* protease member of the apicomplexan clade C APs, has been previously confirmed to be a key maturase involved in erythrocyte invasion by *P. falciparum* merozoites [21]. However, more recent contributions have speculated as to its multifaceted role throughout the malaria parasite life cycle [53]. PfPM IX processes rhoptryassociated protein 1 (RAP1) and rhoptry neck protein 3 (RON3), which are released into the PV during invasion and later aid the parasite's development within the PV [81,82]. The

translation initiation of PM IX has been suggested during *P. berghei* gamete development in mosquitos, where PM IX cleaves merozoite thrombospondin-related anonymous protein (MTRAP), enabling the release of gametes from host erythrocytes [83]. When mosquitoes were infected with PbMTRAP knockout gametocytes, the oocyst formation was aborted. In addition, these gametocytes were unable to form ookinetes (the motile form of zygote) in vitro. Particularly, TRAP is a conserved family of proteins that are involved in the gliding motility of apicomplexan parasites, and are also found in *Babesia* [68,83]. This suggests that piroplasmid APs branching together with PM IX might be involved in the development of zygotes and kinetes within tick tissues, as indicated by the expression of BmASP3 mRNA following the detachment of tick nymphs from the host (Figure 1C).

#### 2.2.4. Clade D

Clade D is the most derived group of apicomplexan APs within our phylogenetic analysis (Figure 2); it consists of distant relatives of the human β-site amyloid precursor protein-cleaving enzyme (BACE) [27]—the major beta secretase generating amyloid-β peptides in the neurons. BACE is also responsible for the generation of the amyloid-β peptides that aggregate in the brains of Alzheimer's patients and, thus, represents a valuable target for drug development [84]. These proteases have a long C-terminal extension with a trans-membrane domain that serves for their anchoring to membranes. The most studied member of this clade is PfPM V—a *P. falciparum* endoplasmic reticulum (ER)-resident protease that has been demonstrated to cleave proteins containing the *Plasmodium* export element RxLxE/Q/D (PEXEL) motif [26]. These effector proteins are then secreted through the PV surrounding the parasite during the intracellular infection and multiplication to host erythrocytes [26,27]. PfPM V has also been shown to be essential for gametocyte development [32]. Thus, specific inhibition of PfPM V activity appears to be a highly convenient therapy targeting both the asexual and sexual stages of malaria [32,85]. Interesting findings have been made with PM V's corresponding *T. gondii* orthologue TgASP5—a Golgi-resident AP that processes TEXEL (*T. gondii* PEXEL-like) motif containing dense-granule proteins (GRAs) [86]. Analogously to its orthologue PfPM V, TgASP5 plays an important role during intracellular parasite survival and multiplication [87], as well as modulation of host cell responses [28]. Another *T. gondii* isoenzyme, TgASP7, which also clusters within clade D, is not expressed in the tachyzoite stage, and its functional role remains unknown [17].

Our phylogenetic analysis revealed piroplasmid clade D APs to be a single subgroup that clusters alongside *Cryptosporidium* and *Plasmodium* PM V. Thus, we propose that they play a similar role in the secretory pathway of piroplasms—the cleavage of PEXEL-like containing proteins in the ER of the cell—although we expect significant differences due to the different strategy of piroplasms in persisting inside infected cells. Shortly upon host cell invasion, *Babesia* and *Theileria* parasites are surrounded by a PVM that is derived from the host cell cytoplasmic membrane, but unlike the *Plasmodium* and *Toxoplasma* PVM, it starts its disintegration after parasite internalization [31,88,89]. It has been discussed previously that the PVM breaks down either because piroplasms are not able to incorporate lipids into the PVM during parasite intracellular growth, or because they have developed this strategy as a more convenient way to transport effector proteins into infected host cells in order to alter their morphology and physiology [90]. The persistence of piroplasms directly in host cell cytoplasm is thus different to *P. falciparum* residing in the PV surrounded by the PVM. In the malaria parasite, the PEXEL-containing protein precursors are cleaved by the ER-resident PM V immediately upon their translation. This facilitates them for export to the PV, PVM, and all the way to the host cell [26]. The PEXEL-motif-containing exportome of *P. falciparum* is estimated to include ~463 proteins [26]. On the other hand, TgASP5 recognizes only several identified TEXEL-motif-harbouring proteins, but also appears to be important for the trafficking of other (non-TEXEL)-secreted proteins during intracellular infection with *T. gondii* tachyzoites [91]. Both *Plasmodium* and *Toxoplasma* use sophisticated protein transporting complexes PTEX (plasmodium translocon of exported proteins) and the MYR translocon, respectively, which are incorporated into the PVM [92,93]. Unsurprisingly, piroplasms do not encode for components of the *Plasmodium* translocon machinery [29], with the single exception of the PTEX complex protein component HSP101, which is expressed across all piroplasms [62].

The PEXEL-like motif (PLM) has also been recognized in various piroplasmid proteins. This supports the above-given hypothesis of functional analogy between the piroplasmid clade D APs and plasmodial PM V enzymes [94]. PLM has been recognized in *B. bovis* variant erythrocyte surface antigens (VESAs) [95] involved in erythrocyte adhesion and antigenic variation of the red blood cell surface—sophisticated parasite-driven mechanisms enabling the infected erythrocytes to evade host immune responses [96]. VESA1-like proteins have also been described from the genomes of other piroplasmid species—the *Babesia* sensu stricto group [97], and the *B. microti*-like group [98]—indicating that VESA processing by clade D APs might be an immunoevasive strategy shared across *Babesia* species. Small open reading frame protein families (SmORFs) and spherical body protein-2 protein family members (SBP2) also contain PLM [94]. While SmORFs are also involved in erythrocyte adhesion [95,96], SBP2 proteins concentrate under the red blood cell cytoplasmic membrane and are believed to alternate the red blood cell surface [99]. Although PLM has not yet been detected in *Theileria* [62,97], the apparent conservation of clade D APs within the phylum suggests the preserved mode of protein processing across Piroplasmida. Further studies on clade D APs should further elucidate the exact role of these enzymes in the trafficking mechanisms of *Babesia* spp. effector proteins affecting infected host cells.

#### 2.2.5. Clade E

This group of apicomplexan APs is clearly determined by our phylogenetic analysis as the sister clade to clade B, which includes PM VII and TgASP6 (Figure 2). PM VII is thus sometimes classified together with clade B PM VI and PM VIII as the transmission-stage plasmepsins [37]. PM VII is produced in the ookinetes of *P. falciparum*, supposedly upon fertilization, as PM VII is not detected in *P. falciparum* gametocytes [23]. Its role remains unknown—it has been proven to be disposable for the parasite, as the PM VII knockout had no effect on any stage of the *P. berghei* life cycle [23]. However, the authors of this study also note that redundancy exists among transmission-stage-expressed plasmepsins—and especially clade B member PM VIII, which shares a similar expression profile with PM VII, may thus compensate for its loss of function [100]. This would probably not be the case for piroplasms, which do not possess direct homologues of PM VIII, but produce a single clade E AP orthologue. If the clade B PM VI orthologues play a role in regulating sporogony in tick salivary glands, the remaining function for clade E piroplasmid APs might thus insist in parasite invasion of tick tissues upon zygote formation. This is in accordance with the obtained expression profile of BmASP6 in stages developing within the vector midgut (Figure 1). However, some role of BmASP6 and other clade E piroplasmid APs in sporogony should not be fully excluded, as BmASP6 mRNA is also present in the salivary glands of infected tick nymphs (Figure 1).

#### 2.2.6. Clade F

This is a diverse group of apicomplexan APs first identified in our 2016 contribution [16], and now once more confirmed by our current phylogenetic analysis (Figure 2). This group is represented by TgASP1 and its coccidian homologues, as well as some APs of gregarines and free-living basal Apicomplexa. TgASP1 is an enzyme associated with the secretory pathway in non-dividing cells, which re-localizes in close proximity to the nascent inner membrane complex (IMC) of daughter cells during replication [17]. However, it's role is non-vital for *T. gondii* [101], and Haemosporida and Piroplasmida genomes apparently do not encode any of the clade F homologues as they might have lost these disposable proteases as a part of their adaptation to the parasitic lifestyle and the evolution of specific intracellular multiplication mechanisms differing from *T. gondii* endodyogeny [102].

#### **3. Conclusions 3. Conclusions**

endodyogeny [102].

*Babesia* and other Piroplasmida encode for several pepsin (cathepsin-D-like) family AP isoenzymes. Their phylogenetic relation to malarial plasmepsins and analogous enzymes from other apicomplexan parasites enable the prediction of their various roles within the lifecycle of these erythrocyte-infecting parasites (Figure 3). These roles are associated with their different protein structures, time-expression profiles, and intracellular localization. As these enzymes have been long considered and recently validated as great therapeutic targets for malaria, they are worthy of scientific attention when proposing novel therapeutic strategies for babesiosis (piroplasmosis). *Babesia* and other Piroplasmida encode for several pepsin (cathepsin-D-like) family AP isoenzymes. Their phylogenetic relation to malarial plasmepsins and analogous enzymes from other apicomplexan parasites enable the prediction of their various roles within the lifecycle of these erythrocyte-infecting parasites (Figure 3). These roles are associated with their different protein structures, time-expression profiles, and intracellular localization. As these enzymes have been long considered and recently validated as great therapeutic targets for malaria, they are worthy of scientific attention when proposing novel therapeutic strategies for babesiosis (piroplasmosis).

nascent inner membrane complex (IMC) of daughter cells during replication [17]. However, it's role is non-vital for *T. gondii* [101], and Haemosporida and Piroplasmida genomes apparently do not encode any of the clade F homologues as they might have lost these disposable proteases as a part of their adaptation to the parasitic lifestyle and the evolution of specific intracellular multiplication mechanisms differing from *T. gondii*

*Pathogens* **2021**, *10*, 1241 12 of 20

**Figure 3:** The role of *Babesia* plasmepsin-like APs marked within the generic life cycle of *Babesia* parasites. Sporozoites are transmitted from the salivary glands of an infected tick to the vertebrate host bloodstream during tick feeding. They invade red blood cells, where they start asexual multiplication (merogony). The cyclic egress and invasion of host red blood cells by haploid merozoites represents the intraerythrocytic cycle, causing the symptoms of babesiosis in vertebrates by time the first sexual stages (gametocytes) occur in the bloodstream. When another naïve tick feeds on the infected host, gametocytes are ingested with the blood meal, mature, and produce gametes in the tick gut lumen (gamogony). Newly developed gametes fuse in a zygote, which passes through the peritrophic matrix into tick gut epithelial cells. Meiotic division and subsequent mitosis give rise to primary kinetes, which invade and multiply in different tick tissues. Newly emerged secondary kinetes undergo multiplication in salivary glands (sporogony). Upon the tick metamorphosis, motile sporozoites are transmitted to the host with the blood meal (transstadial transmission). Secondary kinetes of *Babesia* sensu stricto species are capable of invading and multiplying within adult female ovaries and can be transmitted to the larval progeny (transovarial transmission), while the *Babesia microti*-like group parasites are transmitted solely transstadially [25]. The colour-marked text notes indicate the positions in the *Babesia* life cycle positions where the five aspartyl proteases ASP2, ASP3a, ASP3b, ASP5, and ASP6 supposedly play their herein-deduced roles. White background: **Figure 3.** The role of *Babesia* plasmepsin-like APs marked within the generic life cycle of *Babesia* parasites. Sporozoites are transmitted from the salivary glands of an infected tick to the vertebrate host bloodstream during tick feeding. They invade red blood cells, where they start asexual multiplication (merogony). The cyclic egress and invasion of host red blood cells by haploid merozoites represents the intraerythrocytic cycle, causing the symptoms of babesiosis in vertebrates by time the first sexual stages (gametocytes) occur in the bloodstream. When another naïve tick feeds on the infected host, gametocytes are ingested with the blood meal, mature, and produce gametes in the tick gut lumen (gamogony). Newly developed gametes fuse in a zygote, which passes through the peritrophic matrix into tick gut epithelial cells. Meiotic division and subsequent mitosis give rise to primary kinetes, which invade and multiply in different tick tissues. Newly emerged secondary kinetes undergo multiplication in salivary glands (sporogony). Upon the tick metamorphosis, motile sporozoites are transmitted to the host with the blood meal (transstadial transmission). Secondary kinetes of *Babesia* sensu stricto species are capable of invading and multiplying within adult female ovaries and can be transmitted to the larval progeny (transovarial transmission), while the *Babesia microti*-like group parasites are transmitted solely transstadially [25]. The colour-marked text notes indicate the positions in the *Babesia* life cycle positions where the five aspartyl proteases ASP2, ASP3a, ASP3b, ASP5, and ASP6 supposedly play their herein-deduced roles. White background: part of *Babesia* life cycle occurring in the vertebrate host; grey background: part of *Babesia* life cycle occurring within the tick vector.

#### **4. Materials and Methods**

#### *4.1. B. microti Propagation in Mice*

*B. microti* (Franca) Reichenow (strain Peabody mjr) was obtained from ATCC (ATTC® PRA-99™, USA). Two BABL/c mice, supplied by Charles River Laboratories (VELAZ), were intraperitoneally injected with 150 µL of *B. microti-*infected murine blood (50% parasitemia, 800 <sup>×</sup> <sup>10</sup><sup>6</sup> of infected red blood cells). One mouse was kept under general anaesthesia, and the blood was collected from the carotid artery into sodium citratephosphate-dextrose solution (ratio 1:25, Sigma-Aldrich) on the 6th day post injection (6DPI) according to 50% parasitemia. Similarly, the blood was obtained from the second mouse when parasitemia dropped to 5% on the 10th day post injection (10DPI). The murine experiment was performed repeatedly to confirm the results of expression profiling. All laboratory animals were treated in accordance with the Animal Protection Law of the Czech Republic No. 246/1992 Sb., ethics approval No. 25/2018.

#### *4.2. RNA Isolation from Tick Tissues and Murine Blood Cells*

Pathogen-free *Ixodes ricinus* nymphs were obtained from the in-house breeding facility of the Institute of Parasitology, BC CAS, Ceske Budejovice, Czech Republic. Twenty individuals were placed on one *B. microti*-positive BABL/c mouse in the acute phase of infection (1–4DPI, Figure 1A) and allowed to feed until repletion (around 72 h). The fully fed (FF) nymphs were collected and surface-sterilized by washing in 3% H2O2, 70% ethanol, and distilled water (each wash 30 s). Ticks were separated into two groups: first group of 10 nymphs was dissected immediately (FF stage), while the other 10 individuals were kept at room temperature in a humid chamber until they were dissected on the 6th day post detachment (6DPD). Dissection of tick tissues (salivary glands and midguts) was performed under a stereomicroscope (Olympus) on wax dishes with diethyl pyrocarbonate (DEPC)-treated cold phosphate-buffered saline (PBS), and then transferred into RA1 buffer (NucleoSpin RNA II Kit, Macherey-Nagel, Düren, Germany) supplemented with β-mercaptoethanol (Sigma-Aldrich). Prior to the extraction, the collected tissues were homogenised using an insulin syringe. Total RNA was extracted from the pool of midguts and salivary glands originating from 10 individual ticks via the NucleoSpin RNA II Kit, following the protocol provided by the manufacturer (Macherey-Nagel, Düren, Germany). Murine blood total RNA was isolated using the previously described protocol [103]. Samples were collected from two timepoints at 6 and 10 DPI. The quality and concentration of total RNA were checked by gel electrophoresis and determined using a NanoDrop UV spectrophotometer (Thermo Fisher Scientific; Waltham, MA, USA), and RNA samples were stored at −80 ◦C.

#### *4.3. Quantitative RT-PCR*

Reverse transcription was performed from 1 µg of total RNA isolate using the Transcriptor High-Fidelity cDNA Synthesis Kit (Roche Diagnostics GmbH; Mannheim, Germany). The resulting cDNA was used as a template for the quantitative real-time PCR (qRT-PCR) using a LightCycler 480 (Roche Diagnostics GmbH), the Fast Star Universal SYBR Green Master Mix (Roche Diagnostics GmbH), and according primer pairs BmASP2 forward: 50 - TCCGGCGTCTATTGAAGAGT-30/BmASP2 reverse: 50 -TGAACCGGTGTCAAAA ACAA-30 ; BmASP3 forward: 50 -GGAAGCTTGGGGAGTCTGTA-30/BmASP3 reverse: 50 - TGTGCTCCCTGTGTCGAATA-30 ; BmASP5 forward: 50 -GCCCAAACACCACCAACTAT-3 0/BmASP5 reverse: 50 -CACCAAATGCGAGATACACG-30 ; BmASP6 forward: 50 -GATTGG GCTTCCCAAACAC-30/BmASP6 reverse: 50 -ATCCGCCAGTTGAATCTTTG-30 . All qRT-PCR amplifications were performed in technical triplicates. Relative expressions were calculated using the mathematical model of the ∆Ct method [104] and normalized to *B. microti* actin (GenBank XM\_012791652; primers—BmActin forward: 50 -GGCCTACTCACAGCCCT TTA-30/BmActin reverse: 50 -ACAGGGTTGTAGAGTGTTGGTT-30 ). To express the representation of all isoenzymes as a percentage per time point, the ASP with the highest Ct value was set as 100%, and the values for other ASP isoenzymes were accordingly

recalculated. The algorithm applied later adjusted the values to fit the total of 100% in a pie chart.

#### *4.4. Phylogenetic Analysis*

The dataset used for the phylogenetic analysis comprised 106 AP protein sequences of representatives from the phylum Apicomplexa, and related *Vitrella* and *Chromera* spp. Clade F involving digestive plasmepsins served as an outgroup. All sequences were retrieved from either GenBank or EuPathDB using the blastp and tblastn BLAST algorithms and an Evalue cutoff of 10−<sup>5</sup> . Alignment was constructed in Geneious Prime 2020.1.2. using MAFFT v7.017 [105] with the default parameters for the gap opening penalty (1.53) and the offset value (0.123). The protein sequences were crosschecked for the presence of DTG/DTG or DTG/DSG aspartic protease motifs. Poorly aligned N- (signal peptide included) and C-termini were manually trimmed, which resulted in the final alignment comprising 324 amino acid positions. The phylogenetic tree was reconstructed via maximum likelihood (ML) method in IQ-TREE v1.6.12 [106], using the WAG + F + I + G4 model selected by ModelFinder [107]. Bootstraps were based on 1000 replicates. The tree was visualized in Geneious Prime v2019.0.4 and graphically modified in CorelDRAW graphic suite 2020.

**Supplementary Materials:** The multiple alignment data used to construct the maximum likelihood tree in Figure 2 are accessible online at Mendeley Data, link: http://dx.doi.org/10.17632/ds3f2j32ny.1, accessed on 23 September 2021.

**Author Contributions:** P.Š. and D.S. designed and performed the experiments, performed the analyses, and wrote the original draft; P.B.-S. supervised phylogenetic analysis, its visualization and interpretation, and reviewed and edited the manuscript; M.J. designed the expression profiling experiments, performed qPCR, and reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was primarily supported by the grant Czech Science Foundation (GA CR) project No. 20-05736S, and the ERDF/ESF Centre for research of pathogenicity and virulence of parasites (No.CZ.02.1.01/0.0/0.0/16\_019/0000759). P.Š. was additionally supported by an internal grant from the University of South Bohemia GAJU No. 120/2021/P. D.S. and P.B-S. were additionally supported by the grant number LTAUSA17201 from the Ministry of Education, Youth, and Sports of the Czech Republic.

**Institutional Review Board Statement:** All animal experiments were carried out in accordance with the Animal Protection Law of the Czech Republic No. 246/1992 Sb., ethics approval No. 34/2018, and protocols approved by the responsible committee of the Institute of Parasitology, Biology Centre of the Czech Academy of Sciences.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data are either contained within the manuscript and supporting information or available from the corresponding author on reasonable request.

**Acknowledgments:** We would like to thank Luïse Robbertse, David Hartmann, and Dominika Reichensdörferová for their initial help with *Babesia* parasites, and their previous involvement in apicomplexan-AP-related topics. We would also like to thank to the head of the laboratory of Vector Immunology, IoP BC CAS, Ceske Budejovice, Petr Kopáˇcek, for his continuous support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Pathogens* Editorial Office E-mail: pathogens@mdpi.com www.mdpi.com/journal/pathogens

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-4471-7