Insights into the evolution, virulence and speciation of Babesia MO1 and Babesia divergens through multiomics analyses

Abstract

Babesiosis, caused by protozoan parasites of the genus Babesia, is an emerging tick-borne disease of significance for both human and animal health. Babesia parasites infect erythrocytes of vertebrate hosts where they develop and multiply rapidly to cause the pathological symptoms associated with the disease. The identification of new Babesia species underscores the ongoing risk of zoonotic pathogens capable of infecting humans, a concern amplified by anthropogenic activities and environmental changes. One such pathogen, Babesia MO1, previously implicated in severe cases of human babesiosis in the United States, was initially considered a subspecies of B. divergens, the predominant agent of human babesiosis in Europe. Here we report comparative multiomics analyses of B. divergens and B. MO1 that offer insight into their biology and evolution. Our analysis shows that despite their highly similar genomic sequences, substantial genetic and genomic divergence occurred throughout their evolution resulting in major differences in gene functions, expression and regulation, replication rates and susceptibility to antiparasitic drugs. Furthermore, both pathogens have evolved distinct classes of multigene families, crucial for their pathogenicity and adaptation to specific mammalian hosts. Leveraging genomic information for B. MO1, B. divergens, and other members of the Babesiidae family within Apicomplexa provides valuable insights into the evolution, diversity, and virulence of these parasites. This knowledge serves as a critical tool in preemptively addressing the emergence and rapid transmission of more virulent strains.

Keywords: Babesia MO1; Babesia divergens; Human babesiosis; multiomics; speciation.

Figure 1.

Life cycle of B. MO1 and B. divergens. A. Schematic representation of the life cycle of B. MO1 and B. divergens in vertebrate hosts (humans, cattle, cottontail rabbit) and tick vectors. B. Representative Giemsa-stained light microscopic images of the various stages of B. MO1 propagated in human erythrocytes in vitro. C. Representative Giemsa-stained light microscopic images of the various forms of B. divergens Rouen 87 grown in human erythrocytes in vitro. D. Growth of B. divergens Rouen 87 clones H2 and H6, and B. MO1 clones B12 and F12 in human RBCs in RPMI medium supplemented with 20% fetal bovine serum (FBS) or RPMI medium supplemented with 0.5% Albumax over a course of 4 days. Two independent experiments were performed in triplicates. E. Chromosomal organisation of Babesia MO1. PFGE shows the number and approximate sizes of bands in B. MO1 parental (PA) isolate: ∼5.7 Mb, ∼4.6 Mb, ∼3.5 Mb, ∼3.13 Mb and ∼2.35 Mb; the number and approximate sizes of bands in B. MO1 clones B12, H1, H6, and F1: ∼4.6 Mb (Chromosome I) ∼3.5 Mb (Chromosome II), and ∼2.35 Mb (Chromosome III) and B. MO1 clones F12 and A3: ∼5.7 Mb (Chromosome I), ∼3.5 Mb (Chromosome II), and ∼2.35 Mb (Chromosome III). The experiment was performed in biological duplicates. F. Chromosomal organisation of B. divergens. PFGE shows the number and approximate sizes of bands in B. divergens Rouen 87 parent, clones H6, A6, and H10: ∼4.3 Mb (Chromosome I and Chromosome II), and ∼2.1 Mb (Chromosome III) and B. divergens clones H2, C1 and C7: ∼4.3 Mb (Chromosome I), ∼4.1 Mb (Chromosome II), and ∼2.1 Mb (Chromosome III). Hansenula wingei and Schizosaccharomyces pombe DNA chromosomes were used as DNA markers. The manufacturer’s estimate of the sizes of chromosomes are indicated in Megabase pairs [13] on the right-hand side of panels E and F. The experiment was performed in biological duplicates.

Figure 2.

Apicoplast and mitochondrial genomes of Babesia MO1 and B. divergens. A-B. Graphical circular map of the apicoplast genome of B. MO1. and B. divergens Rouen 1987, respectively. C-D. Linear map of the mitochondrial genome of B. MO1 and B. divergens Rouen 1987, respectively. Orange arrows represent genes encoding proteins involved in the electron transport chain, including cox1, cox3, nad2, and cob. The genes encoding ribosomal RNA (rRNA) are depicted in pink colour. Different tRNA encoding genes are displayed in purple colour.

Figure 3.

Evolutionary analysis of Babesia MO1 genome. A. Upset plot depicting orthogroups between B. MO1 and other apicomplexans. In the upper panel, the percentage of annotated proteins for shared or unique ones from a given organism is presented. In the middle panel, the total number of unique or shared proteins from a given organism is depicted. The lower panel represents the intersection or uniqueness of a given species with horizontal bars at the left side, representing the total number of genes for a given species. B. Heatmap of ANI values between Babesia species and Theileria parva. Higher values (red colour) correspond to greater nucleotide similarity between the genomes.

Figure 4.

Circos synteny plots. The chromosomes of B. MO1 are illustrated on the right semicircle on all circular plots, and the chromosomes of the other organisms are on the left semicircle (A: B. duncani, B: T. parva, C: B. microti, D: B. divergens Rouen 87, E: B. bigemina, F: B. bovis); blue arcs indicate syntenies, red arcs indicate syntenies involved in a reversal; the intensity of the colour is proportional to the level of collinearity; the number after the species’ name refers to the chromosome number (when chromosomes are broken into pieces, fragments).

Figure 5.

Piroplasmida species phylogeny inferred from phylogenomic analysis. A. Species phylogeny obtained by super matrix and super tree phylogenomic approaches. All bootstrap values with super matrix were at 100%. Displayed clade support values are estimated in the case of super tree methods by concordance factors from the source trees of dataset #1/dataset #2. The position of Babesia MO1 was analyzed in relation to the two B. divergens isolates (highlighted in green colour in blue box). B. MO1 from the present study is in red (highlighted in blue box). Hepatocystis sp. (ex Piliocolobus tephrosceles 2019), Plasmodium falciparum 3D7 and P. gallinaceum 8A were taken as outgroup. B. Summary of the genetic exchanges between Piroplasmida species based on patristic distances. A matrix of patristic distances was calculated from the 2499 trees of dataset #1 for all pairs of species. Grey dot: median of the distribution. Comparisons between species of Clade VI, between B. MO1 and species of Clade VI, between B. MO1 and two strains B. divergens, and between two strains B. divergens are shown.

Figure 6.

Transcriptomic profile and epigenomic landmarks of B. MO1. A-B. Logarithms of the TPM counts in B. MO1 clones F12 (panel A) and B12 (panel B) were used as expression values for each gene across the three chromosomes using the R package ggplot2. C-D. RNA-seq data of B. MO1 clones F12 (panel C) and B12 (panel D) as normalized heat maps across the three chromosomes. Chromosomes were divided into 50 kb bins and the average of the log TPM of genes within each bin was calculated. n = 2 biologically independent samples. E. Comparison between epigenetic marks and gene expression in B. MO1 clones F12 and B12. Heat maps were built using normalized log₂H3K9me3 and H3K9ac read counts in addition to the RNA-seq TPM levels of each gene. Read counts for H3K9me3 and H3K9ac were normalized to millions of mapped reads and gene length, whereas TPM was determined by Stringtie. Genes were ranked from high to low TPM highlighting the correlation and anti-correlation between transcript abundance and the H3K9ac3 and H3K9me3 marks, respectively. F-G. Normalized H3K9me3 counts in multigene families, and other genes encoded by B. MO1 clones F12 (panel F) and B12 (panel G) (unpaired t-test with Welch’s correction, P < 0.0001) n = 2 biologically independent samples. H-I. Heterochromatin and euchromatin distribution across the three chromosomes of B. MO1 clones F12 (panel H) and B12 (panel I). Tracks correspond to H3K9ac3 ChIP [1], H3K9me3-ChIP (middle), and IgG control (bottom) and were normalized to millions of mapped reads. n = 2 biologically independent samples.

Figure 7.

Babesia MO1 3D-genome. A-B. Hi-C contact maps coupled with H3K9me3 ChIP-seq tracks (left) of B. MO1 clones F12 and B12 (10-kb kb bins). Tracks are scaled to chromosome lengths. C-D. 3D genome structures of B. MO1 clones F12 and B12 derived from the contact map interactions. Chromosomes one, two, and three correspond to green, pink, and blue sections respectively. Dark green and grey represent the telomeric regions and centromeres.

Figure 8.

Multi-gene families of B. MO1 and their chromosomal localization. A. Plot depicting the unique multigene families (UMGFs) in B. MO1. The blue bars depict the genes localized on one of the three chromosomes, whereas the yellow bars denote the genes found on stray contigs. B. Distribution of B. MO1 vesa1 and vesa2 genes on either chromosomes or stray contigs. C. Localisation of vesa genes and UMGFs members on the three B. MO1 chromosomes (genes localized on unplaced contigs are ignored). Genes denoted on the right side of a chromosome are on the positive strand, whereas those shown on the left side are on the negative strand.