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

Background 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. Results 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. Conclusions 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. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07889-4.

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 a leading cause of infectious morbidity and mortality among the susceptible domestic dog population in South Africa (10). 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 (11,12). 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 proin ammatory cytokine storm (13). Additionally, some individual organ systems have been described in terms of pathology; others are currently being studied (14). 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 in ammation that lead to multiple organ dysfunction (15,16). However, conclusions drawn from clinical case series are elapsed from initial infection to presentation to clinic. 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 RNAsequencing, 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.

Results
Clinical time course correlated with infectious inoculum All 5 experimentally infected subjects developed clinical disease after intravenous inoculation with 10 4 or 10 9 parasites. The rst 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.

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 HCO 3 -, rising lactate) with respiratory compensation (low pCO 2 ) 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 ltration.
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. Creactive protein, a marker of in ammation, peaked on day 3 in the high inoculum cohort and on day 4 in the low inoculum cohort. Together, these biochemical ndings describe the acute onset of a severe hemolytic and in ammatory 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 rst four 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 rst three 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).
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 (link).

Global expression pro les 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. 8,001 genes passed this threshold and were used in subsequent analyses. Pairwise correlation of natural log-transformed counts for ltered host genes was high among a cluster composed of all day 0 samples and day 1 samples from the low inoculum cohort (mean r = .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 = .984).
Principal component analysis (PCA) of natural log-transformed counts of ltered host genes revealed tight clustering of all day 0 samples along with day 1 samples from the low dose cohort (Fig 3B). The rst 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 pro les 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 pro les, 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 pro les 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 were signi cantly upregulated throughout infection and recovery 8,001 ltered genes were tested for differential expression (log 2 fold change > |1|, FDR < 0.05), comparing expression on each experimental day versus day 0 ( Fig 4A) 12 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) Each cluster had a distinct temporal trend represented as the mean (95% CI) log 2 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 identi ed 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 Figs 5C and 5F.

Discussion
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. Prior knowledge of canine babesiosis pathophysiology was based on clinical observational studies in which duration of infection as well as the timing and size of the inoculum were variable and unknown. To address 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, re ected 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 identi ed 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 lifethreatening complications.

Hemolysis
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 (17,18), ALAS2, which catalyzes the rst and rate limiting step in heme synthesis in erythroid cells (19), and UROD, which codes an enzyme in the heme biosynthetic pathway (20). 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 (21). Mutations in these genes are also associated with hereditary spherocytosis (22). 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.
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 (23). 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 (24). Finally, GCLM, which has been shown to be essential to erythrocyte survival during oxidative stress (25) and whose de ciency is associated with hemolytic anemia (26), was also upregulated during infection.
Though there have not been extensive studies correlating gene expression with hemolysis in clinical malaria pathogenesis (27), the transcriptomic pro les reported here share many features in common with the host response to hemolytic anemia, including in sickle cell disease (18,(28)(29)(30)(31).
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 in ammatory response (32,33); GPR84, which serves as a receptor for free fatty acid (34); CD5L, which encodes a regulator of lipid synthesis that in turn regulates in ammatory response mechanisms and T cell activities (35); and ALOX15, which encodes a lipoxygenase that catalyzes the deoxygenation of polyunsaturated fatty acids and has known roles in red cell maturation and the in ammatory immune response (36). 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 (37), and FADS1, which encodes a fatty acid desaturase (38). Previous studies implicated the role of lipid droplets as regulators of the immune response to protozoan infection (39). 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 (40), was in cluster 3 in the low inoculum cohort and cluster 2 in the high. TKT, a key enzyme in the pentose phosphate pathway (41), 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 (42,43). Studies in malaria show glycolysis pathway genes prominently activated early during infection in mouse models and in humans (27,44). 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 in ammatory nature of the disease was re ected 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 (11). 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 re ect the differential expression of speci c 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-in ammatory cytokine (45); PSTPIP2, which encodes a protein that regulates IL-1β (46,47); CASP4, which regulates IL-1β synthesis in macrophages (48); and SIGLEC14, which encodes a receptor that enhances in ammasome 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 (49). 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 (50). Expression of TNFAIP2 is upregulated via the NF-kB pathway and mediates the in ammatory response (51,52). 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 (53,54). 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 signi cantly 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 (55), while IFGGC1 is a member of the IFN-γ-inducible GTPases that are involved in immune response against pathogens (56). SPP1 codes for a cytokine known to upregulate IFN-γ production (57), and APOL2 which mitigates the cytotoxic effects of IFN-γ (58). IFN-γ is a major pro-in ammatory cytokine that has both direct antiparasitic activity and immunoregulatory roles in the response to parasitic infection (59). IFN-γ has a role in protective immunity against infection by B. microti in mice (60) and is involved in the host immune response to malaria in both mice and humans (27,53,54). Our ndings 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 (61), were in cluster 1, indicating that they were signi cantly 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 log 2 fold change threshold for signi cance, excepting the low inoculum cluster 2 gene DLA-12). Additionally, GZMA and PRDX2, genes speci cally expressed by CD8+ T cells and associated with their activity (62,63), 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 and MHC-I and MHC-II play a role in the host immune response to malaria infection (27,54). However, recent studies in human cerebral malaria suggest that CD8+ T cells migrate to and sequester in the brain and contribute to mortality (64,65). 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).
Together, these transcriptomic data implicate the role of viral response genes and pathways in the host immune response to B. rossi infection.

Limitations
The inoculum of 10 9 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. 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 day 4 interventions including blood transfusions and anti-parasite drug. 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.

Conclusions
This study is the rst 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 ndings 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.

Animals
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 ve milliliters of whole venous blood were collected into EDTA vaccutanier tubes from dogs natural 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 con rm 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 quanti ed, and 10 mL of blood was collected into acid citrate dextrose solution and diluted with a culture medium used for invitro babesia culture to create infectious inocula, two doses of 10 4 parasites and three of 10 9 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 ve 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 TM , 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: 1collapsed, 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 re ll 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). A differential count was performed manually by an experienced laboratory technologist. An automated haematology was performed (Advia 2120i, Siemens, Germany) within 1 hour of blood collection. The remaining volume of blood was then centrifuged and aliquoted within 30 minutes 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 speci c 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 minutes, samples were run and the remaining serum volumes were aliquoted within 30 minutes 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 minutes (Rapidpoint 405, Seimens). 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 minutes 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 rst 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 on a single S1 owcell 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 les, then processes them through a sequence of programs to produce the desired output (i.e. raw counts matrix, QC report, differentially expressed genes, etc.). Once conFigd, the pipeline is executed on the NIH high performance computing cluster Biowulf. Pipeliner is publicly available on GitHub [https://github.com/CCBR/Pipeliner].
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 log 2 fold change > |1| were considered signi cantly 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 (87). Reactome-annotated pathway enrichment analysis was performed using the biomaRt (2.44.0) and reactome.db (1.70.0) packages, with signi cance 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 log 2 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 log 2 fold change on each day within each cluster was plotted along with 95% con dence intervals. Data analysis was conducted in R version 4.0.0.