Bovine macrophages transcriptome profiling reveals divergent responses to virulent and attenuated Babesia bovis strains

Bovine macrophages transcriptome profiling reveals divergent responses to virulent and attenuated Babesia bovis strains | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Bovine macrophages transcriptome profiling reveals divergent responses to virulent and attenuated Babesia bovis strains Valenzano Magalí Nicole , Liliana Alvarez , Beatriz Valentini , Silvina Wilkowsky doi: https://doi.org/10.1101/2025.09.16.676518 Valenzano Magalí Nicole 1 Instituto de Agrobiotecnología y Biología Molecular (IABIMO), Instituto Nacional de Tecnología Agropecuaria (INTA), Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET) , de Los Reseros y Nicolás Repetto s/n, Buenos Aires, Hurlingham B1686IGC, Argentina Find this author on Google Scholar Find this author on PubMed Search for this author on this site Liliana Alvarez 1 Instituto de Agrobiotecnología y Biología Molecular (IABIMO), Instituto Nacional de Tecnología Agropecuaria (INTA), Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET) , de Los Reseros y Nicolás Repetto s/n, Buenos Aires, Hurlingham B1686IGC, Argentina Find this author on Google Scholar Find this author on PubMed Search for this author on this site Beatriz Valentini 2 Laboratorio de Inmunología y Parasitología Veterinaria, EEA Rafaela, INTA , RN 34, Km 227, CC 22, 2300, Rafaela, Santa Fe, Argentina Find this author on Google Scholar Find this author on PubMed Search for this author on this site Silvina Wilkowsky 1 Instituto de Agrobiotecnología y Biología Molecular (IABIMO), Instituto Nacional de Tecnología Agropecuaria (INTA), Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET) , de Los Reseros y Nicolás Repetto s/n, Buenos Aires, Hurlingham B1686IGC, Argentina Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: wilkowsky.silvina{at}inta.gob.ar Abstract Full Text Info/History Metrics Preview PDF Abstract Babesia bovis is a tick-borne parasite of major economic impact in the livestock industry. Control strategies rely mainly on the use of acaricides and live attenuated vaccines. Comparative genomic analyses have shown no major differences between virulent and attenuated B. bovis strains, suggesting that studies on the host’s differential response may represent a key step toward clarifying the basis of disease severity and vaccine efficacy. In this study, we analyzed by RNA-seq the differential gene expression in bovine monocyte-derived macrophages after phagocytosis of erythrocytes infected with either the virulent B. bovis S2P strain or the attenuated R1A one. The results revealed a common transcriptional core response of several cellular processes largely centered on lymphocyte related functions, cytokine regulation, and adaptive immune signaling. In turn, the two strains elicited contrasting responses in bovine macrophages, where the virulent strain induced the enrichment of lymphocyte- and antiviral-related pathways, and the attenuated strain led to a stronger pro-inflammatory, chemotactic, and extracellular matrix remodeling signatures. Taken together, these data improved our understanding of the early transcriptional events that develop in macrophages in response to the phagocytosis of red blood cells containing contrasting B. bovis strains. This large dataset could be evaluated in further studies to better characterize the pathogenetic mechanisms of bovine babesiosis and to identify targets for host-directed therapeutic strategies. 1. Introduction Bovine babesiosis, primarily caused by Babesia bovis and B. bigemina , is a tick-borne disease of significant veterinary and economic importance, particularly in tropical and subtropical regions where the Ixodidae tick vectors are prevalent ( Jacob et al., 2020 ). In the mammalian host, these apicomplexan parasites invade and multiply only within bovine red blood cells. Particularly in adults, this leads to severe clinical manifestations that include high fever, hemolytic anemia, hemoglobinuria, and in many cases, cerebral babesiosis associated with sequestration of infected erythrocytes in the brain microvasculature. Of the two main species that cause bovine babesiosis, B. bovis infection can result in high mortality rates due to its capacity to adhere to capillary vessels and obstruct blood circulation ( Bock et al., 2004 ; Brown et al., 2006 ). The immune response to B. bovis involves both innate and adaptive mechanisms. The innate immune response involves the early activation of macrophages and dendritic cells, leading to the production of proinflammatory cytokines. Adaptive response is mainly mediated by opsonizing antibodies and a Th1-type cellular immune response, characterized primarily by the secretion of interferon-γ (IFN-γ) and tumour necrosis factor (TNF-α). These cytokines contribute to parasite control by enhancing splenic macrophage activation for the clearance of infected erythrocytes and by supporting the production of IgG2 ( Brown et al., 2006 ). Current vaccines against B. bovis are based on live attenuated parasites derived through serial passage in splenectomized calves or in vitro culture and have been used successfully in many countries ( de Waal & Combrink, 2006 ). Vaccination with these strains induces robust and long-lasting immunity that closely mimics the protection achieved after natural infection, yet without the occurrence of severe clinical signs ( Ristic, 2018 ). The proposed mechanism suggests that attenuation is the result of the selection of less virulent subpopulations that remain in circulation, while the more virulent parasites are sequestered in the microvasculature and thus lost during serial passages ( Baravalle et al., 2012 ). A previous study comparing the lipid composition of virulent and attenuated B. bovis strains demonstrated that differences in phospholipid profiles modulate the intensity and nature of the host immune response ( Gimenez et al., 2013 ). These findings suggest that signaling mechanisms underlying host–parasite interactions could be important determinants of disease outcome and immune activation. Therefore, understanding how the host differentially responds to virulent and attenuated strains at the molecular level is crucial to elucidating mechanisms underlying disease severity and vaccine efficacy. Recent advances in RNA sequencing (RNA-seq) have enabled the comprehensive characterization of host transcriptional responses to apicomplexan infections ( Menard et al., 2021 ; Terkawi et al., 2018 ). Moreover, transcriptomic analyses have already been successfully applied to study differential gene expression in peripheral blood mononuclear cells during B. bigemina infection, revealing important insights into immune activation and inflammation ( Martínez-García et al., 2025 ). However, to date, no studies have characterized the complete bovine macrophage transcriptome in the context of the immune response induced by B. bovis , nor compared responses to virulent and attenuated strains. In this work, we analyzed the host transcriptional response in bovine monocyte-derived macrophages infected with virulent and attenuated strains of B. bovis . By comparing the gene expression profiles induced by each strain, we aim to better understand the immune responses they trigger in the host and explore the mechanisms by which attenuated parasites elicit protective immunity without triggering severe pathology. 2. Materials and methods 2.1. B. bovis culture The Argentinian pathogenic S2P and the attenuated R1A strains were cultured in vitro using a microaerophilic stationary phase system with normal bovine erythrocytes, following the method outlined by Baravalle et al., 2012 . Parasites were utilized when parasitemia reached 10 %. 2.2. Bovine monocyte-derived macrophages (BMDM) and phagocytosis assays BMDM obtention was performed as described in Valenzano et al., 2025 . For the phagocytosis assays, 2 x 10 6 BMDMs per well were co-cultured with red blood cells (RBCs) infected with B. bovis S2P and R1A strains, and with non-parasitized erythrocytes (nRBCs) as a control, during 18 h at a multiplicity of infection (MOI) of 20:1. This time point was selected to capture the early transcriptional response of macrophages following the phagocytosis of B. bovis -infected erythrocytes ( Diotallevi et al., 2024 ; Ramírez et al., 2012 ). After infection and 3 washes with phosphate-buffered saline (PBS) buffer, BMDMs were harvested and processed. Using these experimental conditions, we found that approximately 70% of the macrophages ingested at least one RBC, confirmed by Giemsa-stained smears ( Valenzano et al., 2025 ). All assays were performed in triplicate. 2.3. RNA extraction For each replicate, cells were harvested using 1 mL of TransZol (TransGen Biotech, Beijing, China) per 6 × 10 6 cells and 200 µL of chloroform (Sigma-Aldrich, St. Louis, MO, USA) was added. The mixture was vigorously shaken for 30 seconds and samples were centrifuged for 15 minutes. The aqueous phase was carefully transferred to a new tube, and an additional 100 µL of chloroform were added. After vigorous shaking for 30 seconds, samples were centrifuged for 10 minutes. The upper aqueous phase was recovered, and 600 µL of isopropanol (Sigma-Aldrich, St. Louis, MO, USA) and 60 µL of sodium acetate (3 M, pH 5.2; Thermo Fisher Scientific) were added. Samples were mixed by gentle inversion and incubated at −70 °C for 2 hours. RNA was pelleted by centrifugation for 20 minutes, and the supernatant was discarded. The pellet was washed with 200 µL of cold absolute ethanol (Sigma-Aldrich, St. Louis, MO, USA), followed by centrifugation for 10 minutes. After discarding the ethanol, pellets were air-dried with the tube caps open for a few minutes. Finally, RNA was resuspended in 20 µL of diethylpyrocarbonate (DEPC)-treated water (Invitrogen, Carlsbad, CA, USA). All centrifugation steps were performed at 12,000 × g for the indicated time at 4 °C. Genomic DNA present in the samples was removed with a DNase I enzyme (1 U/µg RNA; Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions. Sample quality was assessed using the Fragment Analyzer system (Agilent Technologies, Santa Clara, CA, USA) and quantification was performed using both the NanoDrop™ spectrophotometer and the Qubit™ Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). The integrity of the RNA samples was also assessed by electrophoresis on a 1% agarose gel. 2.4. RNA-seq sample preparation and sequencing Libraries were prepared using the TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA) starting from 2 µg of pure RNA following the manufacturer’s instructions. Sequencing was performed at the National Center for Genomics and Bioinformatics, ANLIS-Malbrán, Argentina, on a NovaSeq 6000 system (Illumina, San Diego, CA, USA) to generate paired-end 2×150 bp reads. 2.5. RNA-seq data analysis The experimental design included three replicates for each condition. The quality analysis of the reads was performed using FastQC software ( Andrews, 2010 ). The reads were pre-processed using Trimmomatic ( Bolger et al., 2014 ) to remove the adapters. A comprehensive index was generated using the STAR aligner ( Dobin et al., 2013 ) from the Bos taurus genome ARS-UCD2.0 (GCF_002263795.3), incorporating splice junction annotations with the --sjdbOverhang 100 parameter to improve alignment accuracy of eukaryotic transcripts. Sequencing reads from each sample were mapped to the reference genome using STAR ( Dobin et al., 2013 ) with stringent parameters, allowing a maximum of 20 multiple alignments per read (--outFilterMultimapNmax 20) and were sorted by genomic coordinates (SortedByCoordinate). Transcript assembly and quantification were carried out with StringTie ( Pertea et al., 2015 ) in reference-guided mode, applying strand-specific fragment bias correction (-rf) to generate gene-level abundance estimates. Count matrices were prepared using the prepDE.py script provided by the StringTie pipeline. Differential expression analysis was conducted using DESeq2 ( Love et al., 2014 ). Differentially expressed genes (DEGs) were considered statistically significant if they exhibited a false discovery rate (FDR)-adjusted p-value 1. Heatmaps and principal component analysis (PCA) were constructed in R (version 2025.05.0) using the ggplot2 package ( Wickham, 2011 ). Gene overlap among datasets was visualized with the BioVenn tool ( Hulsen et al., 2008 ). Functional enrichment analyses were performed to identify overrepresented biological functions and pathways in the gene set. Gene Ontology (GO) enrichment analysis was conducted using the enrichGO function from the clusterProfiler R package to identify significantly enriched biological processes ( Yu et al., 2012 ). Pathway enrichment analysis was carried out using enrichKEGG, from the clusterProfiler R package, to detect overrepresented metabolic and signaling pathways according to the KEGG database (Kanehisa and Subramaniam, 2002). Additionally, Reactome pathway analysis was performed using the ReactomePA package to provide complementary insights into biologically relevant pathways ( Milacic et al., 2024 ). All analyses were conducted using default parameters, and pathways with an adjusted p-value < 0.05 were considered statistically significant. Enrichment results were visualized using ggplot2 ( Wickham, 2011 ) generating barplots to illustrate the most significantly enriched terms and pathways. 3. Results 3.1. Sequencing To investigate the transcriptional response of bovine macrophages upon exposure to erythrocytes infected with Babesia bovis strains S2P and R1A, as well as uninfected erythrocytes as a control, we generated paired-end sequencing libraries from the three experimental groups. For clarity, these conditions are hereafter referred to as S2P, R1A, and nRBC, respectively. Between 16151764 and 35361940 reads were obtained per sample ( Table 1 ). Quality assessment indicated consistently high-quality data, with Phred scores of 35.6-35.7. The GC content ranged from 51 % to 54 % with a uniform read length of 150 bp. These results confirmed that the sequencing data were of sufficient quality for downstream analyses. View this table: View inline View popup Download powerpoint Table 1: Summary of read metrics and quality control results. We performed a PCA on the three replicates per condition to assess sample clustering and variability (Supplementary Data 1). One of the control replicates (nRBCs_3) showed a clear deviation from the other two and was therefore excluded. Subsequent analyses were conducted using the remaining two replicates for this condition. 3.2. Differential expression analyses Across the three experimental conditions, 4679 protein-coding genes were aligned against the Bos taurus ARS-UCD2.0 reference genome which contains 21,667 protein-coding genes (Supplementary Data 2). Differential expression analysis revealed that in comparison with nRBC condition, 1654 genes were upregulated and 1618 were downregulated in S2P, meanwhile in R1A, 1496 were upregulated and 1769 were downregulated. The comparison between conditions S2P and R1A showed 700 upregulated genes in S2P relative to R1A, and 371 downregulated genes ( Figure 1 A ). A full list of differentially expressed transcripts for all three comparisons (S2P vs nRBC, R1A vs nRBC, and S2P vs R1A) is provided in Supplementary Data 2. Within protein-coding sequences, no substantial differences were observed in the number of differentially expressed genes (DEGs) between the S2P vs nRBC and R1A vs nRBC comparisons. Venn diagram analysis showed that 1866 transcripts (39.8 %) were differentially expressed in both S2P and R1A conditions vs nRBC, while showing no differential expression between S2P and R1A ( Figure 1 B ). Download figure Open in new tab Figure 1. Differentially expressed genes in BMDMs co-cultures with B. bovis –infected erythrocytes or with nRBCs (A) Number of up- and downregulated genes identified in each experimental condition. (B) Venn diagrams of differentially expressed protein-coding transcripts illustrating shared and unique expression patterns between virulent and attenuated B. bovis strains. 3.3. Common host pathways modulated by virulent and attenuated Babesia bovis strains Among the genes that were differentially expressed in macrophages upon phagocytosis of B. bovis –infected erythrocytes, regardless of whether the parasite was virulent or attenuated, we identified significant enrichment of common immune-related GO categories. These included terms associated with T cell differentiation and activation (e.g., GO:0045619, GO:0032946, GO:0050671), regulation of lymphocyte proliferation and migration (GO:0070665, GO:0045580, GO:0050871, GO:0042102), cytokine-mediated signaling pathways (GO:0045621, GO:0030888), and processes related to antigen processing and presentation (GO:0042130, GO:0031294, GO:0051251). Additional enriched categories involved regulation of adaptive immune responses (GO:0050864, GO:0032944), leukocyte-mediated cytotoxicity (GO:0050670, GO:0070663), and cell–cell adhesion during immune interactions (GO:0008645, GO:0042098, GO:0042129). Collectively, these results highlight that both virulent and attenuated B. bovis strains trigger a common transcriptional program in macrophages that is largely centered on immune-related pathways (Supplementary Data 3). 3.4. Strain-specific host responses to Babesia bovis infection Among the 1071 genes differentially expressed between macrophages co-cultured with the virulent and the attenuated B. bovis strains, 700 genes were found upregulated and 371 downregulated ( Figure 2 ). Download figure Open in new tab Figure 2: Heatmap illustrating DEGs between BMDMs that phagocytosed RBCs infected either the virulent strain (S2P_1, _2 and _3) or the attenuated strain (R1A _1, _2 and _3). The level of gene expression is represented by colors regarding log2FC: the higher score levels of expression are represented in red tones meanwhile the lower values are represented in blue tones. Among the genes upregulated in the S2P condition compared to R1A, we found numerous genes associated with lymphocyte activation and adaptive immune responses, including PRF1, GNLY, ZAP70, IL2RB, RORC, ICOS, ITK, TBX21, GATA3, FOXP3, IRF4, STAT4, and CXCR3. Several co-stimulatory receptors and signaling molecules were also enriched, such as CD3D, CD3E, CD28, CD27, and CD6. GO enrichment analysis highlighted terms related to lymphocyte and leukocyte activation, cytokine-mediated signaling, viral defense, and antigen receptor-mediated pathways ( Figure 3 A ; Supplementary Data 4). Moreover, processes such as response to type II interferon, negative regulation of viral replication, and hematopoiesis were significantly represented, underscoring a transcriptional program associated with immune response. Download figure Open in new tab Figure 3: Functional enrichment analysis of genes up-regulated in BMDMs that phagocytosed erythrocytes infected with the virulent versus attenuated strain. A-Bar plot showing the top enriched GO Biological Process terms. B-Top enriched KEGG pathways. Bars represent the enrichment significance, measured as −log10 adjusted p-value. Pathway analysis further supported this immune activation profile. KEGG analysis identified Th17, Th1, and Th2 enriched cell differentiation pathways, along with the T cell receptor signaling pathway, cytokine–cytokine receptor interactions, hematopoietic lineage commitment, and NF-κB signaling ( Figure 3 B). Similarly, Reactome analysis identified pathways such as cytokine signaling, interleukin-4, interleukin-13, and interferon α/β signaling, and also RUNX1/FOXP3-mediated development of regulatory T cells. These results point out a coordinated activation of the effector and regulatory arms of adaptive immunity (Supplementary Data 4). In contrast, a distinct transcriptional program characterized the genes downregulated in S2P vs R1A conditions. The downregulated genes included several regulators of cytokine and growth factor signaling (TGFBR1, VEGFC, IL18, IL1A, IL1B, IL23R, IL12RB2), chemokines (CXCL8, CXCL2, CXCL3, CXCL5, CXCL12), and structural and extracellular matrix components (COL1A1, COL1A2, COL3A1, COL6A1, COL6A2, LUM, SPARC). Notably, a high number of inflammatory mediators (S100A8, S100A9, S100A12, PTGS2, CCL2) and matrix metalloproteinases (MMP1, MMP2, MMP3, MMP12) were also significantly downregulated. In line with this, GO analysis revealed a marked enrichment of processes related to chemotaxis, leukocyte and granulocyte migration, extracellular matrix organization, actin cytoskeleton dynamics, angiogenesis, and inflammatory response ( Figure 4 A ; Supplementary Data 5). These categories highlight a difference in gene expression of the pathways involved in cellular motility, tissue remodeling, and innate immune effector mechanisms when macrophages ingest the two B. bovis strains. Download figure Open in new tab Figure 4: Functional enrichment of genes down-regulated in macrophages exposed to erythrocytes infected with the virulent versus attenuated strain A-Bar plot showing the top enriched Gene Ontology Biological Process terms. B-Top enriched KEGG pathways. Bars represent the enrichment significance, measured as −log10 adjusted p-value. Pathway analysis reinforced this observation ( Figure 4 B ). KEGG enrichment showed significant association with the IL-17 and TNF signaling pathways, cytokine–cytokine receptor interactions, chemokine signaling, focal adhesion, ECM–receptor interaction, and NF-κB signaling, were downregulated in the virulent condition compared to the attenuated context, indicating that multiple pro-inflammatory and tissue remodeling networks were altered. Consistently, Reactome pathways results included extracellular matrix organization, collagen degradation, integrin interactions, interleukin signaling (IL-4, IL-10, IL-13), and NFE2L2-regulated antioxidant responses, suggesting broad host differences in extracellular remodeling, adhesion, and cytokine-mediated inflammatory signaling between exposure to the virulent and attenuated strains (Supplementary Data 5). In summary, macrophages exposed to the two strains displayed distinct transcriptional programs, with the virulent S2P strain associated with enrichment of lymphocyte- and antiviral-related pathways. In contrast, the attenuated R1A strain induced a response characterized by stronger pro-inflammatory, chemotactic, and extracellular matrix remodeling signatures, highlighting clear differences in innate and adaptive immune engagement between the two strains. 4. Discussion RNA sequencing has become a powerful approach to investigate host-pathogen interactions at the cellular level, since it enables the identification of specific pathways and molecular mechanisms that are modulated during infection, which would be difficult to uncover through targeted approaches. This strategy has been extensively applied to study infections caused by other apicomplexan parasites, including Plasmodium falciparum ( Terkawi et al., 2018 ) and Toxoplasma gondii ( Menard et al., 2021 ), providing valuable insights into how host immune cells detect and respond to parasite invasion. In Babesia bovis , it is well established that attenuated strains can protect cattle against acute disease ( Ristic, 2018 ). Given that common virulence genes were not identified among all virulent B. bovis strains or their attenuated derivatives ( Lau et al., 2011 ) a plausible explanation is that they may generate distinct responses in the host, which in turn influence disease outcome. To test this hypothesis, we designed an experimental model that includes bovine macrophages as central mediators of the innate immune response against Babesia bovis . In natural infections, splenic macrophages recognize erythrocytes carrying damage-associated signals or surface proteins derived from the parasite and subsequently engulf them while presenting foreign antigens through MHC class II molecules. The in vitro model of BMDMs was successfully used by Valenzano et al., 2025 to identify novel B. bovis T-cell epitopes presented during infection, supporting its validity for dissecting host-parasite interactions. Although this interaction system may not perfectly mimic the behaviour of the bovine splenic macrophages, it is a feasible approach due to the lack of well-characterized continuous bovine macrophage cell lines. The comparison between the three conditions allowed the identification of 4,679 DEGs out of the 21,667 protein-coding genes annotated in the Bos taurus genome ( Liu et al., 2009 ). Among these, 1,866 genes were commonly differentially expressed in response to both the virulent and attenuated strains, indicating a core transcriptional response to B. bovis infection. Functional enrichment analysis revealed that these shared DEGs were significantly associated with immune-related GO categories. Collectively, these results indicate that both virulent and attenuated B. bovis strains elicit a conserved macrophage response strongly centered on lymphocyte activation, cytokine regulation, and adaptive immune signaling. The observation that most of the DEGs are related to immune responses, rather than to processes such as cell cycle regulation, proliferation, or migration, as has been reported for T. gondii ( Menard et al., 2021 ), is consistent with the intrinsic life cycle of B. bovis , which only resides and multiplies in the red blood cell. For this reason, our results are consistent with the mechanisms employed by macrophages to process and present parasite antigens, which explains why the transcriptional response observed for both strains is predominantly immune-centered. Among the immune-related genes reported in the natural response of bovines to infection by B. bovis , IFN-γ and TNF-α are particularly relevant molecules given their central roles in the Th1 response against the parasite ( Brown et al., 2006 ; Goff et al., 1998 ). In our dataset both genes were downregulated when comparing both S2P and R1A conditions vs nRBCs, suggesting that infection with B. bovis may suppress their expression in macrophages. In line with this, IL-12 and NF-κB 1 and 2, were consistently reduced in both infections compared to control, which may contribute to reducing the transcription of pro-inflammatory mediators. The observation that the expression of these cytokines is reduced in BMDMs that phagocytised infected erythrocytes is consistent with previous findings in murine macrophages during infection with the related parasites of the genus Plasmodium ( Wu et al., 2015 ). In this study, the authors showed that phagosomal acidification in macrophages suppresses the production of inflammatory cytokines, pointing out that dendritic cells constitute the major source of early proinflammatory mediators. Similar analysis with bovine dendritic cells co-incubated with B. bovis infected erythrocytes could shed light on whether this is a common trait between both apicomplexan parasites. Another possible explanation for the lack of an inflammatory response in macrophages could be the absence of opsonization of the infected erythrocytes. Previous studies in P. falciparum have shown that opsonized infected RBCs stimulate human macrophages to secrete proinflammatory cytokines such as TNF and IL-1β, whereas unopsonized infected RBCs induce a much weaker inflammatory response (Ozarslan et al., 2019; Zhou et al., 2012). In the Babesia cultures used in this study, normal bovine serum from negative donors was employed, meaning that no parasite-specific antibodies were present in the co-cultivation system. Thus, this hypothesis is plausible and requires additional investigation. The results described above were obtained after observing that one of the control replicates displayed a divergent expression profile, leading us to exclude it from further analyses. While several studies have demonstrated that reliable differential expression analyses can still be performed with two replicates when appropriate normalization and statistical methods are applied (Link et al., 2018; Srikumar et al., 2015), these results should be interpreted with caution. Importantly, the primary focus of this study was the comparison between the effects of virulent and attenuated strains on the host immune response, rather than contrasts with the uninfected control, which mitigates the impact of this exclusion on the overall conclusions. The comparison of bovine macrophages exposed to erythrocytes infected with virulent versus attenuated B. bovis , showed clear differences in the expression of genes, highlighting lymphocyte- and antiviral-related pathways for the former and pro-inflammatory, chemotactic, and extracellular matrix remodeling signatures for the latter. Consistently, the higher number of upregulated genes in the virulent strain condition may suggest that the host response is stronger with that phenotype. A similar pattern has been reported for the apicomplexan parasite T. gondii , where infection with a highly virulent strain elicited a more pronounced transcriptional response in host cells than infection with a less virulent counterpart ( Menard et al., 2021 ). However, a recent RNA-seq study of the bovine immune responses upon infection with B. bigemina showed that attenuated strain induces a stronger immune activation than the virulent one ( Martínez-García et al., 2025 ). These apparent discrepancies can certainly be explained by differences in the host cell model (animal vs in vitro cell cultures) and the timing of the response assessed. While our study focuses on macrophages at an early stage of infection, the B. bigemina study analyzed peripheral blood mononuclear cells at the peak of disease, at day 8 after inoculation. This reinforces the importance of carefully comparing host responses across multiple cell types and infection stages, as they may reveal complementary aspects of the mechanisms driving virulence and attenuation, thereby adding significant value to the overall understanding of host-parasite interactions. 5. Conclusions In this article, we investigated by RNA-seq the transcriptional responses in monocyte-derived macrophages that internalized RBCs infected with a virulent or an attenuated B. bovis strain. Aside from a common transcriptional response eliciting the upregulation or downregulation of several immune cellular processes, the different strain s also triggered a significantly different response in bovine macrophages, where several functional pathways appear activated or deactivated in response to the different phenotypes of the strains used. Collectively, these data improved our understanding of the early transcriptional events that occur in macrophages in response to B. bovis -infected RBCs and could be considered in further studies to better characterize the pathogenetic mechanisms of the disease as well as the comprehension of how live vaccines exert their effects. Supplementary Data captions Supplementary Data 1: PCA of RNA-seq samples. Each point represents one replicate (n = 3) for each experimental condition. Axes indicate the percentage of variance explained by the first and second principal components. Clustering replicates reflects the similarity of gene expression profiles within each condition. Supplementary Data 2: Lists of differentially expressed transcripts for all three pairwise comparisons: S2P vs uninfected, R1A vs uninfected, and S2P vs R1A. Supplementary Data 3: Gene Ontology (Biological Process) enrichment analysis of DEGs shared between S2P and R1A conditions compared to uninfected nRBCs. Supplementary Data 4: Functional enrichment of genes up-regulated S2P vs R1A A-Gene Ontology (Biological Process) enrichment analysis B-KEGG enrichment analysis C-Reactome analysis. Supplementary Data 5: Functional enrichment of genes down-regulated S2P vs R1A A-Gene Ontology (Biological Process) enrichment analysis B-KEGG enrichment analysis C-Reactome analysis. Acknowledgements The authors would like to thank Lic. Pablo Vera for his assistance with library preparation and sequencing, and Dr. Julio Di Rienzo for his helpful suggestions on statistical analyses. Abbreviations RBC Red Blood Cell IFN-γ Interferon Gamma TNF-α Tumor Necrosis Factor Alpha RNA-seq RNA Sequencing BMDM Bone Marrow-Derived Macrophage PBS Phosphate-Buffered Saline MOI Multiplicity of Infection DEPC Diethyl Pyrocarbonate FDR False Discovery Rate PCA Principal Component Analysis GO Gene Ontology References ↵ Andrews , S. ( 2010 ). FastQC - A quality control tool for high throughput sequence data . http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ . Babraham Bioinformatics . ↵ Baravalle , M. E. , Thompson , C. , Valentini , B. , Ferreira , M. , Torioni de Echaide , S. , Christensen , M. F. , & Echaide , I. ( 2012 ). Babesia bovis biological clones and the inter-strain allelic diversity of the Bv80 gene support subpopulation selection as a mechanism involved in the attenuation of two virulent isolates . Veterinary Parasitology , 190 ( 3–4 ). doi: 10.1016/j.vetpar.2012.06.037 OpenUrl CrossRef ↵ Bock , R. , Jackson , L. , De Vos , A. , & Jorgensen , W. ( 2004 ). Babesiosis of cattle . Parasitology , 129 . doi: 10.1017/S0031182004005190 OpenUrl CrossRef ↵ Bolger , A. M. , Lohse , M. , & Usadel , B. ( 2014 ). Trimmomatic: A flexible trimmer for Illumina sequence data . Bioinformatics , 30 ( 15 ). doi: 10.1093/bioinformatics/btu170 OpenUrl CrossRef PubMed Web of Science ↵ Brown , W. C. , Norimine , J. , Knowles , D. P. , & Goff , W. L. ( 2006 ). Immune control of Babesia bovis infection . Veterinary Parasitology , 138 ( 1–2 ), 75 – 87 . doi: 10.1016/j.vetpar.2006.01.041 OpenUrl CrossRef PubMed ↵ de Waal , D. T. , & Combrink , M. P. ( 2006 ). Live vaccines against bovine babesiosis . Veterinary Parasitology , 138 ( 1–2 ), 88 – 96 . doi: 10.1016/J.VETPAR.2006.01.042 OpenUrl CrossRef PubMed ↵ Diotallevi , A. , Bruno , F. , Castelli , G. , Persico , G. , Buffi , G. , Ceccarelli , M. , Ligi , D. , Mannello , F. , Vitale , F. , Magnani , M. , & Galluzzi , L. ( 2024 ). Transcriptional signatures in human macrophage-like cells infected by Leishmania infantum, Leishmania major and Leishmania tropica . PLoS Neglected Tropical Diseases , 18 ( 4 ). doi: 10.1371/journal.pntd.0012085 OpenUrl CrossRef ↵ Dobin , A. , Davis , C. A. , Schlesinger , F. , Drenkow , J. , Zaleski , C. , Jha , S. , Batut , P. , Chaisson , M. , & Gingeras , T. R. ( 2013 ). STAR: ultrafast universal RNA-seq aligner - Supplementary Data . Bioinformatics (Oxford, England) , 29 ( 1 ). ↵ Gimenez , G. , Belaunzarán , M. L. , Poncini , C. V. , Blanco , F. C. , Echaide , I. , Zamorano , P. I. , Lammel , E. M. , González Cappa , S. M. , & Isola , E. L. ( 2013 ). Babesia bovis: lipids from virulent S2P and attenuated R1A strains trigger differential signalling and inflammatory responses in bovine macrophages . Parasitology , 140 ( 4 ). doi: 10.1017/S003118201200193X OpenUrl CrossRef ↵ Goff , W. L. , Johnson , W. C. , & Cluff , C. W. ( 1998 ). Babesia bovis immunity: In vitro and in vivo evidence for IL-10 regulation of IFN-γ and iNOS . Annals of the New York Academy of Sciences , 849 . doi: 10.1111/j.1749-6632.1998.tb11046.x OpenUrl CrossRef PubMed ↵ Hulsen , T. , de Vlieg , J. , & Alkema , W. ( 2008 ). BioVenn - A web application for the comparison and visualization of biological lists using area-proportional Venn diagrams . BMC Genomics , 9 . doi: 10.1186/1471-2164-9-488 OpenUrl CrossRef PubMed ↵ Jacob , S. S. , Sengupta , P. P. , Paramanandham , K. , Suresh , K. P. , Chamuah , J. K. , Rudramurthy , G. R. , & Roy , P. ( 2020 ). Bovine babesiosis: An insight into the global perspective on the disease distribution by systematic review and meta-analysis . Veterinary Parasitology , 283 . doi: 10.1016/j.vetpar.2020.109136 OpenUrl CrossRef Kanehisa , M. , & Subramaniam . ( 2002 ). The KEGG database . Novartis Foundation Symposium , 247 . doi: 10.1002/0470857897.ch8 OpenUrl CrossRef PubMed Web of Science ↵ Lau , A. O. , Kalyanaraman , A. , Echaide , I. , Palmer , G. H. , Bock , R. , Pedroni , M. J. , Rameshkumar , M. , Ferreira , M. B. , Fletcher , T. I. , & Mcelwain , T. F. ( 2011 ). Attenuation of virulence in an apicomplexan hemoparasite results in reduced genome diversity at the population level . www.454.com ↵ Liu , Y. , Qin , X. , Song , X. Z. H. , Jiang , H. , Shen , Y. , Durbin , K. J. , Lien , S. , Kent , M. P. , Sodeland , M. , Ren , Y. , Zhang , L. , Sodergren , E. , Havlak , P. , Worley , K. C. , Weinstock , G. M. , & Gibbs , R. A. ( 2009 ). Bos taurus genome assembly . BMC Genomics , 10 . doi: 10.1186/1471-2164-10-180 OpenUrl CrossRef PubMed ↵ Love , M. I. , Huber , W. , & Anders , S. ( 2014 ). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biology , 15 ( 12 ). doi: 10.1186/s13059-014-0550-8 OpenUrl CrossRef PubMed ↵ Martínez-García , G. , Estrada , K. , Lira-Amaya , J. J. , Santamaria-Epinosa , R. M. , Lopez-Arellano , M. E. , Sciutto-Conde , E. L. , Rojas-Martinez , C. , Alvarez-Martínez , J. A. , Sánchez-Flores , A. , & Figueroa-Millán , J. V. ( 2025 ). Comparative Analysis of Immune Response Genes Induced by a Virulent or Attenuated Strain of Babesia bigemina . International Journal of Molecular Sciences , 26 ( 2 ). doi: 10.3390/ijms26020487 OpenUrl CrossRef ↵ Menard , K. L. , Bu , L. , & Denkers , E. Y. ( 2021 ). Transcriptomics analysis of Toxoplasma gondii-infected mouse macrophages reveals coding and noncoding signatures in the presence and absence of MyD88 . BMC Genomics , 22 ( 1 ). doi: 10.1186/s12864-021-07437-0 OpenUrl CrossRef PubMed ↵ Milacic , M. , Beavers , D. , Conley , P. , Gong , C. , Gillespie , M. , Griss , J. , Haw , R. , Jassal , B. , Matthews , L. , May , B. , Petryszak , R. , Ragueneau , E. , Rothfels , K. , Sevilla , C. , Shamovsky , V. , Stephan , R. , Tiwari , K. , Varusai , T. , Weiser , J. , … D’Eustachio , P. ( 2024 ). The Reactome Pathway Knowledgebase 2024 . Nucleic Acids Research , 52 ( D1 ). doi: 10.1093/nar/gkad1025 OpenUrl CrossRef PubMed ↵ Pertea , M. , Pertea , G. M. , Antonescu , C. M. , Chang , T. C. , Mendell , J. T. , & Salzberg , S. L. ( 2015 ). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads . Nature Biotechnology , 33 ( 3 ). doi: 10.1038/nbt.3122 OpenUrl CrossRef PubMed ↵ Ramírez , C. , Díaz-Toro , Y. , Tellez , J. , Castilho , T. M. , Rojas , R. , Ettinger , N. A. , Tikhonova , I. , Alexander , N. D. , Valderrama , L. , Hager , J. , Wilson , M. E. , Lin , A. , Zhao , H. , Saravia , N. G. , & McMahon-Pratt , D. ( 2012 ). Human Macrophage Response to L. (Viannia) panamensis: Microarray Evidence for an Early Inflammatory Response . PLoS Neglected Tropical Diseases , 6 ( 10 ). doi: 10.1371/journal.pntd.0001866 OpenUrl CrossRef PubMed ↵ Ristic , M. ( 2018 ). Babesiosis of domestic animals and man . In Babesiosis of Domestic Animals and Man . doi: 10.1201/9781351070027 OpenUrl CrossRef ↵ Terkawi , M. A. , Takano , R. , & Kato , K. ( 2018 ). Differential gene expression profile of human neutrophils cultured with Plasmodium falciparum-parasitized erythrocytes . Journal of Immunology Research , 2018 . doi: 10.1155/2018/6709424 OpenUrl CrossRef ↵ Valenzano , M. N. , Parker , R. , Nicastri , A. , Valentini , B. , Gravisaco , M. J. , Muñiz , X. F. , Eirin , M. E. , Montenegro , V. N. , Nielsen , M. , Ternette , N. , & Wilkowsky , S. E. ( 2025 ). Identification of T-cell epitopes of the intracellular parasite Babesia bovis by immunopeptidomic analysis of BoLA-II presented peptides . Vaccine , 61 . doi: 10.1016/j.vaccine.2025.127369 OpenUrl CrossRef ↵ Wickham , H. ( 2011 ). ggplot2 . Wiley Interdisciplinary Reviews: Computational Statistics , 3 ( 2 ). doi: 10.1002/wics.147 OpenUrl CrossRef ↵ Wu , X. , Gowda , N. M. , & Gowda , D. C. ( 2015 ). Phagosomal acidification prevents macrophage inflammatory cytokine production to malaria, and dendritic cells are the major source at the early stages of infection: Implication for malaria protective immunity development . Journal of Biological Chemistry , 290 ( 38 ). doi: 10.1074/jbc.M115.671065 OpenUrl Abstract / FREE Full Text ↵ Yu , G. , Wang , L. G. , Han , Y. , & He , Q. Y. ( 2012 ). ClusterProfiler: An R package for comparing biological themes among gene clusters . OMICS A Journal of Integrative Biology , 16 ( 5 ). doi: 10.1089/omi.2011.0118 OpenUrl CrossRef PubMed Web of Science View the discussion thread. Back to top Previous Next Posted September 16, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Bovine macrophages transcriptome profiling reveals divergent responses to virulent and attenuated Babesia bovis strains Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Bovine macrophages transcriptome profiling reveals divergent responses to virulent and attenuated Babesia bovis strains Valenzano Magalí Nicole , Liliana Alvarez , Beatriz Valentini , Silvina Wilkowsky bioRxiv 2025.09.16.676518; doi: https://doi.org/10.1101/2025.09.16.676518 Share This Article: Copy Citation Tools Bovine macrophages transcriptome profiling reveals divergent responses to virulent and attenuated Babesia bovis strains Valenzano Magalí Nicole , Liliana Alvarez , Beatriz Valentini , Silvina Wilkowsky bioRxiv 2025.09.16.676518; doi: https://doi.org/10.1101/2025.09.16.676518 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Immunology Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41913) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13372) Ecology (19889) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40365) Molecular Biology (17163) Neuroscience (88540) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15130) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9818) Zoology (2269)