Skip to main content

Experimental Babesia rossi infection induces hemolytic, metabolic, and viral response pathways in the canine host



Babesia rossi is a leading cause of morbidity and mortality among the canine population of sub-Saharan Africa, but pathogenesis remains poorly understood. Previous studies of B. rossi infection were derived from clinical cases, in which neither the onset of infection nor the infectious inoculum was known. Here, we performed controlled B. rossi inoculations in canines and evaluated disease progression through clinical tests and whole blood transcriptomic profiling.


Two subjects were administered a low inoculum (104 parasites) while three received a high (108 parasites). Subjects were monitored for 8 consecutive days; anti-parasite treatment with diminazene aceturate was administered on day 4. Blood was drawn prior to inoculation as well as every experimental day for assessment of clinical parameters and transcriptomic profiles. The model recapitulated natural disease manifestations including anemia, acidosis, inflammation and behavioral changes. Rate of disease onset and clinical severity were proportional to the inoculum. To analyze the temporal dynamics of the transcriptomic host response, we sequenced mRNA extracted from whole blood drawn on days 0, 1, 3, 4, 6, and 8. Differential gene expression, hierarchical clustering, and pathway enrichment analyses identified genes and pathways involved in response to hemolysis, metabolic changes, and several arms of the immune response including innate immunity, adaptive immunity, and response to viral infection.


This work comprehensively characterizes the clinical and transcriptomic progression of B. rossi infection in canines, thus establishing a large mammalian model of severe hemoprotozoal disease to facilitate the study of host-parasite biology and in which to test novel anti-disease therapeutics. The knowledge gained from the study of B. rossi in canines will not only improve our understanding of this emerging infectious disease threat in domestic dogs, but also provide insight into the pathobiology of human diseases caused by Babesia and Plasmodium species.

Peer Review reports


Hemoprotozoal parasites cause significant morbidity and mortality in both humans and animals [1,2,3,4,5,6]. Parasites of the genus Babesia are the most common blood transfusion-transmitted pathogen in the United States and are the second most ubiquitous tick-borne pathogen in animals [7,8,9]. World-wide, human babesiosis is emerging as the tick vector expands its territory and encounters more people than ever before [4, 9].

Babesia rossi, a common parasite in wild jackals, is known to cause the most severe disease of all the Babesia species infecting canines [10], and is a leading cause of infectious morbidity and mortality among the susceptible domestic dog population in South Africa [11]. Over the course of several decades, observational studies of naturally-occurring B. rossi infections in canines have led to clinical and pathological descriptions of the disease [12, 13]. For one, clinical presentation and laboratory measurables are used to objectively classify the disease as uncomplicated or complicated. Predictors of severity and poor outcome are largely known. The host response to infection consists of pro-inflammatory cytokine storm [14]. Additionally, some individual organ systems have been described in terms of pathology; others are currently being studied [15]. Understanding the pathogenesis of severe babesiosis in a large mammal may provide insights into human severe malaria, which shares key pathophysiologic features such as hemolysis and inflammation that lead to multiple organ dysfunction [16, 17]. However, conclusions drawn from clinical case series are unreliable because the days that have elapsed from initial infection to presentation to clinic are unknown. As such, the clincal and transcriptomic time course of the disease was unable to be characterized.

To advance our understanding of the temporal dynamics of the host response to B. rossi infection, we performed a controlled infection with identical replicates in which both the infectious inoculum and time of infection were known. We inoculated purpose-bred canines with clinical isolates of B. rossi at two dose inocula and recorded clinical and laboratory measures of disease severity and progression on a daily basis. We extracted RNA from whole blood collected at six timepoints and subjected it to RNA-sequencing, differential gene expression analysis, pathway enrichment analysis, and temporal cluster analysis to identify genes and pathways that were differentially expressed as well as groups of genes that had similar expression trajectories over time. These results characterize the host response during the onset and treatment of severe babesiosis in canines and identify both known and novel pathways in the host response to severe hemoprotozoal infection.


Clinical time course correlated with infectious inoculum

All 5 experimentally infected subjects developed clinical disease after intravenous inoculation with 104 or 108 parasites. The first signs of overt clinical disease (fever, reduced habitus, and appetite loss) appeared on day 3 in the high inoculum cohort and on day 4 in the low inoculum cohort (Fig. 1). All subjects were treated on day 4 with diminazene aceturate (3.5 mg/kg intramuscularly), blood transfusion, and supportive care. One subject in the high inoculum cohort died suddenly on day 4 despite all attempts to treat it. The remaining subjects all made a full recovery and were homed as pets after completion of the experiment.

Fig. 1
figure 1

Clinical parameters recorded throughout the experiment reveal inoculum-dependent disease course and severity. Data are provided for each individual (dots) and for high and low inoculum cohort means (lines). ALT: alanine aminotransferase; BPM: beats per minute

Hematological parameters indicated severe hemolytic anemia

Parasitemia was detectable on day 2 in both the high and low inoculum cohorts. Hemoglobin concentration began to fall after day 2 in the high inoculum cohort and after day 3 in the low inoculum cohort. Plasma hemoglobin, an indicator of intravascular hemolysis, was elevated by day 3 in the high inoculum cohort. The absolute reticulocyte count was low for the degree of anemia present and peaked during peak infection (day 3 in the high inoculum and day 4 in the low inoculum cohort), then began to rise again after treatment had been administered. Both the total white blood cell count and immature neutrophil (band form) count rose steeply after day 5 in the high inoculum cohort. Severe thrombocytopenia developed by day 3 in the high inoculum cohort and by day 4 in the low inoculum cohort (Fig. 1).

Biochemical parameters revealed impact on multiple organ systems

On days 3 and 4, signs of metabolic acidosis (low pH, low HCO3, rising lactate) with respiratory compensation (low pCO2) developed in the high inoculum cohort. Serum urea and serum creatinine peaked on days 3 and 4, respectively, indicating an acute but transient decline in glomerular filtration. Serum alanine aminotransferase (ALT) peaked on day 5, indicating acute hepatocellular injury. Blood glucose fell precipitously in the high inoculum cohort and remained low until day 4 when treatment was administered. Cortisol, a stress hormone, peaked on day 4 in the high inoculum cohort. Thyroid hormone (T4) fell in both cohorts and was well below normal in the high inoculum cohort on days 4 and 5. C-reactive protein, a marker of inflammation, peaked on day 3 in the high inoculum cohort and on day 4 in the low inoculum cohort. Together, these biochemical findings describe the acute onset of a severe hemolytic and inflammatory disease affecting multiple organ systems.

Cytokine concentrations increased during infection and recovery

GM-CSF (granulocyte monocyte chemotactic factor), IL-6, macrophage chemotactic protein, and TNF (tumor necrosis factor) remained low during the first 4 days following inoculation, and then rose steeply after day 4 in one subject in the high inoculum cohort. Keratinocyte chemotactic-like cytokine and IL-10 remained low during the first 3 days following inoculation, peaked in concentration on day 4, then returned to low levels following treatment. IL-8 followed this pattern in the high inoculum cohort but not the low (Fig. 2).

Fig. 2
figure 2

Plasma cytokine concentrations increased with parasitemia or after treatment. Data are provided for each individual (dots) and for high and low inoculum cohort means (lines)

Quality control assessment reported high quality sequencing data

To understand the effect of B. rossi inoculation on the host transcriptomic response over time, blood was drawn into Tempus tubes on experimental days 0, 1, 3, 4, 6, and 8 from each subject. High-quality RNA (mean RINe 9.4 ± 0.29) was extracted from each sample, reverse transcribed, and subjected to 100 bp paired-end next generation sequencing on a NovaSeq instrument. The number and percentage of reads mapped to the canine genome was inversely related to parasitemia. On day 3 in the high inoculum cohort and on day 4 in both cohorts (when parasitemias were greatest), the number of mapped reads (range 4.0–118.6 million) and percentage of aligned reads (mean 10.3%) was lower than on all other days, when the number of mapped reads per sample ranged from 133.6 to 260.7 million and on average 82.5% aligned to the canine reference genome (Table S1). For all samples, the range of GC content was 47–51%. The sequence data are available in the NCBI database system (

Global expression profiles correlated with parasite density

To be included in the analysis, a gene transcript had to be present in two or more subjects in the same inoculum on the same day at 1 count per million (CPM) or greater. Eight thousand one genes passed this threshold and were used in subsequent analyses. Pairwise correlation of natural log-transformed counts for filtered host genes was high among a cluster composed of all day 0 samples and day 1 samples from the low inoculum cohort (mean r = 0.985; Fig. 3a). This implies a high degree of dependence of gene expression on parasitemia and a relatively low level of gene expression variation between subjects. A second cluster composed of samples from day 1 in the high inoculum cohort and day 3 in the low inoculum cohort (timepoints that coincided with the onset of parasitemia in each cohort) was also highly correlated (mean r = 0.984).

Fig. 3
figure 3

Correlation and cluster analyses of gene expression reveal primary sources of variation. A Heatmap showing pairwise correlation of gene expression across all samples. Raw counts were transformed using the nautral log. B Principal component analysis (PCA) of all samples. PCs were determined based on the natural log-transformed raw counts. Number label indicates day post-inoculation

Principal component analysis (PCA) of natural log-transformed counts of filtered host genes revealed tight clustering of all day 0 samples along with day 1 samples from the low inoculum cohort (Fig. 3b). The first principal component, which explained 92.37% of expression variation, distinctly separated samples by inoculum during post-treatment days 6 and 8. The second principal component, which explained 2.84% of expression variation, separated samples by day of infection. Subjects that received the low inoculum underwent changes in transcriptomic profiles that were distinct for days 3, 4, 6, and 8 after infection, with day 8 returning close to day 0. Subjects that received the high inoculum underwent similar changes in transcriptomic profiles, but the changes were evident at earlier timepoints. High inoculum samples separated on PC1 on days 6 and 8, with both samples on the far left corresponding to the subject that recorded the highest levels of parasitemia. In congruency with the observed progression of clinical signs and symptoms prior to treatment, the changes in transcriptional profiles among subjects that received the low inoculum appear to lag 1–2 days behind those that received the high inoculum as shown on PC2.

Differentially expressed genes (DEGs) were significantly upregulated throughout infection and recovery

Eight thousand one filtered genes were tested for differential expression (log2 fold change > |1|, FDR < 0.05), comparing expression on each experimental day versus day 0. In the low inoculum cohort, 1614 unique genes were differentially expressed on at least one experimental day (Fig. 4a). Twenty-six genes were differentially expressed on every experimental day after day 1 (on which only 4 genes were differentially expressed). In the high inoculum cohort, 2455 unique genes were differentially expressed on at least one experimental day (Fig. 4b). Fifty-two genes were differentially expressed on every experimental day after inoculation. Table 1 reports the numbers of DEGs that were up- and down-regulated on each day in each inoculum.

Fig. 4
figure 4

Pairwise differential gene expression for each experimental day versus day 0 in the low (A) and high (B) inoculum cohorts. Upregulated genes (log2(FC) > 1) are displayed in red; downregulated genes (log2(FC) < 1) are displayed in blue. Genes that did not meet the thresholds for significance (|log2(FC)| > 1; FDR < 0.05) are displayed in gray. Vertical dashed lines indicate the log2(FC) threshold; horizontal dashed line indicates the FDR threshold. Labeled DEGs comprise the top 20 with the highest absolute log2(FC) on each day in each cohort

Table 1 Up- and down-regulation of DEGs in the low and high inoculum cohorts throughout the experimental course

Twelve genes were differentially expressed on every experimental day after inoculation in both inoculum cohorts: RSAD2, PSTPIP2, GCLM, OBSCN, OAS1, ISG15, DDR1, IL7R, CCR7, CHST2, APOL2, and TCF7 (Table S2).

Gene ontology (GO)- and Reactome-annotated functional enrichment of DEGs

DEGs from each day in each of the high and low inoculum cohorts were separately subjected to pathway enrichment analyses using the GO and Reactome databases (Fig. S1). In the low inoculum cohort, GO pathway analysis identified 23 pathways on day 1, 221 pathways on day 3, 129 pathways on day 4, 208 pathways on day 6, and 243 pathways on day 8 to be significantly enriched. Reactome analysis identified 20 pathways on day 1, 29 pathways on day 3, 80 pathways on day 4, 68 pathways on day 6, and 56 pathways on day 8 as significantly enriched. In the high inoculum cohort, GO pathway analysis identified 281 pathways on day 1, 210 pathways on day 3, 282 pathways on day 4, 217 pathways on day 6, and 181 pathways on day 8 as significantly enriched. Reactome analysis identified 37 pathways on day 1, 48 pathways on day 3, 92 pathways on day 4, 115 pathways on day 6, and 95 pathways on day 8 as significantly enriched.

Cluster analysis revealed temporal expression trajectories of DEGs

Separate hierarchical cluster analyses were performed on the 1614 DEGs in the low inoculum and the 2455 DEGs in the high inoculum based on the similarity of each gene’s expression trajectory over time. Assessment of clusters using silhouette analysis [18] identified three robust clusters in each inoculum cohort, with silhouette widths (si) of 0.60, 0.57, and 0.50 (weighted mean si = 0.57) and si of 0.28, 0.49, and 0.51 (weighted mean si = 0.40) in the high inoculum (Fig. 5a, d).

Fig. 5
figure 5

Hierarchical clustering on the log2(FC) of DEGs reveals three distinct temporal trajectories in each the low (A, B, C) and high (D, E, F) inoculum cohorts. A, D Silhouette plots that demonstrate how well each gene fits in its respective cluster. Silhouette width (si) < 0: gene is better suited to a neighboring cluster; si > 0: gene is best suited in its assigned cluster. B, E Temporal trajectories of genes by cluster, shown as the average edgeR-calculated log2(FC) on each day compared to day 0 for each cluster. Cluster number is shown above each plot; color corresponds to silhouette plot color. Bars indicate 95% confidence intervals. C, F Constituent genes of each cluster that were differentially expressed on every experimental day in the high inoculum and every day except day 1 in the low inoculum cohort

Each cluster had a distinct temporal trend represented as the mean (95% CI) log2 fold change of all genes in a cluster over time (Fig. 5b,e). In both cohorts, cluster 1 consisted of genes that were downregulated as infection progressed and rose toward baseline after treatment. In cluster 2, genes were upregulated as infection progressed and downregulated after treatment. Cluster 3 trajectories differed slightly between the cohorts. In the low inoculum, genes were upregulated during infection onset (days 1 and 3), downregulated slightly during peak infection (day 4), and upregulated after treatment (day 6) (Fig. 5b). In the high inoculum, cluster 3 consisted of genes that were upregulated during infection onset (day 1) and continued to increase expression levels after treatment (day 6) (Fig. 5e). Representative gene members of each cluster were identified as those that were differentially expressed at all time points in the high inoculum group and at all time points from day 3 onwards in the low inoculum group; these are listed in Fig. 5c and f.


Our understanding of the pathophysiology of canine babesiosis has, until now, been based on clinical observational studies of natural infection in a broad range of canine hosts who had been exposed to different infectious inocula at variable times in the past. To constrain these sources of clinical heterogeneity, we performed controlled inoculations at two doses in susceptible, purpose-bred laboratory beagles and characterized the infection in terms of clinical manifestations and gene expression changes over time. In doing so, we learned that experimental inoculation with wild B. rossi parasites consistently induces the rapid onset of symptomatic canine babesiosis. Both the rate of onset and the severity of disease were related to the size of the infectious inoculum. Disease severity, reflected through behavior, vital signs, hematology, biochemistry, acid-base status, and hormone levels, increased rapidly with the rising parasitemia and responded to treatment.

To further understand the progression of this hemoprotozoal infection which transpires entirely within the hematological system, we monitored host gene expression in circulating whole blood at frequent intervals. The results of transcriptomic analysis identified differential expression of pathways involving (1) response to hemolysis (including erythropoiesis and iron homeostasis), (2) metabolism, and (3) distinct immunological responses including innate immunity, adaptive immunity, and viral response genes. These highlight host pathways that could be targeted to slow disease progression and prevent life-threatening complications. Understanding the pathogenesis of severe babesiosis in canines can not only help us to improve the treatment of canine babesiosis, but also provide insights into the pathophysiology of human babesiosis and the related parasitic disease, severe malaria.


Anemia is a hallmark of babesiosis in humans, horses, cattle, and canines [7]. In this experimental model, massive hemolysis, as indicated by the rapid rise in plasma hemoglobin, along with inadequate erythropoiesis caused rapid and severe anemia that required transfusion to correct. After an early rise in reticulocyte counts, the reticulocytosis was suppressed as parasitemia climbed and remained inadequate for several days after treatment with an anti-parasite drug, implying that the parasite or the host response to infection was limiting erythropoiesis or blocking the release of reticulocytes from the bone marrow into circulation.

At a transcriptional level, enrichment of erythrocyte-related genes and pathways was prominent during peak parasitemia and also during recovery. During peak infection in both inoculum cohorts, oxygen transport (GO:0015671), iron ion homeostasis (GO:0055072), erythrocytes take up oxygen and release carbon dioxide (R-CFA-1247673), heme biosynthesis (R-CFA-189451), and metabolism of porphyrins (R-CFA-189445) were all among the top 10 most enriched pathways as measured by negative log FDR (Fig. S1).

Many genes with functions related to erythopoesis and heme biosynthesis were in cluster 2, which increased in expression as infection progressed and returned to baseline during recovery. This includes KLF1, a key transcription factor regulating erythropoiesis and hemoglobin expression [19, 20], ALAS2, which catalyzes the first and rate limiting step in heme synthesis in erythroid cells [21], and UROD, which codes an enzyme in the heme biosynthetic pathway [22]. In the same cluster, SPTB and SPTA1 were differentially expressed; these genes are expressed in erythroid precursors and code for a constituent of the erythrocyte cytoskeleton [23]. Mutations in these genes are also associated with hereditary spherocytosis [24]. The rise and fall in erythroblast transcripts parallels the reticulocyte count, which peaked during peak infection (day 3 in the high inoculum and day 4 in the low inoculum cohort) and began to rise again after treatment. This association is reflected quantitatively by high correlation values between erythroblast transcript expression trajectories and the reticulocyte count (Table S3), and suggests that the changes in expression of these genes may be due to reticulocyte population dynamics.

Other genes related to hemolysis were in cluster 2. This includes HP, the gene that codes for haptoglobin, a protein that binds free hemoglobin from lysed red cells and is a marker for hemolysis [25]. Additionally PRDX2, which encodes an antioxidant enzyme that stabilizes hemoglobin under oxidative conditions, was differentially expressed in both cohorts. Mutations in this gene are strongly associated with hemolytic anemia [26]. Finally, GCLM, which has been shown to be essential to erythrocyte survival during oxidative stress [27] and whose deficiency is associated with hemolytic anemia [28], was also upregulated during infection.

Though there have not been extensive studies correlating gene expression with hemolysis in clinical malaria pathogenesis [29], the transcriptomic profiles reported here share many features in common with the host response to hemolytic anemia, including in sickle cell disease [20, 30,31,32,33].


Both high and low dose inoculation with B. rossi was associated with enrichment of the Reactome pathway metabolism (R-CFA-1430728) during peak infection. In the low inoculum, long-chain fatty acid metabolic process (GO:0001676) was enriched on days 4 and 6 and cholesterol metabolic process (GO:0008203) was enriched on days 3 and 4. Both pathways were also enriched in the high inoculum on days 3 and 4 (Fig. S1).

Many lipid metabolism-related genes were in cluster 2 in both inoculum cohorts. This includes: PLA2G4A, which encodes a phospholipase involved in membrane lipid remodeling and biosynthesis of lipid mediators of the inflammatory response [34, 35]; GPR84, which serves as a receptor for free fatty acid [36]; CD5L, which encodes a regulator of lipid synthesis that in turn regulates inflammatory response mechanisms and T cell activities [37]; and ALOX15, which encodes a lipoxygenase that catalyzes the deoxygenation of polyunsaturated fatty acids and has known roles in red cell maturation and the inflammatory immune response [38]. Several lipid metabolism-related genes that were in cluster 2 in the low inoculum cohort were in cluster 3 in the high, indicating that they continued to increase expression levels throughout infection and recovery. This includes both ECHDC3, known to play a major role in fatty acid biosynthesis [39], and FADS1, which encodes a fatty acid desaturase [40]. Previous studies implicated the role of lipid droplets as regulators of the immune response to protozoan infection [41]. The high prevalence of lipid metabolism genes among differentially expressed genes implicates these pathways in the host response to B. rossi infection (Fig. 5).

Genes related to glucose metabolism were also differentially expressed throughout the infection and recovery. In the high inoculum cohort, glucose metabolic process (GO:0006006) was enriched on days 6 and 8. GATM, which encodes the rate-limiting step in creatine metabolism [42], was in cluster 3 in the low inoculum cohort and cluster 2 in the high. TKT, a key enzyme in the pentose phosphate pathway [43], was in cluster 2 in both cohorts, indicating an increase in glycolytic activity. Hypoglycemia and hyperlactatemia are a feature of both babesia and malaria infections [44, 45]. Studies in malaria show glycolysis pathway genes prominently activated early during infection in mouse models and in humans [29, 46]. Thus, canine babesiosis and human malaria appear to share features of both lipid and carbohydrate metabolic pathway activity.

Immune response

Both the clinical and transcriptomic data indicated a robust immune response to infection in the high and low inoculum cohorts. The inflammatory nature of the disease was reflected by the inoculum-dependent changes in the fever curve, CRP level, as well as the white cell, neutrophil, and band counts – the last of which is known to predict poor outcomes [12]. While some cytokines reached a maximum at peak parasitemia (IL-8, IL-10, KC-like) others rose only after the infection was treated (TNF, IL-6, MCP-1). Multiple immune pathways were differentially regulated in response to both high and low inocula. Defense response (GO:0006952) was enriched in the high inoculum on day 1 and both inoculum cohorts on day 3. Immune system/response (R-CFA-168256; GO:0006955) was enriched in the high inoculum on days 1 and 8, and both inoculum cohorts on days 3–6. Innate immune system/response (R-CFA-168249; GO:0006955) was enriched on every experimental day in the high inoculum cohort. Additionally, cytokine signaling in immune system (R-CFA-1280215) was enriched during infection onset in both inoculum cohorts (day 1 in the high and day 3 in the low inoculum) These pathways are composed of hundreds of constituent genes, which primarily followed cluster 2 trajectories. These pathways reflect the differential expression of specific genes involved in innate immunity, adaptive immunity, and viral response (Figs. 4, S1).

Innate immunity

Genes associated with positive regulation of interleukin-1 beta (IL-1β) production (GO:0032731) were in cluster 2. This includes: IL1B, which encodes the pro-inflammatory cytokine [47]; PSTPIP2, which encodes a protein that regulates IL-1β [48, 49]; CASP4, which regulates IL-1β synthesis in macrophages [50]; and SIGLEC14, which encodes a receptor that enhances inflammasome activation and macrophage IL-1β release. Furthermore, genes associated with NF-kB signaling were also constituents cluster 2, including FABP5, ACOD1, and TNFAIP2. FABP5 encodes a fatty acid binding-protein that promotes the activation of NF-kB signaling [51]. ACOD1 serves a multifaceted role in the innate immune response to pathogen infection (bacteria and viruses) and in cytokine signaling (TNF and interferons), and is associated with both Toll-like receptor and NF-kB signaling pathways [52]. Expression of TNFAIP2 is upregulated via the NF-kB pathway and mediates the inflammatory response [53, 54]. Studies of malaria infection in humans also implicate the role of IL-1β signaling and NF-kB pathways in the innate immune response to infection [55, 56]. Thus, Babesia and Plasmodium induce similar innate immune responses across different host species.

Pathways associated with the production and regulation of interferon-gamma (GO:0032649) were significantly enriched in the host response to B. rossi infection, especially during infection onset (Fig. S1). IFN-γ signaling-related genes, including IRGM (IFI1), IFGGC1, SPP1, and APOL2, were in cluster 2 in both inoculum cohorts. IRGM plays a key role in regulating IFN-γ-induced host defense to protozoa [57], while IFGGC1 is a member of the IFN-γ-inducible GTPases that are involved in immune response against pathogens [58]. SPP1 codes for a cytokine known to upregulate IFN-γ production [59], and APOL2 which mitigates the cytotoxic effects of IFN-γ [60]. IFN-γ is a major pro-inflammatory cytokine that has both direct antiparasitic activity and immunoregulatory roles in the response to parasitic infection [61]. IFN-γ has a role in protective immunity against infection by B. microti in mice [62] and is involved in the host immune response to malaria in both mice and humans [29, 55, 56]. Our findings suggest IFN-γ has a similar role in protective immunity against B. rossi infection in canines.

Adaptive immunity

Genes that code the MHC class II protein complex, namely DLA-DOA and DLA-DOB [63], were in cluster 1, indicating that they were significantly downregulated in both inoculum cohorts. In contrast, MHC class I genes DLA-12, DLA-64, and DLA-88 were upregulated during infection (though they did not meet the log2 fold change threshold for significance, excepting the low inoculum cluster 2 gene DLA-12). Additionally, GZMA and PRDX2, genes specifically expressed by CD8+ T cells and associated with their activity [64, 65], were cluster 2 DEGs in both cohorts.

CD8+ and CD4+ T cells, whose T-cell receptors recognize antigens presented by MHC-I and MHC-II respectively, have been previously implicated in protective immunity against B. microti infection in mice and B. bovis infection in cattle. Previous studies in mice and humans suggest that both MHC-I and MHC-II play a role in the host immune response to malaria infection [29, 56]. However, recent studies in human cerebral malaria suggest that CD8+ T cells migrate to and sequester in the brain and contribute to mortality [66, 67]. Our data reveal downregulation of genes associated with MHC class II CD4+ T cells and potential induction of MHC class I CD8+ T cells in the host response to B. rossi, raising the question of what cells are presenting antigens and whether the CD8+ and CD4+ responses are immunopathologic or contribute to parasite killing and adaptive immunity.

Viral response

Interestingly, many viral response genes were differentially regulated during B. rossi infection, particularly during infection onset. On day 1 in the high inoculum cohort, pathways including positive regulation of defense response to virus by host (GO:0002230), antiviral mechanisms by IFN-stimulated genes (R-CFA-1169410), and ISG15 antiviral mechanism (R-CFA-1169408) were enriched. The latter two pathways were also enriched on day 3 in the low inoculum cohort, along with response to interferon-alpha (GO:0035455) and cellular response to interferon-beta (GO:0035458). (Fig. S1).

This is also reflected in the significant upregulation of constituent viral response genes throughout the infection in both inoculum groups. All of the following genes were identified as cluster 2 DEGs: RSAD2, ISG15, DDX58, DDX60, OAS1, OAS2, OASL1, MX1, PLSCR1, ICAM1, TRIM10, KLRD1, LGALS9, and LTF. In general these genes are related to type I interferon (IFN-I) antiviral signalling and have been associated with a range of viruses including influenza, west Nile virus, hepatitis C, and dengue virus [68,69,70,71,72,73,74,75,76,77,78,79,80,81,82]. Interestingly, previous studies in malaria and leishmaniasis models suggest that type I interferons may also suppress anti-parasitic immunity [83,84,85,86]. CD8+ T cells are also associated with antiviral activity [64, 65, 87, 88]. Together, these transcriptomic data implicate the role of viral response genes and pathways in the host immune response to B. rossi infection.

In the context of Plasmodium infection, parasite DNA and RNA have been shown to bind to receptors that activate downstream production of IFN-Is. For one, toll-like receptor 9 (TLR9) recognizes parasite DNA [89, 90]. Cyclic GMP-AMP synthase (cGAS) also binds to parasite DNA in the cytoplasm and catalyzes the synthesis of cyclic guanosine monophosphate-adenosine monophosphate (cGAMP), which acts as a second messenger to activate the downstream signal adaptor, stimulator of interferon genes (STING) [89, 91, 92]. Furthermore, major cytosolic receptor for sensing malaria parasite RNA is melanoma differentiation-associated protein 5 (MDA5) [92, 93]. Parasite RNA also possibly stimulates the retinoic acid-inducible gene I (RIG-I) pathway [94], though this has yet to be confirmed in vivo [89]. It is likely that Babesia parasite DNA and RNA similarly activate receptors that activate downstream production of IFN-Is. Of note, it is possible that different sensors are activated by parasite DNA and RNA at different developmental stages, and IFN-I subtypes (IFN-α and IFN-β) have different biological and immunological roles in viral infections [89]. Future experiments will explore the mechanistic roles of IFN-Is in B. rossi infection.


The inoculum of 109 parasites resulted in an acute and severe disease course. This rapid disease course may be more representative of transfusion-associated babesiosis than natural tick-transmitted infections which would evolve more slowly. Future experiments will utilize a lower infectious inoculum, which may be more representative of natural infections. B. bovis and Plasmodium falciparum are known to cytoadhere to the vascular endothelium in severe babesiosis and malaria [95,96,97,98], respectively, while P. vivax and P. knowlesi can cause severe disease without evidence of cytoadherence [99,100,101,102]. The available evidence suggests that B. rossi causes severe disease without cytoadherence to vascular endothelium [103,104,105]. The low-mapped reads on day 4 in the high and low inoculum cohorts and day 3 in the high appear to be associated with the high levels of parasitemia recorded on those days. The lower percentage of mapped reads indicates lower overall sequencing accuracy and the potential presence of contaminating DNA, thus increasing the likelihood of spurious results on these days. Additionally, one animal in the high inoculum died on the morning of day 4, which may have affected the expression means reported on days 4, 6, and 8. Transcriptomic and clinical results reported on days 6 and 8 may have been further affected by the administration of an anti-parasite drug diminazene and blood transfusions on day 4. While the primary effect of the drug is to reduce parasitemia, it may also have a weak anti-inflammatory effect [106]. Furthermore, both cohorts were small (low n = 2; high n = 3); however, there was a high degree of consistency between subjects of the same cohort with regard to clinical observations and transcriptomic data.

It is likely that cellular population dynamics account for many of the changes in gene expression observed, as exemplified by the correlation of reticulocyte count with erythroblast transcripts. When pure cellular reference populations are available, cell-type deconvolution can be applied to understand the relative contributions of changes in gene expression versus changes in cell abundance. Another approach would be to measure gene expression at the single cell level. Either of these techniques could be employed to identify transcripts that changed in expression levels within a cell population as opposed to transcripts whose changes in expression were associated with changing cellular populations.


This study is the first to describe canine host gene expression induced by B. rossi infection. The rate of onset and disease severity were inoculum-dependent, as shown by clinical data including vital signs, hematological parameters, parasitemia, acid-base chemistry, and biochemical markers. The temporal dynamics of differential gene expression and pathway enrichment implicate response to parasite burden, hemolysis, induction of lipid metabolism, and a robust role of different immunological pathways including those related to innate immunity, adaptive immunity, and response to viral infection. These findings offer a fundamental understanding of the host transcriptomic response to B. rossi, which can serve as a model of human hemolytic diseases caused by intracellular parasites, such as babesiosis and malaria.



All methods were carried out in accordance with relevant guidelines and regulations. Six purpose bred male beagles were born and raised in the Onderstepoort Veterinary Animal Research Unit. The animals were all microchipped, vaccinated and dewormed according to international guidelines and neutered at 6 months of age. Stringent ectoparasite control was applied by means of the isolation of the housing facility and 3 monthly dosing with an isoxazoline compound. All dogs were proven to be free of any blood borne parasites by PCR and reverse line blot. All dogs were housed communally in a play enriched environment and socialized with animal care givers each day. They were trained through a food reward system to become accustomed to being held for blood collection and blood pressure monitoring on a daily basis for several months before the beginning of the experiment.

Creation of a Babesia rossi cryopreservate

Three to five milliliters of whole venous blood were collected into EDTA vacutainer tubes from dogs naturally infected with Babesia rossi and presented for veterinary care at the Onderstepoort Veterinary Academic Hospital and cryopreserved in dimethylsuphoxide. All samples collected this way were subjected to PCR and reverse line blot to confirm a mono-infection with B. rossi. One of these samples was later used to infect a splenectomized Beagle dog to raise the infectious inoculum required for the infection of the experimental dogs.

Experimental infection

Animals were cared for in compliance with the ARRIVE (Animal Research: Reporting of In Vivo Experiments) guidelines on an experimental protocol approved by the University of Pretoria Animal Ethics Committee (V003–18). One of the 6 Beagles was splenectomized at 7 months of age, allowed to recover for 3 weeks and then infected intravenously with one of the wild-type B. rossi cryopreserved samples by intravenous injection. This canine was monitored closely, and once a patent parasitemia of ~ 2% was detected on blood smear of central venous blood, the parasitaemia was accurately quantified, and 10 mL of blood was collected into acid citrate dextrose solution and diluted with a culture medium used for in-vitro babesia culture to create infectious inocula, two doses of 104 parasites and three of 108 parasites. The splenectomized dog was immediately treated and made a full recovery without becoming obviously ill. These inocula were immediately injected intravenously into two canines that received the low infectious dose and three that received the high infectious dose. These five canines were then followed carefully through frequent physical observation and daily blood collections.

Biospecimen collection

Blood was collected daily into EDTA tubes for complete blood counts and the harvesting of plasma for cytokine determination, clot tubes for serum biochemistry, heparinized blood collected anaerobically for venous blood gas, lactate and glucose determination, citrated tubes for coagulation studies and RNALater tubes (Tempus), for the transcriptomic study reported here. Baseline data was generated twice at a 2-week interval before any infections took place to generate a baseline.

Clinical and laboratory data collection

Clinical examinations and all blood samples were collected from the jugular vein with 21G vacutainer needles (Precision Glide™, UK) between 8:00 AM and 10:00 AM daily. The following observational data were collected once every day for each experimental subject: habitus and mental status (scored 1–4: 1-collapsed, 2-very weak, just able to stand, 3-able to stand but depressed, 4-normal interaction); appetite (scored 1–4: 1-no interest in food, 2- slight interest with minimal ingestion, 3-partial ingestion of a full meal, 4-normal ingestion of the whole meal); rectal temperature; pulse (palpated on the femoral artery); respiratory rate counted by physical observation; mucous membrane color and capillary refill time; thoracic auscultation; abdominal palpation; stained blood smear examination; blood pressure.

Haematology was performed on a venous sample collected into an EDTA Vacutainer Brand Tubes (Beckton Dickinson Vacutainer Systems, UK) from the jugular and was run on an ADVIA 2120 (Siemens, Munich, Germany). within an hour of collection A differential count was performed manually by an experienced laboratory technologist. Siemens, Germany) within 1 h of blood collection. The remaining volume of blood was then centrifuged and aliquoted within 30 min of collection for storage of EDTA plasma at -80 °C.

Serum biochemistry samples were collected in Serum Vacutainer Brand Tubes (Beckton Dickinson Vacutainer Systems, UK), and the CRP (using canine specific immunoturbidimetric CRP method), albumin (using a colorimetric assay with bromocresal green) were determined. Glucose (using the point of care AlphaTRAK 2) and the automated hexokinase method was run on an automated chemistry analyser (Cobas Integra© 400Plus). This sample was allowed to clot for 10 min, samples were run and the remaining serum volumes were aliquoted within 30 min of collection for storage at -80 °C.

Blood gas samples were collected anaerobically into a commercially prepared heparinized syringe (BD A-Line, arterial blood collection syringe, Becton, Dickinson and Company, UK) using a 21G needle from the jugular. Lactate was obtained directly from a fresh blood sample using a Lactate Pro 2 hand held lactate reader and from the venous blood gas analysis, analyzed within 20 min (Rapidpoint 405, Seimens).

The cytokines measured included: GM-CSF, IFN-γ, IL-2, IL-6, IL-7, IL-8, IL-15, IL-10, IL-18, TNF-α, IP 10, KC and MCP-1. Concentrations were determined in duplicate by fluorescent-coded magnetic beads (MagPlex-C; MILLIPLEX. MAP Kit, Canine Cytokine Magnetic Bead Panel, 96-Well Plate Assay, CCYTO-90 K, Millipore, Billerica, MA), based on the Luminex xMAP technology (Luminex 200, Luminex Corporation, Austin, TX). Two quality controls were included in the plate as internal quality controls. The assay was performed according to the manufacturer’s instructions.

Blood pressure was measured using the Vet HDO® MDPro and the HDO management software. The following protocol for blood pressure measurement was standardized: The environment was isolated, quiet and away from other animals. The same environment was used every day and the canines were trained and conditioned for 4 weeks prior to the onset of the experiment. No sedation was used, and the canines were allowed 5–10 min to relax and become accustomed to the environment. Each subject was gently restrained in right lateral recumbency. The cuff width was approximately 40% of the circumference of the cuff site and the cuff size was recorded at every measurement. The cuff was placed around the tail base. All blood pressure measurements were performed by the same person every day. The canines were calm and relatively motionless as they were all well-conditioned to the process prior to the start of the experimental period. The first measurement was discarded. Five consecutive and consistent (< 20% variability in systolic values) were recorded and the average of the values was calculated to obtain the blood pressure measurement.

RNA isolation

Blood from each subject was collected into Tempus Blood RNA Tubes on days 0, 1, 3, 4, 6, and 8 and subsequently stored at -20 °C. Samples were shipped on dry ice to the United States, where they remained at -20 °C prior to RNA extraction. RNA was isolated using Qiagen RNeasy Mini Kit. RNA quality was controlled on a Bioanalyzer Instrument. All RNA samples had an RNA Integrity Number (RIN) of ≥9.7. Two independent biological replicates were selected for each time point for subsequent RNA sequencing.

RNAseq library preparation and sequencing

PolyA-selected mRNA libraries were generated and all samples were indexed and pooled before sequencing in a single S1 flowcell on an Illumina NovaSeq.

Data analysis

Reads were trimmed and mapped to the canine reference genome using the Center for Cancer Research Collaborative Bioinformatics Resource (CCBR) Pipeliner. Pipeliner provides open access to some of the next-generation sequencing data analysis pipelines used by CCBR. It allows users to select a set of genomes and data files, then processes them through a sequence of programs to produce the desired output (i.e. raw counts matrix, QC report, differentially expressed genes, etc.). Specifically, Pipeliner employs the reference-mapping tool STAR, which is a widely used and accepted tool and allows for mismatches and sequence divergence well beyond that expected amongst canine breeds [107]. Once configured, the pipeline is executed on the NIH high performance computing cluster Biowulf. Pipeliner is publicly available on GitHub [].

Raw gene counts were converted to counts per million (CPM). Genes with fewer than 1 CPM in two or more treatment group samples (inoculum + day) were excluded from analysis. Differential expression analysis was run for each experimental day compared to day 0 in each inoculum cohort using the edgeR package (3.30.3), and p-values were corrected using the Benjamini-Hochberg (BH) method. Genes with FDR < 0.05 and log2 fold change > |1| were considered significantly differentially expressed. Gene Ontology (GO)-annotated pathway enrichment analysis was performed using the topGO package (2.40.0); package authors discourage multiple testing correction as the p-values returned are already conservative [108]. Reactome-annotated pathway enrichment analysis was performed using the biomaRt (2.44.0) and reactome.db (1.70.0) packages, with significance assigned according to a hypergeometric distribution and p-values corrected using the BH method. Hierarchical cluster analyses were performed separately for each inoculum cohort, and included all genes that were differentially expressed on at least one experimental day. The Pearson correlation between genes was calculated using the edgeR-calculated log2 fold change across experimental days. The Euclidean distance between those correlation values was then calculated and clustering was performed using Ward’s method. Cluster robustness was assessed using silhouette analysis as part of the factoextra R package (1.0.7). To visualize cluster temporal trajectories, the mean edgeR-calculated log2 fold change on each day within each cluster was plotted along with 95% confidence intervals. Data analysis was conducted in R version 4.0.0.

Availability of data and materials

The datasets supporting the conclusions of this article are available on the NCBI Gene Expression Omnibus (GEO) database as Series GSEI67201 []. The code generated for this analysis is available at the GitHub repository Babesia_RNAseq [].



Alanine aminotransferase


Counts per million


Differentially expressed genes


Fold change


False discovery rate


Granulocyte monocyte chemotactic factor


Gene Ontology


National Center for Biotechnology Information


Principal Components Analysis


RNA integrity number


Tumor necrosis factor


  1. Elahi R, Islam A, Hossain MS, Mohiuddin K, Mikolon A, Paul SK, et al. Prevalence and diversity of avian haematozoan parasites in wetlands of Bangladesh. J Parasitol Res. 2014;2014:493754.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Abubakar II, Tillmann T, Banerjee A. Global, regional, and national age–sex specific all-cause and cause-specific mortality for 240 causes of death, 1990–2013: a systematic analysis for the Global Burden of Disease Study 2013. Lancet. 385(9963):117–71.

  3. Gorenflot A, Moubri K, Precigout E, Carcy B, Schetters TP. Human babesiosis. Ann Trop Med Parasitol. 1998;92(4):489–501.

    Article  CAS  PubMed  Google Scholar 

  4. Kjemtrup AM, Conrad PA. Human babesiosis: an emerging tick-borne disease. Int J Parasitol. 2000;30(12–13):1323–37.

    Article  CAS  PubMed  Google Scholar 

  5. Maharana BR, Tewari AK, Saravanan BC, Sudhakar NR. Important hemoprotozoan diseases of livestock: challenges in current diagnostics and therapeutics: an update. Vet World. 2016;9(5):487–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Schoeman JP. Canine babesiosis. Onderstepoort J Vet Res. 2009;76(1):59–66.

  7. Homer MJ, Aguilar-Delfin I, Telford SR, Krause PJ, Persing DH. Babesiosis. Clin Microbiol Rev. 2000;13(3):451–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Levin AE, Krause PJ. Transfusion-transmitted babesiosis: is it time to screen the blood supply? Curr Opin Hematol. 2016;23(6):573–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Schnittger L, Rodriguez AE, Florin-Christensen M, Morrison DA. Babesia: a world emerging. Infect Genet Evol. 2012;12(8):1788–809.

    Article  PubMed  Google Scholar 

  10. Penzhorn BL. Why is southern African canine babesiosis so virulent? An evolutionary perspective. Parasit Vectors. 2011;4(1):51.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Collett MG. Survey of canine babesiosis in South Africa. J S Afr Vet Assoc. 2000;71(3):180–6.

    Article  CAS  PubMed  Google Scholar 

  12. Leisewitz AL, Goddard A, Clift S, Thompson PN, de Gier J, Van Engelshoven JMAJAJ, et al. A clinical and pathological description of 320 cases of naturally acquired Babesia rossi infection in dogs. Vet Parasitol. 2019;271:22–30.

    Article  PubMed  Google Scholar 

  13. Eichenberger RM, Riond B, Willi B, Hofmann-Lehmann R, Deplazes P. Prognostic markers in acute Babesia canis infections. J Vet Intern Med. 2016;30(1):174–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Leisewitz A, Goddard A, de Gier JD, Engelshoven JV, Clift S, Thompson P, et al. Disease severity and blood cytokine concentrations in dogs with natural Babesia rossi infection. Parasite Immunol. 2019;41(7):e12630.

    Article  CAS  PubMed  Google Scholar 

  15. Henning A, Clift SJ, Leisewitz AL. The pathology of the spleen in lethal canine babesiosis caused by Babesia rossi. Parasite Immunol. 2020;42(5):e12706.

    Article  PubMed  Google Scholar 

  16. Krause PJ, Daily J, Telford SR, Vannier E, Lantos P, Spielman A. Shared features in the pathobiology of babesiosis and malaria. Trends Parasitol. 2007;23(12):605–10.

    Article  CAS  PubMed  Google Scholar 

  17. Maegraith B, Gilles HM, Devakul K. Pathological processes in Babesia canis infections. Z Tropenmed Parasitol. 1957;8(4):485–514.

    CAS  PubMed  Google Scholar 

  18. Rousseeuw PJ. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J Comput Appl Math. 1987;20:53–65.

    Article  Google Scholar 

  19. Huang J, Zhang X, Liu D, Wei X, Shang X, Xiong F, et al. Compound heterozygosity for KLF1 mutations is associated with microcytic hypochromic anemia and increased fetal hemoglobin. Eur J Hum Genet. 2015;23(10):1341–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Perkins A, Xu X, Higgs DR, Patrinos GP, Arnaud L, Bieker JJ, et al. Krüppeling erythropoiesis: an unexpected broad spectrum of human red blood cell disorders due to KLF1 variants. Blood. 2016;127(15):1856–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Muckenthaler MU, Rivella S, Hentze MW, Galy B. A Red Carpet for Iron Metabolism. Cell. 2017;168(3):344–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Franken ACW, Lokman BC, Ram AFJ, Punt PJ, van den Hondel CAMJJ, de Weert S. Heme biosynthesis and its regulation: towards understanding and improvement of heme biosynthesis in filamentous fungi. Appl Microbiol Biotechnol. 2011;91(3):447–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Qin L, Nie Y, Zhang H, Chen L, Zhang D, Lin Y, et al. Identification of new mutations in patients with hereditary spherocytosis by next-generation sequencing. J Hum Genet. 2020;65(4):427–34.

    Article  CAS  PubMed  Google Scholar 

  24. Wang X, Liu A, Lu Y, Hu Q. Novel compound heterozygous mutations in the SPTA1 gene, causing hereditary spherocytosis in a neonate with coombs-negative hemolytic jaundice. Mol Med Rep. 2019;19(4):2801–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Shih AWY, McFarlane A, Verhovsek M. Haptoglobin testing in hemolysis: measurement and interpretation. Am J Hematol. 2014;89(4):443–7.

    Article  CAS  PubMed  Google Scholar 

  26. Emberesh M, Giger Seu K, Emberesh S, Trump L, Risinger M, Zhang W, et al. Peroxiredoxin II (PRDX2) Is a Novel Candidate Gene for Congenital Dyserythropoietic Anemia. Blood. 2018;132(Supplement 1):3605.

    Article  Google Scholar 

  27. Föller M, Harris IS, Elia A, John R, Lang F, Kavanagh TJ, et al. Functional significance of glutamate–cysteine ligase modifier for erythrocyte survival in vitro and in vivo. Cell Death Differ. 2013;20(10):1350–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Almusafri F, Elamin HE, Khalaf TE, Ali A, Ben-Omran T, El-Hattab AW. Clinical and molecular characterization of 6 children with glutamate-cysteine ligase deficiency causing hemolytic anemia. Blood Cells Mol Dis. 2017;65:73–7.

    Article  CAS  PubMed  Google Scholar 

  29. Lee HJ, Georgiadou A, Otto TD, Levin M, Coin LJ, Conway DJ, et al. Transcriptomic studies of malaria: a paradigm for investigation of systemic host-pathogen interactions. Microbiol Mol Biol Rev. 2018;82(2):e00071–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Griffin PJ, Sebastiani P, Edward H, Baldwin CT, Gladwin M, Gordeuk V, et al. The Genetics of Hemoglobin A2 Regulation in Sickle Cell Anemia. Am J Hematol. 2014;89(11):1019–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Hamda CB, Sangeda R, Mwita L, Meintjes A, Nkya S, Panji S, et al. A common molecular signature of patients with sickle cell disease revealed by microarray meta-analysis and a genome-wide association study. PLoS One. 2018;13(7):e0199461.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Peralta R, Low A, Kim A, Murray S, Guo S, Freier S, et al. Targeting BCL11A and KLF1 for the treatment of sickle cell disease and β-thalassemia in vitro using antisense oligonucleotides. Blood. 2013;122(21):1022.

    Article  Google Scholar 

  33. Saraf SL, Viner M, Rischall A, Raslan R, Shah BN, Zhang X, et al. HMOX1 and acute kidney injury in sickle cell anemia. Blood. 2018;132(15):1621–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Gubern A, Casas J, Barceló-Torns M, Barneda D, de la Rosa X, Masgrau R, et al. Group IVA phospholipase A2 is necessary for the biogenesis of lipid droplets. J Biol Chem. 2008;283(41):27369–82.

    Article  CAS  PubMed  Google Scholar 

  35. Paloschi MV, Lopes JA, Boeno CN, Silva MDS, Evangelista JR, Pontes AS, et al. Cytosolic phospholipase a 2 -α participates in lipid body formation and PGE 2 release in human neutrophils stimulated with an l -amino acid oxidase from Calloselasma rhodostoma venom. Sci Rep. 2020;10(1):10976.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Montgomery MK, Osborne B, Brandon AE, O’Reilly L, Fiveash CE, Brown SHJ, et al. Regulation of mitochondrial metabolism in murine skeletal muscle by the medium-chain fatty acid receptor Gpr84. FASEB J. 2019;33(11):12264–76.

    Article  CAS  PubMed  Google Scholar 

  37. Wang C, Yosef N, Gaublomme J, Wu C, Lee Y, Clish CB, et al. CD5L/AIM regulates lipid biosynthesis and restrains Th17 cell pathogenicity. Cell. 2015;163(6):1413–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Ivanov I, Kuhn H, Heydeck D. Structural and functional biology of arachidonic acid 15-lipoxygenase-1 (ALOX15). Gene. 2015;573(1):1–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Sharma NK, Key C-CC, Civelek M, Wabitsch M, Comeau ME, Langefeld CD, et al. Genetic regulation of Enoyl-CoA hydratase domain-containing 3 in adipose tissue determines insulin sensitivity in African Americans and Europeans. Diabetes. 2019;68(7):1508–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. He Z, Zhang R, Jiang F, Zhang H, Zhao A, Xu B, et al. FADS1-FADS2 genetic polymorphisms are associated with fatty acid metabolism through changes in DNA methylation and gene expression. Clin Epigenet. 2018;10(1):113.

    Article  Google Scholar 

  41. Vallochi AL, Teixeira L, Oliveira K, Maya-Monteiro CM, Bozza PT. Lipid Droplet, a Key Player in Host-Parasite Interactions. Front Immunol. 2018;9:1022.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Kazak L, Chouchani ET, Lu GZ, Jedrychowski MP, Bare CJ, Mina AI, et al. Genetic depletion of adipocyte creatine metabolism inhibits diet-induced thermogenesis and drives obesity. Cell Metab. 2017;26(4):660–671.e3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Qin Z, Xiang C, Zhong F, Liu Y, Dong Q, Li K, et al. Transketolase (TKT) activity and nuclear localization promote hepatocellular carcinoma in a metabolic and a non-metabolic manner. J Exp Clin Cancer Res. 2019;38(1):154.

    Article  PubMed  PubMed Central  Google Scholar 

  44. English M, Wale S, Binns G, Mwangi I, Sauerwein H, Marsh K. Hypoglycaemia on and after admission in Kenyan children with severe malaria. QJM Int J Med. 1998;91(3):191–7.

    Article  CAS  Google Scholar 

  45. Keller N, Jacobson LS, Nel M, de Clerq M, Thompson PN, Schoeman JP. Prevalence and risk factors of hypoglycemia in virulent canine babesiosis. J Vet Intern Med. 2004;18(3):265–70.

    Article  PubMed  Google Scholar 

  46. Sexton AC, Good RT, Hansen DS, D’Ombrain MC, Buckingham L, Simpson K, et al. Transcriptional profiling reveals suppressed erythropoiesis, up-regulated glycolysis, and interferon-associated responses in murine malaria. J Infect Dis. 2004;189(7):1245–56.

    Article  CAS  PubMed  Google Scholar 

  47. Essayan DM, Fox CC, Levi-Schaffer F, Alam R, Rosenwasser LJ. Biologic activities of IL-1 and its role in human disease. J Allergy Clin Immunol. 1998;102(3):344–50.

    Article  Google Scholar 

  48. Lukens JR, Gross JM, Calabrese C, Iwakura Y, Lamkanfi M, Vogel P, et al. Critical role for inflammasome-independent IL-1β production in osteomyelitis. Proc Natl Acad Sci U S A. 2014;111(3):1066–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Yao Y, Yu H, Liu Y, Xu Q, Li X, Meng X, et al. PSTPIP2 inhibits the inflammatory response and proliferation of fibroblast-like Synoviocytes in vitro. Front Pharmacol. 2018;4:9.

    Google Scholar 

  50. Cheung KT, Sze DM-Y, Chan KH, Leung PH-M. Involvement of caspase-4 in IL-1 beta production and pyroptosis in human macrophages during dengue virus infection. Immunobiology. 2018;223(4–5):356–64.

    Article  CAS  PubMed  Google Scholar 

  51. Senga S, Kobayashi N, Kawaguchi K, Ando A, Fujii H. Fatty acid-binding protein 5 (FABP5) promotes lipolysis of lipid droplets, de novo fatty acid (FA) synthesis and activation of nuclear factor-kappa B (NF-κB) signaling in cancer cells. Biochim Biophys Acta Mol Cell Biol Lipids. 2018;1863(9):1057–67.

    Article  CAS  PubMed  Google Scholar 

  52. Wu R, Chen F, Wang N, Tang D, Kang R. ACOD1 in immunometabolism and disease. Cell Mol Immunol. 2020;17(8):822–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Chen C-C, Liu H-P, Chao M, Liang Y, Tsang N-M, Huang H-Y, et al. NF-κB-mediated transcriptional upregulation of TNFAIP2 by the Epstein-Barr virus oncoprotein, LMP1, promotes cell motility in nasopharyngeal carcinoma. Oncogene. 2014;33(28):3648–59.

    Article  CAS  PubMed  Google Scholar 

  54. Jia L, Shi Y, Wen Y, Li W, Feng J, Chen C. The roles of TNFAIP2 in cancers and infectious diseases. J Cell Mol Med. 2018;22(11):5188–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Hu W-C. Microarray analysis of PBMC after Plasmodium falciparum infection: molecular insights into disease pathogenesis. Asian Pac J Trop Med. 2016;9(4):313–23.

    Article  CAS  PubMed  Google Scholar 

  56. Ockenhouse CF, Hu W, Kester KE, Cummings JF, Stewart A, Heppner DG, et al. Common and divergent immune response signaling pathways discovered in peripheral blood mononuclear cell gene expression patterns in Presymptomatic and clinically apparent malaria. Infect Immun. 2006;74(10):5561–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Henry SC, Daniell XG, Burroughs AR, Indaram M, Howell DN, Coers J, et al. Balance of Irgm protein activities determines IFN-γ-induced host defense. J Leukoc Biol. 2009;85(5):877–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Premzl M. Comparative genomic analysis of eutherian interferon-γ-inducible GTPases. Funct Integr Genomics. 2012;12(4):599–607.

    Article  CAS  PubMed  Google Scholar 

  59. Psallidas I, Stathopoulos GT, Maniatis NA, Magkouta S, Moschos C, Karabela SP, et al. Secreted phosphoprotein-1 directly provokes vascular leakage to foster malignant pleural effusion. Oncogene. 2013;32(4):528–35.

    Article  CAS  PubMed  Google Scholar 

  60. Liao W, Goh FY, Betts RJ, Kemeny DM, Tam J, Bay B-H, et al. A novel anti-apoptotic role for apolipoprotein L2 in IFN-γ-induced cytotoxicity in human bronchial epithelial cells. J Cell Physiol. 2011;226(2):397–406.

    Article  CAS  PubMed  Google Scholar 

  61. Hunt NH, Grau GE. Cytokines: accelerators and brakes in the pathogenesis of cerebral malaria. Trends Immunol. 2003;24(9):491–9.

    Article  CAS  PubMed  Google Scholar 

  62. Igarashi I, Suzuki R, Waki S, Tagawa Y-I, Seng S, Tum S, et al. Roles of CD4+ T cells and gamma interferon in protective immunity against Babesia microti infection in mice. Infect Immun. 1999;67(8):4143–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Wagner JL. Molecular Organization of the Canine Major Histocompatibility Complex. J Hered. 2003;94(1):23–6.

    Article  CAS  PubMed  Google Scholar 

  64. Kiniry BE, Hunt PW, Hecht FM, Somsouk M, Deeks SG, Shacklett BL. Differential Expression of CD8+ T Cell Cytotoxic Effector Molecules in Blood and Gastrointestinal Mucosa in HIV-1 Infection. J Immunol. 2018;200(5):1876–88.

    Article  CAS  PubMed  Google Scholar 

  65. Geiben-Lynn R, Kursar M, Brown NV, Addo MM, Shau H, Lieberman J, et al. HIV-1 antiviral activity of recombinant natural killer cell enhancing factors, NKEF-A and NKEF-B, members of the Peroxiredoxin family. J Biol Chem. 2003;278(3):1569–74.

    Article  CAS  PubMed  Google Scholar 

  66. Rénia L, Grau GER, Wassmer SC. CD8+ T cells and human cerebral malaria: a shifting episteme. J Clin Invest. 2020;130(3):1109–11.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Riggle BA, Manglani M, Maric D, Johnson KR, Lee M-H, Neto OLA, et al. CD8+ T cells target cerebrovasculature in children with cerebral malaria. J Clin Invest. 2020;130(3):1128–38.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Xia M, Gonzalez P, Li C, Meng G, Jiang A, Wang H, et al. Mitophagy enhances oncolytic measles virus replication by mitigating DDX58/RIG-I-like receptor signaling. J Virol. 2014;88(9):5152–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Miyashita M, Oshiumi H, Matsumoto M, Seya T. DDX60, a DEXD/H box helicase, is a novel antiviral factor promoting RIG-I-like receptor-mediated signaling. Mol Cell Biol. 2011;31(18):3802–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Haller O, Staeheli P, Kochs G. Interferon-induced mx proteins in antiviral host defense. Biochimie. 2007;89(6–7):812–8.

    Article  CAS  PubMed  Google Scholar 

  71. Kristiansen H, Gad HH, Eskildsen-Larsen S, Despres P, Hartmann R. The oligoadenylate synthetase family: an ancient protein family with multiple antiviral activities. J Interferon Cytokine Res. 2011;31(1):41–7.

    Article  CAS  PubMed  Google Scholar 

  72. Verhelst J, Hulpiau P, Saelens X. Mx proteins: antiviral gatekeepers that restrain the uninvited. Microbiol Mol Biol Rev. 2013;77(4):551–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Dong B, Zhou Q, Zhao J, Zhou A, Harty RN, Bose S, et al. Phospholipid scramblase 1 potentiates the antiviral activity of interferon. J Virol. 2004;78(17):8983–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Othumpangat S, Noti JD, McMillen CM, Beezhold DH. ICAM-1 regulates the survival of influenza virus in lung epithelial cells during the early stages of infection. Virology. 2016;487:85–94.

    Article  CAS  PubMed  Google Scholar 

  75. Carthagena L, Bergamaschi A, Luna JM, David A, Uchil PD, Margottin-Goguet F, et al. Human TRIM gene expression in response to interferons. PLoS One. 2009;4(3):e4894.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Ozato K, Shin D-M, Chang T-H, Morse HC. TRIM family proteins and their emerging roles in innate immunity. Nat Rev Immunol. 2008;8(11):849–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Bongen E, Vallania F, Utz PJ, Khatri P. KLRD1-expressing natural killer cells predict influenza susceptibility. Genome Med. 2018;14:10.

    Google Scholar 

  78. Dapat IC, Pascapurnama DN, Iwasaki H, Labayo HK, Chagan-Yasutan H, Egawa S, et al. Secretion of Galectin-9 as a DAMP during dengue virus infection in THP-1 cells. Int J Mol Sci. 2017;18(8):1664.

    Article  Google Scholar 

  79. Mengshol JA, Golden-Mason L, Arikawa T, Smith M, Niki T, McWilliams R, et al. A crucial role for Kupffer cell-derived galectin-9 in regulation of T cell immunity in hepatitis C infection. PLoS One. 2010;5(3):e9504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Berlutti F, Pantanella F, Natalizi T, Frioni A, Paesano R, Polimeni A, et al. Antiviral properties of Lactoferrin—a natural immunity molecule. Molecules. 2011;16(8):6992–7018.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Kurokawa C, Iankov ID, Galanis E. A key anti-viral protein, RSAD2/VIPERIN, restricts the release of measles virus in infected cells. Virus Res. 2019;263:145–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Perng Y-C, Lenschow DJ. ISG15 in antiviral immunity and beyond. Nat Rev Microbiol. 2018;16(7):423–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. González-Navajas JM, Lee J, David M, Raz E. Immunomodulatory functions of type I interferons. Nat Rev Immunol. 2012;12(2):125–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Haque A, Best SE, De Oca MM, James KR, Ammerdorffer A, Edwards CL, et al. Type I IFN signaling in CD8-DCs impairs Th1-dependent malaria immunity. J Clin Invest. 2014;124(6):2483–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Haque A, Best SE, Ammerdorffer A, Desbarrieres L, de Oca MM, Amante FH, et al. Type I interferons suppress CD4+ T-cell-dependent parasite control during blood-stage Plasmodium infection. Eur J Immunol. 2011;41(9):2688–98.

    Article  CAS  PubMed  Google Scholar 

  86. Kumar R, Bunn PT, Singh SS, Ng SS, Montes de Oca M, De Labastida Rivera F, et al. Type I Interferons Suppress Anti-parasitic Immunity and Can Be Targeted to Improve Treatment of Visceral Leishmaniasis. Cell Rep. 2020;30(8):2512–2525.e9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Kulinski JM, Tarakanova VL, Verbsky J. Regulation of antiviral CD8 T-cell responses. Crit Rev Immunol. 2013;33(6):477–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Welsh RM, Bahl K, Marshall HD, Urban SL. Type 1 interferons and antiviral CD8 T-cell responses. PLoS Pathog. 2012;8(1):e1002352.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. He X, Xia L, Tumas KC, Wu J, Su X-Z. Type I interferons and malaria: a double-edge sword against a complex parasitic disease. Front Cell Infect Microbiol. 2020;10:594621.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Pichyangkul S, Yongvanitchit K, Kum-arb U, Hemmi H, Akira S, Krieg AM, et al. Malaria blood stage parasites activate human plasmacytoid dendritic cells and murine dendritic cells through a Toll-like receptor 9-dependent pathway. J Immunol. 2004;172(8):4926–33.

    Article  CAS  PubMed  Google Scholar 

  91. Yao X, Wu J, Lin M, Sun W, He X, Gowda C, et al. Increased CD40 expression enhances early STING-Mediated type i interferon response and host survival in a rodent malaria model. PLOS Pathog. 2016;12(10):e1005930 Gazzinelli RT, editor.

    Article  PubMed  PubMed Central  Google Scholar 

  92. Yu X, Cai B, Wang M, Tan P, Ding X, Wu J, et al. Cross-regulation of two type I interferon signaling pathways in Plasmacytoid dendritic cells controls anti-malaria immunity and host mortality. Immunity. 2016;45(5):1093–107.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  93. Liehl P, Zuzarte-Luís V, Chan J, Zillinger T, Baptista F, Carapau D, et al. Host-cell sensors for Plasmodium activate innate immunity against liver-stage infection. Nat Med. 2014;20(1):47–53.

    Article  CAS  PubMed  Google Scholar 

  94. Wu J, Tian L, Yu X, Pattaradilokrat S, Li J, Wang M, et al. Strain-specific innate immune signaling pathways determine malaria parasitemia dynamics and host mortality. Proc Natl Acad Sci. 2014;111(4):E511–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. O’Connor RM, Long JA, Allred DR. Cytoadherence of Babesia bovis-infected erythrocytes to bovine brain capillary endothelial cells provides an in vitro model for sequestration. Infect Immun. 1999;67(8):3921–8.

    Article  PubMed  PubMed Central  Google Scholar 

  96. Sondgeroth KS, McElwain TF, Allen AJ, Chen AV, Lau AO. Loss of neurovirulence is associated with reduction of cerebral capillary sequestration during acute Babesia bovis infection. Parasit Vectors. 2013;6(1):181.

    Article  PubMed  PubMed Central  Google Scholar 

  97. Hermand P, Cicéron L, Pionneau C, Vaquero C, Combadière C, Deterre P. Plasmodium falciparum proteins involved in cytoadherence of infected erythrocytes to chemokine CX3CL1. Sci Rep. 2016;6(1):33786.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Sherman IW, Eda S, Winograd E. Cytoadherence and sequestration in Plasmodium falciparum: defining the ties that bind. Microbes Infect. 2003;5(10):897–909.

    Article  CAS  PubMed  Google Scholar 

  99. Totino PR, Lopes SC. Insights into the Cytoadherence phenomenon of Plasmodium vivax: the putative role of phosphatidylserine. Front Immunol. 2017;20:8.

    Google Scholar 

  100. Korir CC, Galinski MR. Proteomic studies of Plasmodium knowlesi SICA variant antigens demonstrate their relationship with P. falciparum EMP1. Infect Genet Evol. 2006;6(1):75–9.

    Article  CAS  PubMed  Google Scholar 

  101. Mueller I, Galinski MR, Baird JK, Carlton JM, Kochar DK, Alonso PL, et al. Key gaps in the knowledge of Plasmodium vivax, a neglected human malaria parasite. Lancet Infect Dis. 2009;9(9):555–66.

    Article  CAS  PubMed  Google Scholar 

  102. Maguire JD, Baird JK. The “non-falciparum” malarias: the roles of epidemiology, parasite biology, clinical syndromes, complications and diagnostic rigour in guiding therapeutic strategies. Ann Trop Med Parasitol. 2010;104(4):283–301.

    Article  CAS  PubMed  Google Scholar 

  103. Leisewitz A, Turner G, Clift S, Pardini A. The neuropathology of canine cerebral babesiosis compared to human cerebral malaria. Malar J. 2014;13(Suppl 1):P55.

    Article  PubMed Central  Google Scholar 

  104. Piercy SE. Hyper-acute canine babesia (tick fever). Vet Rec. 1947;59(44):612.

    CAS  PubMed  Google Scholar 

  105. Basson JG, Pienaar PA. Canine babesiosis: a report on the pathology of three cases with special reference to the “cerebral” form. J S Afr Vet. 1965;36(3).

  106. Jacobson LS, Reyers F, Berry WL, Viljoen E. Changes in haematocrit after treatment of uncomplicated canine babesiosis: a comparison between diminazene and trypan blue, and an evaluation of the influence of parasitaemia. J S Afr Vet Assoc. 1996;67(2):77–82.

    CAS  PubMed  Google Scholar 

  107. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.

    Article  CAS  PubMed  Google Scholar 

  108. Alexa A, Rahnenfuhrer J. topGO: Enrichment Analysis for Gene Ontology. In: Bioconductor version: Release (3.12); 2020. Available from:

    Google Scholar 

Download references


We would like to thank the CCR Sequencing Facility at Frederick National Laboratory for Cancer Research for performing the sequencing.


Funding for this project was from the National Research Foundation, South Africa, Grant number CPRR16042516163064 (AG, JS, AL). This research was supported in part by the Intramural Research Program of the NIH (RS, AB, SB, JL, HA). Open Access funding provided by the National Institutes of Health (NIH).

Author information

Authors and Affiliations



HA and AL designed the study and applied for funding. AG, JS, and AL cared for the animals, conducted the inoculations, and collected the clinical data and biospecimens. RS and SB extracted the RNA and submitted the material for sequencing. RS, AB, JL, HA, and AL analyzed and interpreted the data. RS, HA, and AL wrote the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Andrew Leisewitz or Hans Ackerman.

Ethics declarations

Ethics approval and consent to participate

The experimental protocol was approved by the University of Pretoria Animal Ethics Committee (V003–18). All methods were carried out in accordance with relevant guidelines and regulations.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

A summary of MultiQC data. Shown are the number of reads mapped (in millions), percent of these reads aligned to the canine reference genome, and percent GC content of these reads. Samples corresponding to low inoculum day 4 and high inoculum days 3 and 4 were not included in the mean.

Additional file 2: Table S2.

DEGs on every day in both inoculum cohorts. Genes that were differentially expressed on every experimental day in both cohorts (except day 1 in the low inoculum) are shown, along with their corresponding cluster.

Additional file 3: Table S3.

Pearson correlation values of erythroblast transcript trajectory with reticulocyte count. The expression level of each gene (in CPM) was correlated with the reticulocyte count (109/L) on each day. This allowed for correlation based on expression trajectory through time.

Additional file 4: Fig. S1.

Top 10 Gene Ontology (GO)- and Reactome-annoted pathways by FDR on each day in each inoculum cohort (A – low; B – high). x-axis: percentage of genes in pathway that were identified as differentially expressed on the respective day; y-axis: pathway name; point color: absolute number of genes associated with the pathway that were identified as differentially expressed on the corresponding day. Note that color does not indiciate up- or down-reulgation of genes in pathway.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Smith, R.L., Goddard, A., Boddapati, A. et al. Experimental Babesia rossi infection induces hemolytic, metabolic, and viral response pathways in the canine host. BMC Genomics 22, 619 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: