Single-cell sequencing reveals unexpected genetic diversity among Bodo spp. flagellates and their bacterial endosymbionts

11 Bodo is a cosmopolitan genus of free living bacterivorous single-celled flagellates in 12 the class Kinetoplastea. Members of genus B odo are considered the closest free-living 13 relatives to the parasitic lineages Tr ypan os oma and Leis h mania , the causative agents of the 14 human diseases sleeping si ckness, Ch agas disease, and leishmaniasis. Currently, a single 15 genome exists for the one formally d escribed species in the genus, Bodo s altans . Previous 16 studies on B. s alt ans have show n t hat it is dependent on an endosymbiot ic bacterium from 17 the order Holosporales, “Candidat us Bodocaedibacter vickermanii”. U sing single cell-18 sequencin g, we i solat ed, sequenced, and assembled genomes f or seven uncult ured Bodo 19 spp. cells from a single freshwater sample from Royal Leamington Spa, UK. By usin g 20 comparat ive genomics, we show that t hese seven cell s represent t hree potentially novel 21 Bodo species and exhibit unexpected levels of diversity at the genome level. Our results 22 indicate that Small Subunit (SSU) rDNA sequencing, of ten used to classif y Bo do fla gella tes, is 23 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 2 insufficient for determining species delimitation in this genus. In addition, w e recovered a 24 Holosporales bacterium genome from all seven Bodo spp. cells. S urprisingly, these seven 25 bacterial endosymbionts also represent three p otentially novel species and one novel genus 26 of Ho losporales bacteria. This diversity would be indistinguishable in routin ely-used SSU 27 ribosomal DN A ( rDNA) metabarcodin g or bulk sequencing pipelines, thus demonstrating the 28 utility of using single- cell sequencin g to reveal the level o f genomic diversity within lineages 29 of microbial eukaryotes and their cobionts. 30 31

32 Single-cell sequencing, Kinetoplastea, Bodo saltans , protist, Holosporales, environmen tal 33 sequencin g 34 35

36 Bodo is a genus of heterotrophic free-living bi-flagellated pro tists common in fresh 37 and bracki sh water and soil. They are members of the Kinetoplastea, a class of parasitic and 38 free-living pro tists distingui shed by the presence of a large mass of mitochondrial DNA 39 known as kinetoplast D N A, or kDNA (1). Class Kinetoplastea is divided into tw o subclasses, 40 the Prokinetoplastina and the Metakinetoplastina, the latter containing fo ur orders, 41 Eubodonida, Parabodonida , Neobodonida and the Trypanosmat ida (2,3). Order 42 Trypanosomatida includes the parasites responsible for human diseases such as sleeping 43 sickness, Chagas’ disease, and leishmaniases. Phylogenies based on Small Subunit (SSU) 44 ribosomal DN A ( rDNA) and protein sequences place order Eubodonida, and its one genus, 45 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 3 Bodo, as the sister clade to the Trypanosmatida, and Bodo is considered the closest free-46 living lineage to the parasitic Trypanosomatida lineages (4–6). 47 Traditionally, Bodo species were identified and distinguished by morphology, e.g. cell 48 size and shape, length of the f lagella, position of the nucleus and kinetoplast, and 49 ultrastructural featu res visible by light and electron m icroscopy (7,8). Later, phylogenies 50 based on molecular data showed the genus Bodo was paraphylet ic and its members were 51 distributed in to three of the aforementioned four orders, with only one clade containing 52 three species, B . saltans , B. edax, and B. unci natus r emaining genus B odo (2,4,9). Further, it 53 is suggested that B. edax and B. uncin atus are not tru e species and are likely isolates of B. 54 saltans ( 2,10). Despite the important phylogenic position of genus Bodo, B . saltans remains 55 the only formally described species in the genus, and a single genome from bulk culture 56 currently exists for B. saltans strain Lak e Konstanz (11). 57 Holosporales are alphaproteobacterial t hat are widespread obligate endosymbionts 58 of eukaryotes, particularly protists (12,13). Holosporales form complex associations with 59 their eukaryotic hosts, including infectious parasitic species (14) , some that confer 60 competitive advan tages to their hosts ( 15), and some that may increase host f itness under 61 certain conditions (16,17) . Recent work has shown that B. s alta ns Lake Konst anz harbours 62 an endosymbiotic bacterium from the order Holosporales, “ Candida t us (C a.) 63 Bodocaedibacter vickermanii” (18). B. saltans appears dependent on its endosymbiont as 64 antibiotic treatmen t results in rapid cell death, and it is hypothesized that this dependence 65 is due to three pu tative addictive toxin/antitoxin systems encoded in the endosymbionts 66 genome (18) . These systems render the endosymbiont essential to its host as t hey encode a 67 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 4 long-lived toxin mo lecule alongside an antitoxin with a shorter half- l if e, and loss of the 68 antitoxin results in the activation o f the toxin and host cell dea th (19). 69 Here, we use single- cell sequencing t o assemble genomes from seven uncultured 70 cells from one environmental sample that were iden tified as belonging to genus Bodo. 71 Genomic comparisons show that the seven genomes form t hree clades, pot entially 72 represent ing three novel Bo do species that diverge significantly from B. s a l t ans and each 73 other. In addition, we recover a Holosporales bacterial genome f rom each raw assembly and 74 show that these putative endosymbionts also represent three poten tially novel species, 75 forming three clades that appear congruent with the phylogeny of the hosts. These results 76 highlight the utility of using single-cel l sequencing and comparative genomics to reveal the 77 diversity within and between uncultured populations of microbial eukaryotes and their 78 cobionts. 79 80

81 Single-cell sequencin g and genome assembly 82 Seven Bodo spp. cells w ere identified on tw o 96-well plates containing single cells 83 isolated by Fluorescence-Activated Cell Sorting (FACS) from an environmental sample. After 84 DNA amplification and short-read sequencing, between 20, 738, 413 – 48,141,972 paired-end 85 reads were generated for each cell (Table 1). The cells were named after the well in which 86 they were deposited and are hereafter referred to as A8, A10, B2, B7, F10, G10 and H10. 87 Each raw assembly w as curated by bi nning and taxonomic assignment to remove scaffolds 88 of bacterial and non-target origin, and the curated eukaryotic assemblies r ange in size f rom 89 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 5 29,029,105 base pairs ( bp) to 36,294, 242 bp compared to t he bulk- cult ure assembly for B. 90 saltans lake Konstanz at 39,862,120 bp (Table 1). All assemblies are fragmented, each in 91 more than 2,619 contigs with contig N50s ranging from 22,220 bp to 31, 876 bp (Table 1). All 92 assemblies range from 54.84 % t o 56. 02 % G C, a hi gher % GC than B. saltans lake Konstanz 93 at 51.79 % (Table 1). Protein annotat ions of each genome result ed in 11,121-13,166 protein 94 annotations per single cell genome, 30-40 % less than B. saltans lake Konstanz at 18,190 95 proteins ( Table 1) . Completen ess anal ysis of the protein sets using BUSCO v5 with the 96 Euglenozoa odb10 database gave 72. 3 %-80.8 % complete BUSCOs per assembly, compared 97 to B. saltans lake Konst anz at 88.5 % complete BUSCOs (Table 1). 98 99 Ce l l Identifi e r raw pa i r ed - en d rea ds A s sem b ly s ta t i s t ics No. p r otei ns a nno ta te d B US C O - Eu g l eno zoa o db 1 0 Len gt h No. con ti gs N50 GC % Co m pl et e % S D F M A8 20,814,8 14 30,867,202 2,619 26,848 54.94 11,505 72.3 72 0 5 22 A10 25,608,2 25 31,501,931 2,828 25,310 54.9 11,858 76.2 75 1 3 21 B2 20,738,4 13 36,294,242 3,436 26,706 55.02 13,166 80.8 80 1 6 13 B7 48,141,9 72 30,770,982 3,056 22,521 56.02 11,736 77.0 76 1 4 19 F1 0 21,853,2 21 29,029,105 2,665 24,738 56.02 11,121 75.4 75 0 6 18 G10 33,945,9 03 33,462,613 3,452 22,220 55 12,168 74.6 73 2 5 20 H10 21,728,8 38 34,908,636 3,605 22,252 54.97 12,867 80.0 79 2 5 15 B. saltans 39,862,120 2,256 31,876 51.79 18,190 88.5 112 3 12 3 100 Table 1 - Bodo spp. sequencing and assembly statistics. Number of pair ed-end r eads 101 generated for each single cell assembly, and general assembl y and annotation statistics. The 102 assembly statistics for B. s alt ans lake Konstanz genome are show n as a comparison. BU SCO 103 completeness f or protein sets, S = complete single copy, D = complete dup licated, F = 104 fragmented, M = missing. Euglenozoa odb10 includes 130 BU SCOs. 105 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 6 The seven Bodo cells form three clades 106 BLASTn searches against the GeneBank Nucleotide (nt ) database using the SSU rDNA 107 sequences from each assembly returned three B odo spp. as top hits (Table 2). For A8 a nd 108 A10, the top hit is to ATCC isolate Bo do edax (ATCC30903), originally isolat ed from the Czech 109 Republic (5). For B7 and F10, the best hit is to B. s alt ans strain HFCC309 i solated from a 110 eutrophic pond in G ermany (20). For B2, G10 and H10, the top hit is to a B. saltans isolate 111 from soil in Malaysia ( 10). In all cases, the query sequences are more than 99% identical to 112 their top hits. Figure 1A shows the pairwise identities of the SSU rDN A nucleotide sequences 113 from th e seven single cell assemblies. A8 and A10 are 100% identical. B7 and F10 are 99.82% 114 identical, with one m ismatch between the two sequences, a one nucleotide indel o f a G 115 residue in a string of G residues, however, due to the fragmented nature o f the assemblies 116 the SSU rDNA sequence from F10 is truncated, limiting f ull comparison of these sequences 117 (Table 2 and Additional File 1). B2 and H10 are 100% identical, with G 10 99.92 % identical to 118 B2 and H10 w ith a one nucleotide polymorphism, a C to T transition in G 10. All th ree SSU 119 rDN A sequences from B2, H10 and G10 are truncated, missing ~700 nucleotides f rom their 120 5’ ends (Table 2 and Additional File 1). A8/A10 differ fro m B2/G10/H10 by one nucleotide, 121 an A/G transition, while B7/F10 differ more substantially from all the others, with pairwise 122 identities ~98% (Figure 1A) . 123 A Maximum-Likelihood ( ML) tree of SSU rDN A sequences of Bodo and related 124 species available in G eneBank places our seven cells onto two clades, all within the B. 125 saltans lineage (Figure 1B). A8/A10 are sister to B2/G 10/H10 and are part of a radiation 126 that includes B. edax , B. unic in at us and three o ther B. s alt a ns isolates. B7/F10 are si ster to 127 this is clade, along w ith two other B. s alt ans isolates. B . s altan s lake Konst anz branches with 128 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 7 other B. s alt a ns isolates on a separate clade that diverges closer to the ro ot of th e tree 129 (Figure 1B). The seven single cell Bod o genomes have SSU rDNA pairwise identities ranging 130 between 94.7 and 95.7 % with B. saltans lake Kon stnaz and appear closer to t he SSU rDNA 131 sequence f rom B. edax, with pairwise identities ranging between 98.0 and 99.6 % identity 132 (Figure 1B) . 133 134 Ass e m bl y SS U leng th To p hit Bl as t n su bj ec t ac ce ss i o n T op hit B l as t n su bj ect or g anis m T op hit B l as t n quer y cov er age T op hit B l a s t n pai rwis e id enti t y T op hit B l as t n bi t s c or e A8 2,147 AY028451 Bodo edax 99.81% 99.60% 3829.84 A10 2,147 AY028451 Bodo edax 99.81% 99.60% 3829.84 B7 2,147 DQ207572 Bodo saltans 98.04% 99.30% 3724.34 F10 1,124 DQ207572 Bodo saltans 97.86% 99.40% 1951.63 B2 1,393 AY490226 Bodo saltans 99.86% 99.70% 2492.64 G10 1,263 AY490226 Bodo saltans 100% 99.80% 2266.32 H10 1,409 AY490226 Bodo saltans 99.86% 99.80% 2521.49 135 Table 2 – BLAST results for SSU rDNA sequences. Top hits resulting from a BLASTn search of 136 the nt database f or the SSU rDNA sequence recovered from each single cell assembly. 137 138 Extensive genomic diversity among t he Bodo spp. 139 We next built a Maximum Likelihood (ML) tree from 488 single-copy orthogroups 140 shared by B. s a l t ans lake Konstanz and our seven single cell Bodo asse mblies, using 141 Per kinsela sp. as an outgroup ( Figure 2A). Here, the seven single-cell genomes split into 142 three clades with A8/A10 form ing a si ster clade to B2/G 10/H10 and B7/F10 forming another 143 clade branching closer to the root o f the tree. B. sal t a ns lake Konstanz forms its own clade 144 closest to the root of the tree. N ext, we cal culated average pairwise amin o acid identities 145 (AAI) for all eight Bodo assemblies, wi th each pairwise comparison averagi ng AAI over a 146 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 8 minimum of 7,419 and a maximum of 10,181 reciprocal best hits (Addition al File 2). The AAI 147 within each clade is high (>97%) but varies substantially betw een the clades (Figure 2B). B. 148 saltans lake Konstanz i s equally di stant from all seven single cell assemblies, with pairwise 149 AAI between 51.3-51.7 % in all comparisons. B7/F10 are equally distant from A8/A10 and 150 B7/G 10/H10 with pairwise AAI between 62.1-62.6%. A8/A10 and B7/G10/H10 appear 151 closer, w ith pairwise AAI between 84.8-85.0%. 152 We also compared t he Average N ucleotide Identity (ANI) and the Aligned Fraction 153 (AF) of each Bodo spp. genome using the too l skani (21) (Figure 2C). Here, a result of 0.00 154 for AN I and/or AF suggests that sequence similarity at the nucleotide level is too low for 155 pairwise comparison using this method (21), e.g. B. s a l t ans lake Konstanz is too distant to 156 any of the single-cell Bod o genomes to compare using this metric, with values of 0.00 in 157 every pairwise comparison. However, ANI and AF are relat i vely high within each of the three 158 single-c ell Bodo clades ( >98 % AN I and > 80 % AF). Between the three clad es B7/F10 is too 159 distant from the ot hers to compare, while A8/A10 and B2/G10/H10 show an ANI of ~84%, 160 but an AF of ~14-16% indicating that only a small portion of the genome is close enough for 161 comparison using AN I . The AF values shown in Figure 2C are the average of two AF values 162 calculated f or each pairwise compari son. A matrix with bo th AF values calculated for each 163 pairwise comparison is in Additional File 3. 164 Next , w e compared how many orthogroups are shared among the single-cell 165 assemblies and B. s alt ans lake Konstanz ( Figure 2D). The majority of ortho groups, 55.5% are 166 shared by all eight assemblies. The two clades that branch most closel y in the phylogenetic 167 trees, A8/A10 and B2/G10/H10 share the highest percentage of orthogroups, 83.4%. B7/F10 168 share 66.1% with A8/A10 and 68% wi th B2/G 10/H10. B. saltans lake Konst anz shares 62.8 % 169 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 9 with B7/F10. 61. 4 % wit h A8/A10 and 63.7% wit h B7/G10/H10. Each clade has a small 170 percentage of orthogroups that are unique, w ith B. s alt ans lake Konstanz having the most 171 unique orthogroups. 172 To compare t he functional complim ent of each genome, we mapped each protein 173 set to the Pfam database. In each case, less than 50% of proteins in each annotation contain 174 one or more Pfam domains (Figure 2E) . W e tabulat ed the Pfams found in each prot ein 175 annot at ion and ran a Principal Components Analysi s (PCA) on that table. The PCA separates 176 the genomes into t hree clust ers congruent with t he phylogeny. Principle Component ( PC) 1 177 separates B. s alt ans lake Konstanz f rom all others. While PC 2 separates B2 and F10 f rom 178 A8/A10 and B7/G 10/H10, w hi ch clust er together (Figure 2F, Additional File 4). 179 Taken together, the phylogeny, AAI, AN I, AF and Pf am clustering indicate that the seven 180 single cell Bodo assemblies form t hree clades: clade A8/A10, clade B7/F10 and clade 181 B2/G 10/H10. All three clades are potentially novel Bodo species, and all appear equally 182 distant from B. salta ns lake Konstanz. 183 184 All Bodo spp. harbour Holosporales bacterial endosymbionts 185 Previous st udies show that B . saltans lake Konst anz is dependent on an 186 endosymbiotic alphaprot eobacterium f rom t he order Holosporales, “ Ca. B. vickermanii” 187 (18). We went back t o our binned raw genome assemblies and were able t o identify a single 188 bin in five of the genomes, 2 bins in B2 and 4 bins in B7 that were > 80 % complete 189 according to CheckM (Additional File 5). CheckM analysi s of these genomes showed that all 190 were at least 89 % complete, and all but one had contamination lev els < 5% (Table 3). 191 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 10 Ass e m b l y C h e c k M An not ati o n Cell Le n g t h ( b p) No . c ont ig s N50 GC % Com p- le t en e ss % C ont - am ina t ion str ai n he te r o- ge neit y no. ge nes no. tR NA s rDNA s A8 1,591,871 98 28,883 37.82 90.32 0 0 1,591 39 (20) 23s, 16s A1 0 1,758,263 139 26,374 37.27 90.32 0 0 1,591 41 (19) 0 B2 2,174,152 145 31,250 34.85 93.55 26.34 10.34 2,011 55 (19) 0 G1 0 1,604,499 119 27,931 38.14 89.01 2.2 100 1,363 39 (20) 23s, 16s H1 0 1,693,415 108 27,858 37.65 94.62 0 0 1,463 39 (17) 0 F10 1,279,012 83 33,152 41.61 94.62 0 0 1,234 35 (18) 0 B7 1,292,737 53 55,440 41.89 95.6 0 0 1,220 38 (19) 0 “Ca. B. vi cker ma nii” 1,391,311 1 NA 40.61 94.62 0 0 1,214 192 Table 3 – Assembly statistics, CheckM results for Holosporales genomes. G eneral assembly 193 and annot at ion st atist ics for each Holosporales genomes. No. tRNAs shows the tot al number 194 of tRN As, and the number of dif feren t amino acids in brackets. Values from “ Ca. B. 195 vickermanii” are shown where appropriat e for comparison. 196 197 GTDB-TK classified one bin from B7 and the bin in F10 in th e order Ho losporales, 198 family “ Ca. Paracaedibacteraceae”, genus U BA6184, as the same genus as “ Ca. B. 199 vickermanii” (Additional File 5). Six other bin, including a single bin from B2 were classif ied 200 in order Paracaedi bacterales , family U BA11393, genus J AGOTX01, whi ch currently contains 201 a single bacterial met agenom e assembled genome (MAG) isolat ed from wastewater ( 22) 202 (Additional Fil e 5) . The additional bin s f rom B2 and B7 were not classified as 203 alphaproteobacterial and were no t analysed further (Additional File 5). Order 204 Para c aedibactera l es is considered a heterotypic synonym of order Holosporales ( 12). 205 Recently it has been proposed that due to phylogenetic nesting of order H olosporales, it be 206 down-ranked to suborder Ho losporin eae, which w ould include the families H ol o s por ac eae , 207 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 11 “ Caedimonidace a ” and “Ca. Paracaedibacteracea” ( 13,23). U nder this ranki ng, all our 208 genomes would fall w ithin suborder H olos por ineae , two falling in fam ily “Ca. 209 Paracaedibacteracea” and the rest as unclassified H olos por ineae . 210 To better place the bod o Ho losporlaes, we built a ML tree fro m 24 single-copy 211 ort hogroups shared among Holosporales genomes and MAG s available in GeneBank and 212 using two alphaprot eobact erial as an out group (Figure 3A). This t ree places the seven 213 genomes into two lineages. F10/B7 branch on a clade wit h “ Ca. B. vickermanii” and two 214 MAG s classified as family “ Ca. Paracaedibacteracea” assembled from metagenomes of 215 wast ewat er samples (24, 25). A8/A10 and B2/H10/G 10 form another clade with t hree MAG s 216 assembled from m etagenomes of water samples, including wastewater (22,26, 27). “ Ca. 217 Finniella inopinata”, an endosymbiont of the amoebaflagellat e V iri dirap t or invadens 218 branches at the base of this clade (23,28). W e calculated pairwise AAI for all the genomes in 219 the tree (Figure 3A, Add itional File 6). The pairs B7 and F10, A8 and A10, and G10 and H10 220 all have a pairwise A AI greater than or equal to 95 %, considered a cutoff for species 221 delimitation in the literature (29,30) , indicating that they are most likely the same species. 222 B2 is more distant from G10 and H10 with AAI values of 58.8 and 59.62, respectively 223 (Additional file 6). However, CheckM result s show that t his genome cont ains contaminat ion, 224 which could be interfering with these analyses. The AAI boundary for genus delimitatio n 225 varies depending on the taxa under consideration (31–34). Here, we used a cutof f of 55%, as 226 was used in Midha et al. (2021) w hen describing “ Ca. B. vi ckermanii”. This cutoff place clade 227 B7/F10 as a novel species in the same genus as “ Ca. B. vi ckermanii”, while clade A8/A10 and 228 G10/H10 form t wo novel species in a novel genus. 229 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 12 We also compared AN I and AF values f or all pairwise comparisons between the Bodo spp. 230 endosymbionts and “ Ca. B. vickermanii” (Figure 3B). AN I and AF are hig h w ithin clades (ANI 231 >97%, AF >85%), but, bet ween the clades, only A8/A10 and H10/G10 are close enough to 232 compare, with AN I values of 85.42-86.31 % and AF at 27.73-29. 52 %, indicating that like 233 their hosts, these genomes share high sequence similarity only over a small proportion o f 234 their genomes. The B2 genome shows high ANI ( > 99% ) w it h the genomes from H10 and 235 G10, but low AF (<25%) indicating that this genome is only similar to those genomes over a 236 small proportion of the genome. How ever, contamination in this genome may mean that 237 this result is an artifact. The AF values shown in Figure 3B are the average of two AF values 238 calculated f or each pairwise compari son. A matrix with bo th AF values calculated for each 239 pairwise comparison is in Additional File 7. 240 Finally, we compared t he number of shared ort hogroups between the Holosporales 241 genomes and “ Ca . B. vickermanii” (Figure 3C). 23. 2% of orthogroups are shared by all 242 genomes. The genomes t hat are closest in t he phylogeny, A8/10 and B2/G10/H10, share t he 243 highest proportion of orthogroups at 52.6 %. G enomes B2/G10/H10 have the highest 244 proportion of un ique orthogroups at 12.35, ho w ever, the contamination in the B2 genome 245 is likely to be inflating this number. 246 In summary, B. saltans and each novel Bodo clade presented in this paper harbour a 247 unique species of Holosporales endosymbiont . These seven endosymbionts split into t wo 248 genera (genus “ Ca. Bodocaedibacter” and a novel genus) t hat appear congruent with the 249 phylogeny of the host and further support the conclusion that these are three novel Bodo 250 species. 251 252 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 13 The Bodo Holosporales differ in their metabolic capacity. 253 We compared the secretion systems and met abolic modules present in t he bact erial 254 genomes and compared them to “ Ca. B. vickermanii” (Figure 4) . Genus “ Ca . 255 Bodocaedibacter” (B7/F10) and the novel genus described in this study (A8/A10 and 256 B2/G 10/H10) encode a SEC-SRP and a Type VI secretion system (18). In addition, the n ovel 257 genus encodes two out the fou r proteins of the TAT secretion system, while a partial Type IV 258 conjugal transfer pilus assembly protein systems i s spottily distributed across the taxa. 259 All seven genomes have a partially complete (>= 66% complete) pat hway for lysine 260 biosynthesis, which is absent in all seven hosts ( <=22% complete) (Figure 4 ). In addition, the 261 novel genus encodes several met abolic pat hways absent in genus “ Ca. Bodocaedibacter”, 262 notably including heme biosynthesi s ( Figure 4). Heme is an essential nutrient, and to date, 263 no complete heme pat hway has been described for a Kinet oplast ea member, which either 264 require an exogenous source of heme, or contain endosymbiotic bact eria that produce it 265 (35). Figure 4 only shows KEGG modules that were at least 50% complet e in one or m ore 266 genomes. A full list of all the KEGG modules, including those at less than 50% completeness 267 is presented in Additional File 8. 268 Midha et al. (2021) reported that “ Ca. B. vickermanii” encodes three putative 269 polymorphic t oxin/anitoxin system s which may be responsible for its host ’s dependency 270 (18). Polymorphic toxin/antitoxin systems are typically composed of a large multi-domain 271 protein contain ing an N - terminal secretion signal f ollowed by a toxin domain, a protective 272 immunit y protein, and multiple casset tes encoding alt ernative toxic domains and associat ed 273 immunit y proteins (19). The t oxin genes can cont ain a homologous repeat region t hat 274 enables recombination betw een the full-length toxin gene and the alternative toxins. We 275 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 14 identified th ree genes belonging to one putative toxin/antitoxin system in both B7 and F10, 276 but none in the other five genomes from the novel genus. These genes show some of t he 277 characteristics of a polymorphic toxin /antitoxin system including a homolo gous region. 278 How ever, the pu tative large mul tido main protein in B7 d oes not contain an N-terminal 279 signal pept ide, and the one put ative orphan module in bot h B7 and F10 contains a t oxin but 280 no cognate anti-toxin (Additional File 9). Therefore, the l ocus’ identity and f unction as a 281 polymorphic toxin/antitoxin system is tentative. In addition, the incompleteness of the 282 single cell genomes presented here raises the possibility that key proteins and/or pathways 283 are missing from these assemblies due to under sampling. 284 285 Genus Bodo encodes many unique protein families, but core metabolism is generally 286 conserved with parasitic Kinetoplastea 287 We next compared ort hogroups distribut ed across the Kinet oplast ea lineage (Figure 288 5). For each kinetoplastid genome we tallied the total numb er of orthogrou ps in that 289 genome, and the number of orthogroups unique to that genome (Figure 5). Species on long 290 branches or located at the base of lin eages have the highest number of unique orthogroups, 291 with B. saltans lake Konstanz having t he most unique orthogroups. We next calculated t he 292 tot al number of orthogroups and the number of unique orthogroups in t hree kinetoplastid 293 genera, Bodo, Leis h mania , and Trypanosoma (Figure 5). Consistent with previous studies 294 showing gene number reduct ion in parasit ic lineages compared to B. salta ns (11), genus 295 Bodo has the greatest number o f orth ogroups, and the highest proportion of unique 296 ort hogroups, 53%, compared t o 35% and 5.5 % in Tr y panosoma and L e is hma n i a, 297 respectively. 298 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 15 Finally, we compared t he completeness of metabolic pathways present in B odo spp. 299 and kinetoplastid species available in the Kyot o Encyclopedia of Genes and Genomes (KEGG) 300 database (Figure 6). The four Bodo spp. are nearly identical in their metabo lic capacity, with 301 almost no KEGG modules mi ssing from one lineage and present in another (Figure 6). 302 Exceptions are Arginine biosynthesis and Ribof lavin biosynthesi s, whi ch are missing in B. 303 saltans lake Konstanz and partially present (25-80% complete) in the o ther species. The four 304 Bodo spp. also share most of their metabolic capacity w ith the parasitic lineages, as 305 previously reported (11). However, we do observe some pat hways that are complete or 306 nearly complete (>=80% complet e) in the Bod o lineage and absent or nearly absent ( 309 beta-alanin e, thymine => 3-am inoisobutanoate) . Conversely, the G lyoxylate cycle is partial 310 (60 % complete) in the Trypan osoma and L eishmania, and near absent in all Bodo spp., 311 consistent w ith previous findings in B. saltans (36), while Nucleotide sugar biosynthesis 312 (galactose => UDP- galactose) is partially complete (50 %) in all Leishmania and T. c r uz i, but 313 absent in all Bo d o spp. Figure 6 only shows modules that are at least 60% complete in at 314 least one or more genomes. The f ull list of all KEGG modules identified in each genome, 315 including those less than 50% i s in Additional File 10. 316 317 Discu ssion 318 In this study, we show that single-cell sequencing reveals unexpected levels of 319 diversity among uncultured candidat e Bodo spp. cells f rom a single enviro nmental sample. 320 Met abarcoding, part icularly of SSU rDNA, and met agenomics have become the norm for 321 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 16 assessing the diversity of microbial eukaryotes in environmen tal samples. For protists, 322 traditionally these methods involve clustering reads and/or amplicons into Operational 323 Taxonomic units (OTUs) based on similarity thresholds, usually at around 97-99% identity, or 324 into Amplicon Sequence Variants (ASVs) produced b y t he DADA2 pipeline (37). The seven 325 single cells presented here all have a pairw ise SSU rD N A identity greater th en 99% and 326 would therefore have b een clustered into one O TU in some metabarcoding pipelines. 327 Instead, by generating single-cell genomes, we show that these seven cells form th ree 328 distinct clades, with each clade potentially representing a novel species. ANI has been used 329 to compare distances within and between closel y-related species of eukaryotic microbes 330 including microsporidia, yeast and some protists, with an interspecies ANI cutoff of 95% 331 considered appropriate (38–40) . Here, we found that genomes within each of the three 332 putative B odo species have an ANI > 98%, supporting the conclusion that each clade 333 represent s a Bodo species. However, ANI is not useful when comparing more dissimilar 334 genomes, where AAI becomes a more useful met ric (41,42). While more commonly used to 335 compare prokaryotic genomes, AAI has also been applied to eukaryotes. A comparison of 336 1,196 Human and Mouse prot ein sequences show ed an AAI of 85% (43), while recent 337 studies in f ungal linages found AAI values for members of the same species were often 338 >97% and demonstrated an est imated a genus boundary threshold of ~70-75% AAI f or the 339 family Hypoxylaceae (Ascomycota) ( 44,45). W hile further w ork is needed to explore 340 appropriate species and genus cutoffs for AAI and ANI within and between protist lineages, 341 our results showing that AAI with in each Bodo spe cies w as > 97% suggest that the species 342 thresholds of 95- 98% used in bacteria and fungi may be appropriate here too. 343 All three novel Bo do species presented here are more similar by SSU rDN A identity 344 to B. edax ATCC strain 30903 than they are t o B. s alt ans lake Konstanz. However, as no 345 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 17 other molecular sequence dat a exists from t his organism, we are unable t o determine if any 346 or all are indeed B. ed ax . Callahan et al. (2002) mentions that B. edax contains 347 endosymbiotic bacteria and differs morphologically f rom B. saltans by lacking 348 mastigonemes on the anterior flagellum (5) , yet several authors have proposed that B. edax 349 be reclassified as an isolate of B. salta ns (2,10). Given the extensive genetic diversity 350 observed among the B odo spp. present ed in this paper, more work is needed to determine 351 if B. edax is a true species. Importantly , this study sugg ests that the SS U rDN A locus is 352 insufficient for species delimitation within genus Bodo as it does not reflect the genomic 353 diversity within the group. 354 Furt her support ing our conclusion t hat the three Bo do clades represent three 355 species i s the f inding that each putative species harbors a unique species of Holosporlaes 356 bacterium. The three novel bacterial species split into two genera, one species belonging to 357 genus “ Ca. Bodovickermani, and the other two species forming a novel genus. The 358 phylogeny of the host Bod o spp. appears congruent wit h the phylogeny of the 359 endosymbionts, however more sampling covering greater taxonomic dist ribution is needed 360 to test if this congruence indicates co -evolution, or is the result of oth er processes ( 46,47). 361 Infectivity has been demonstrated in the Ho losporaceae family of Holosporales, 362 endosymbionts of various Param ecium ciliates ( 16,48), and some Holosporales have 363 demonstrated the ab ility to invade n ovel hosts experimentally (49,50). However, the 364 distribut ion and genetic differences observed between the Holosporales genomes 365 presented here suggest that these endosymbionts are and have been associated with their 366 respective host s for a long period. 367 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 18 The nature of the symbioses between H olos por ales bacteria and their hosts is 368 complex and varied (13). Midha et al. ( 2021) hypothesi zed that B. saltans is dependent on 369 its endosymbiont due to the presence of three putative addiction p olymorphic toxin/anti-370 toxin systems encoded in the “ Ca. B. vickermanii” genome ( 18). These systems typically 371 encode a long-lived toxin molecule alongside an antitoxin with a shorter h alf-life, so that if 372 the system is lost, the toxin is activated and becomes lethal. We found evidence for a 373 putative toxin/anti- toxin system in the two genomes belonging to genus “ Ca . 374 Bodovickermani”, but not in five gen omes in the novel genus. However, th e identity and/or 375 functionality of th e toxin/anti-tox in systems described in this study is speculative, as the 376 proteins and loci differ from those in “ Ca. B. vickermanii” by encoding fewer orphan 377 modules. The lack of a similar system in the genomes of the n ovel Holosporales genus 378 presented in t his study suggests that alt ernative processes may be maintaing this symbiosis. 379 All kineto plasts including B. saltans lack the biosynthet ic pat hway required for heme 380 biosynthesis (36,51) , as do the Bodo spp. presented in this study. The Kinetoplastids 381 Angom o nas and Str igomonas harbour endosymbiotic beta-proteobacteria which provide 382 their hosts with nutrien ts including heme, essential ammino acids, and vitamins (35, 52,53). 383 We did f ind partially complete pathw ays for heme biosynthesi s ( >=70% complet e) in the 384 novel genus, and a partially complet e pathway for Lysine biosynthesi s (>=66% complete) in 385 all the B odo Holosporlaes genomes. Theref ore, the two H olosporales genera presented here 386 could demonstrate two alternative symbiotic strategies; with genus “ Ca. Bodovickermani” 387 becoming addictive to its host through the toxin/antitoxin system, while the novel genus is 388 possibly a source of heme and/or other essential nut rients. All seven Holosporales genomes 389 presented in this paper encode a Type VI secretion system, con served throughout 390 Holosporales and hypot hesized to play a role in host-endosymbiont int eractions (13). 391 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 19 392 Previous st udies have shown that B. s altans sha res the majority of its meta bolic 393 pathways w ith the parasitic lineages Trypanosoma and leis hma ni a ( 11,36). Our results show 394 that the metabolic capacit y of B. s alt ans is shared with the novel Bodo spp. presented here. 395 We also show t hat the Bod o genus contains an abundance of unexplored and uncategorized 396 protein d iversity. The lineage contains far more unique o rthogroups and proteins than the 397 parasitic lineages. However, as most of these protein sequences have no known homology 398 to functional do mains in databases such as UniProt and KEGG; the functional implications of 399 this unique gene repertoire remains unknown. 400 401

402 Most of our knowledge of protist genomic diversity and their symbionts comes from 403 studies of species that are culturable. However, we know that environmen tal samples 404 contain a wealth of underexplored d iversity. Thi s study uses single-cell sequencing and 405 comparat ive genomics to show that seven B odo spp. cells f rom a single environmen tal 406 sample represent three po tentially novel species, each harbouring a novel and unique 407 species of bacterial endosymbiont. C omparing these data, we demonstrate the varied 408 nature of these symbioses and show that single-cell sequencing i s a power ful meth od for 409 exploring the diversity of uncultured protists and their cobionts. 410 411 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 20

412 Sample collection 413 Surf ace water (~1 m) was collected fr om t he River Leam (52.287295, - 1.547563), 414 Royal Leamington Spa (UK) in August 2022. Init ially, t he sample was prefiltered t o remove 415 larger debris and then con cent rated using diff erent polycarbonat e filt er pore sizes 416 (Millipore) t o obtain concent rated subsamples of different protist size ranges (10 t o 40 µm 417 and 0.8 t o 5 µm). The subsamples were supplement ed with 2- 3 autoclaved barley grains to 418 support heterotrophic/mixotrophic growth via bact eria increase and incubat ed for one 419 week prior to cell sorting. Single-cell nuclei were stained with 1xSybr G reen for 10 minut es 420 and sorted int o 96- w ell microplates (pre-filled wit h 5ul autoclaved/sterile filtered media), 421 using fluorescence- activated cell sorting (FACS; flow rate=1) and selecting against 422 chlorophyll a- negative cells while selectiong for SybrGreen-positive cells. After cell sorting, 423 10 µl of RLT-plus lysi s buffer (Qiagen) was added to t he wells and the plat e was frozen at -424 80°C until further processing. 425 426 Whole genome amplification, library preparation, and sequencing 427 A modif ied G&T- seq protocol (54) was carried out as follows. U sing a magnet ic 428 separator, Dynabeads MyOne Streptavidin C1 (Invitrogen) beads were washed according to 429 the manufactur er’s guidance and then incubated with 2 × Binding &Wash buffer (10 mM 430 Tris-HCl pH 7.5, 1 mM EDTA, 2 M NaCl) and Biotinylated Oligo- dT primer (IDT, 5’-431 /BiotinTEG / AAG CAG TGG TAT CAA CG C AG A GTA CTT TTT TTT TTT TTT TTT TTT TTT TTT T TT 432 TVN-3’) at 100 μM for 30 minutes at room temperature on a rotator. The oligo-t reated 433 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 21 beads were washed four t imes in 1 × Binding & Wash buffer (5 mM Tr is- HCl pH 7.5, 0.5 mM 434 EDTA, 1 M NaCl) and then suspended in 1 × SuperScript II First Strand Buff er (Invit rogen) 435 supplemented with SUPERaseIn R N ase Inhibit or (Invitrogen) to a final concentration of 1 436 U/μl. The lysate w as thawed on ice. 10 μl of prepared oligo-dT beads was added to each well 437 containing 12 μl cell lysate using a Dragonf ly Discovery liquid dispenser ( SPT Labtech). The 438 lysat e plat e was sealed and incubated on a ThermoMixer C (Eppendorf) w it h a heated lid at 439 21°C for 20 minutes shaking at 1000 rpm. Using a Fluent 480 liquid handling robot ( Tecan) 440 and a Magnum FLX magnet ic separator (Alpaqua) , the lysat e super natant was transferr ed to 441 a new plate, and t he beads wer e was hed twice in a cust om wash buffer ( 50 mM Tr is-HCl pH 442 8.3, 75 mM KCl, 3 mM MgCl2, 10 mM DTT, 0. 5% Tween-20). The supernatant f rom the 443 washes was added t o the left-over cell lysate containing t he genomic DNA which w as stored 444 at -20°C overnight. The mRN A was not used for this study. The r emaining cell lysat e was 445 thawed and subjected to a 0.6 × vols Ampure XP clean- up with 80% et hanol. The bead-446 bound gD N A was i sothermally amplified for 3 hours at 30°C t hen 10 minut es at 65° C using a 447 miniaturised (1/5 vol s) Repli-g Single-Cell assay ( Q iagen). The amplified gDN A was cleaned 448 up with 0.8 × vols Ampure XP and 80% ethanol, then elu ted in 10 mM Tris-HCl. 449 Sequencing libraries for t his project were const ructed by t he Techni cal Genomics 450 Group at the Earlham Inst itute, N orwich, U K. Initial libraries w ere const ruct ed for shallow 451 depth sequencing as follows: gD NA was quantified by f luorescence ( Quant- iT HS-DNA, 452 Invitrogen) on an Infinite Pro 200 plate reader (Tecan) t hen normalised to a f ina l 453 concentrat ion of 0.2 ng/μl in 10 mM Tris-HCl. The Mosquito HV and D ragonfly Di scover y 454 liquid handling inst rument s (SPT Labtech) were used to prepare miniaturised (1/12.5 vols) 455 Next era XT ( Illumina) dual-indexed sequencing libraries as follows: A tagmentation 456 mastermix comprising two part s TD buffer (Illumina) and one part ATM (Illumina) was 457 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 22 prepared. Using the Dragonfly Discovery, 1.2 µl tagmentat ion mast ermix was dispensed per 458 well to a 384-well skirted PCR pl ate (Eppendorf) . U sing the Mosqui to HV, 0.4 µl of 459 normalised DN A was transferred to each well containing tagment ation mastermix. The plate 460 was sealed and spun down then incubated at 55°C for 10 minutes on a ther mal cycler. Using 461 the Dragonfly Discovery, 0.4 µl 0.2% SDS was dispensed to each reaction. The plat e was 462 sealed, spun down and then incubated at room temperature for 5 m inut es. Using the 463 Dragonfly Discovery 1. 2 µl NPM (Illumina) was added to each react ion. Using t he Mosquito 464 HV, 0.8 µl index primers pairs (i5 + i7) at a concentration of 0.5 µM was added to each 465 reaction ensuring a unique combination in each well. The libraries were amplified under the 466 following thermal cycler condit ions: 72° C for 3 minutes, 95°C for 30 seconds, 12 cycles (95°C 467 for 10 seconds, 55°C for 30 seconds, 72° C for 60 seconds), 72°C for 5 minutes, 4°C hold. The 468 librar ies w er e pooled and cleaned up using 0.8 × vols Ampure XP and 80% et hanol. The 469 librar y pools wer e eluted in 20 µl 10mM Tr is-HCl and assessed using a Bioanalyzer HS D N A 470 assay (Agilent) , HS DN A Q ubit assay (Invitrogen) and finally an Illumina Library 471 Quantification Kit assay (KAPA). These r eads f rom these libraries were used f or t axonomic 472 assessment of the cells. 473 Aft er t axonomic assessment deep sequencing libraries f or t he Bodo cells were 474 const ructed using t he KAPA High Throughout Library Prep Kit (Roche Part N o: 475 KK8234/07961901001). W her e possi ble, 1µg of genomic DNA was sheared t o 450bp usin g 476 the Covaris ML230 Sonicator ( Covaris) and the ends of t he DNA were repaired; 3' t o 5' 477 exonuclease activity removed t he 3' overhangs and the polymerase activity filled in the 5' 478 overhangs creating blunt ends. A single ‘A’ nucleot ide was added t o the 3’ ends of t he blunt 479 fragments to allow for the ligation of barcoded adapters ( 6bp - Perkin Elm er NEX TFLEX DNA 480 Barcodes 1-48 ( NOVA- 514101/2/3/4) ) or ( 12bp - Perkin Elmer N EXTFLEX-HT (NOVA-481 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 23 51474/ 5/6/7)) at a concent ration of 6µM prior to a 0. 8x clean up using Beckman Coulter 482 AMPure XP beads (A63882). The size of the libraries was estimated using an Agilent High 483 Sensit ivit y DNA chip (5067- 4626) and the concent rations w ere quant ified by fluorescence 484 with a High Sensitivity Q ubit assay ( ThermoFisher Q32854). 485 Bot h shallow and deep sequencing l ibraries wer e sequenced on either an Illumina 486 NovaSeq 6000 (cells A8, A10, B2, F10, G10, H10) with SP Reagent Kit v1.5 kit (300 cycles), , 487 using Real- tim e A nalysi s (RTA) ( version 3.4.4) and Control Software (version 1.7.5) or an 488 Illumina NovaSeq X Plus (cell B7) wit h 10B Reagent Kit (300 Cycle) using RTA (version 4.6.7) 489 and Cont rol Software (version 1.2.2) to produce 150 bp paired-end reads. The result ing BCL 490 files were converted to f astq with bcl2f astq ( version 2.0) . 491 492 Genome assembly, curation and classification 493 Sequencing reads were trimmed using Trim Galore (version 0.6.6) 494 (ht tps://github. com/FelixKrueger/Tri mGalore ) with Cutadapt (version 3.4) (55) . G enome 495 assemblies were generated using SPAdes (version 3.15. 5) (56) with single-cell mode enabled 496 (-- sc) and k-mer sizes 21, 33, 55, and 77. 497 Scaf folds less than 1,000 nt long were discarded from the assembly. Each assembly 498 was manually curated and cont aminant scaf folds/bins w ere removed using a combinat ion of 499 metagenomic binning with MetaBAT2 ( 57) based on t etra- nucleotide frequencies and 500 taxonomic classification with CAT (v5.2) (58), Blobtools (v 1.1.1) (59) and Tiara ( v1.0.1) (60) 501 and EukRep (v0.6.6)( 61). Bins that were classified as majority eukaryotic by all four 502 classifiers and majority “Euglenozoa” by Blobtools were retained, as well as any unbinned 503 scaff old classified as “Euglenozoa” by Blobtoo ls or CAT. 504 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 24 The Holosporales bins w ere first identif ied from the taxon omic classif ications 505 described above. However, all raw assembly bins were assessed for complet eness using 506 CheckM (v1.1.2) ( 62) , and all bins > 70% complete were furth er classified GTDB- TK ( v2.3.2) 507 on KBase. In all the cases, only one bin in each assembly was classified as 508 Alphapro teobacteria by GTDB-TK, and w as the bin retained as the Ho losporlaes 509 endosymbiont MAG. The retained scaff old set s were manually checked whi ch resulted in 510 the removal of one scaff old from t he B2 Holosporales bin that was a cont aminate from 511 another well on the 96 well plate. 512 The general statistics, classification summaries, CheckM results for all g enome bins, 513 and the GTDB- TK classifications are tabulated in Additional File 5. 514 The SSU rDNA BLASTs were done t hrough Geneious Prime. In each case, the subject 515 sequence with the highest bitscore was considered the top hit. 516 517 Gene prediction annotation 518 Eukaryotic gene prediction was done using Companion (v2.2.11) (63) at the 519 WebServer, with default sett ings and using Bodo sal t ans as the ref erence. rDN As were 520 annot at ed with Barrnap (v0.9). For one sample (B7) , t he SSU rDN A sequence was annotated 521 as t wo overlapping fragment s, which were manually merged into one sequence by 522 alignment in Geneious Prime. Four of t he seven SSU rDNA sequences are truncated (See 523 Table 1 and Additional File 1). Three of these, B2, G 10, and H10 w ere manually extended by 524 alignment of the scaffolds and truncated SSU rDNA sequences annotated b y Barrnap in 525 Geneious Prime and ext racting the longest region with shared homology from each scaff old. 526 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 25 All three sequences were extended by ~70 bp. For G10, the SSU rDNA annotation w as 527 truncated due to end of th e scaff old, w hile B2 and H10 are truncated due to loss of 528 homology, possibly due to miss-assembly of those scaf folds. Holosporales gene predictions 529 were done using PROKKA (v1.14. 6) ( 64). 530 Completeness of the genomes and protein sets was assesses with BU SCO ( v5.3.2) 531 (65). 532 533 Phylogenetic analyses of SSU rDNA sequences 534 SSU rDNA sequences were collect ed from GeneBank, filtering f or sequences that 535 were at least 1,000 nucleotides long ( Additional File 9). The sequences were aligned with 536 MAFFT (v7.520) using the FFT-NS-i model (66). The alignment was trimmed manually in 537 Geneious Prime to remove all positio ns with less than 50% coverage. The untrimm ed 538 alignment used to build the tree is in Additional file 1, while th e sequence accessions are in 539 Additional file 11. The Maximum-Likelihood (ML) tree was built with IQ-TREE (v2.3. 2) (67) 540 using the TIM3e+I+R3 model, which was t he best fit model determined by ModelFinder (68), 541 with 1,000 non- parametric bootstrap replicates, and rooting at two outgro up species 542 Dima s tigella tr ypa nif orm i s and Rhynchomonas nas ut a . The tree was visualized and plotted 543 with its associated distance matrix using Interactive Tree Of Life ( iTOL) (v7.2)(69) . 544 545 Phylogenomic analyses 546 For t he Kinetoplastea, prot ein sets for related species were collect ed from GeneBank 547 (Additional File 12) . O rthologous gene sets were identif ied using Orthofind er ( v2.5.4) (70). 548 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 26 Protein sequences for shared single copy orthogroups ( n = 488) were aligned w it h MAFFT 549 (v7.520) using the FFT- N S- 2 model (66), and trimmed with TrimAI (v2.0) using the gappyout 550 option ( 71). The trimmed alignments were concat enated with AMAS (72) int o a matrix with 551 234,727 amino acid sites. The ML tree w as produced f rom a partition ed analysis undertaken 552 using IQ-TREE ( v2. 3.2) w ith a partitioning scheme that merged the 488 pro teins into 20 553 part itions, wit h root ing at the outgroup Perkin sela sp.. Model selection was performed by 554 ModelFinder, and 1,000 non-parametric boot st raps w ere run. 555 For t he Holosporales, prot ein sets were collect ed from GeneBank for related species 556 and several closel y related MAG s (Addit ional File 10). For one M AG, G CA_002422845.1, a 557 protein set w as not available on G eneBank, so an annotation w as done using PRO K KA 558 (v1.14.6). Orthologous gene sets were identified using O rthof inder (v2.5.4). Protein 559 sequences from shared single copy orthogroups (n = 24) were aligned with MAFFT ( v7.520) 560 using the FFT- NS-2 model and t rimm ed wit h TrimAI (v2.0) using the gappyout option. The 561 trimm ed alignments were concatenated with AM AS into a matrix w ith 8,004 amino acid 562 sit es. The ML t ree was produced from a partitioned analysis undertaken using I Q-TREE 563 (v2.3.2) with a partitioning scheme that merged the 24 pro teins into 6 partitions, with 564 rooting at the out group containing Temper atibacter mar in us and Kor diimonas pumila . 565 Model selection was performed by ModelFinder, and 1000 non-parametric bootstraps were 566 run. 567 The Kinetoplastea tree (Figure 5) is th e Species Tree from All G enes (STAG) inferred 568 species tree produced by O rthofinder. It is inf erred from 331 gene trees from the 569 ort hogroups where all species are present. 570 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 27 All trees were visualized and plotted either using the Interactive Tree of Life (iTOL) or 571 in R using the package ggtree ( v3.14. 0). 572 573 Functional analyses 574 Pfams were assigned to each protein set using Int erProScan (v5. 52-86.0) using 575 option - aap Pfams and with the Pre- calculated match lookup service disabl ed (----disable-576 precalc). Pfam counts were tabulat ed in R. A PCA was run on the count matrix for Bod o spp. 577 after removal o f Pfam with no variation, which reduced the dataset f rom 2938 variables to 578 2545. The PCA was run in R using the function prcomp with scaling and centering. 579 KEGG identifiers were assigned t o each protein set using KofamKOALA or kofamscan (v1.3.0) 580 (73). The complet eness of each KEGG module was cal culat ed using kegg-pathway-581 completeness-tool (v1.3.0) ( https://gi thub.com/EBI- Metagenomics/kegg-pathways-582 completeness-tool ). 583 Venn diagrams showing the proportion of shared orthogroups were generat ed in R 584 using the package ggven (v0.1.10). For an orthogroup to be missing from a species, it must 585 be missing from all genomes in that species. For an orthogroup to be present, it can be 586 present in one or more genomes in that species. 587 The presence of bacterial secretion systems was assessed manually by visualizing the 588

of KofamKO ALA on the KEGG website using the KEGG mapper tool Reconstruct. The 589 number of unique components present in each protein set that mapped to each secretions 590 system was counted and compared t o the tot al number of unique component s list ed for 591 each system. The one exception was for the Type VI system, which was annot ated using a 592 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 28 combination of KEG G Reconstruct, BLASTp and psiblast, where the Type VI proteins from Ca . 593 B. vickermanii w ere used as queries to search a database of the protein sets from all seven 594 Holosporales MAG s. The query sequence identifiers, and the significant hit subject 595 identifiers are listed in Additional File 13. 596 The toxin anti/toxin systems were investigated manually using BLASTp and psiblast, 597 where t he proteins of Ca. B. vickermanii systems were used as queries to search a database 598 of the protein sets from all seven Holosporales MAGs. N terminal signal peptides were 599 identified using SignalP v 6. 0 at the webserver. 600 601 AAI, ANI, and AF analyses 602 Pairwise amino acid identity (AAI) was calculated f or the Bo d o spp. using the aai.rb 603 ruby script and BLAST+ v(2.16.0) that is part of the Enveomics Collection at 604 https: //github.com/lmrodriguezr/enveomics ( 74). For the Holosporales, pairwise AAI w as 605 calculated using FastAAI (34) . AN I and AF were calculated using the tool Skani (v0.2.2) (21). 606 607 Figures 608 All p lots were generated in R and edited for publication using Illustrator. 609 610 Declarations 611 Availability of data and materials 612 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 29 The raw reads and annotat ed assemblies f or Bodo and Bodo Holosporales have been 613 deposited in the European N ucleotid e Archive under accession PRJEB97217. The assembly 614 sequence f iles, protein sequence files and gff files used in this study have b een deposited in 615 Zenodo at 10.5281/zenodo.16948154. 616 617 Competing interests 618 The authors declare that there are no conflicts of interest. 619 620 Funding 621 This work was funded by the Wellcome Trust Darw in Tree of Life Awards (218328 622 and 226458), and by the Biotechnology and Biological Sciences Research Council (BBSRC), 623 part of UK Research and Innovat ion, through the Earlham Institute’s Core Capabilit y G rant s 624 (BB/CCG1720/1 and BB/CCG 2220/1), it s Strat egic Programme Grant Decoding Biodiversit y 625 (BBX011089/1) and its constituent Work Package 2 BBS/E/ER/230002B, its N ational 626 Capability BBS/E/T/000PR9816 (NC1 - Supporting EI’s I SP s and the UK Communit y with 627 Genomics and Single Cell Analysis), and Transf ormat ive Genomics, National Bioscience 628 Research I nfrastruct ure (BBS/E/ER/23N B0006). TAR is supported by a Royal Societ y 629 Universit y Resear ch Fellowship (URF/R/191005). Part of this work was delivered with 630 support for t he physical HPC inf rastructure and data center delivered via the N BI Comput ing 631 infrastructure f or Science (CiS) group. 632 633 Authors' contributions 634 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 30 Sally D. W arring, Conceptualization, Data curation, Formal analysi s, Investigation, 635 Methodology, Resources, Visualization, W riting | Jamie McG owan, Con ceptualization, Data 636 curation, Investigation, Methodology, Resources, Software, Writing |Estelle S . Kilias, 637 Invest igation, Met hodology, Resources, Writing | James Lipscombe, Invest igation, 638 Methodology, Resources, Writing | El isabet Alacid, Investigation, Methodology, Resources, 639 Writing | Tom Barker, Investigation, Methodology | Leah Catchpole, Investigation, 640 Methodology, Project administration, Writing | Seanna McTaggart, Funding acquisition, 641 Project administration, Writing | Karim Gharbi, Funding acquisition, Methodology, 642 Resources, Supervision, Writing | Thomas A. Richards, Conceptualization, Funding 643 acquisition, Resources, Supervision, W riting | David Swarbreck, Conceptualization, Funding 644 acquisition, Supervision | Neil Hall, Conceptualization, Funding acquisition, Resources, 645 Supervision, W riting. 646 647

648 We would like to acknowledge the members of t he Technical Genomics and Core 649 Bioinformatics groups at the Earlham Institute, and note the specific contributions of Chris 650 Watkins, Sacha Lucchini, Kendall Baker, and N eil Shearer. W e also acknowledge the work 651 delivered via the Laboratory Managers and Resear ch Computing Groups at EI who manage 652 and deliver High Performance Computing at EI. 653 654

655 1. Adl SM, S imp son A G B, Lan e CE, L uke š J, Bas s D, Bow se r SS, et al . The Rev i s e d Cla s s i f i catio n of 656 Eukary ot e s . J E ukaryot Mic ro biol. 20 12 S ep t ;5 9 (5) :429–514 . 657 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 31 2. Mor e ira D, Lóp e z - Ga rcía P, Vic kerman K . An updat ed view of kine t op la stid phylog eny using 658 environme n t a l s e qu ence s a nd a c lo se r o utgroup: p ropo sal for a n e w c las s i fic ati o n of th e clas s 659 Kinetopla s tea . I n t J Sys t Evol Microbiol . 2004 S ept 1;54(5 ):18 61 –75. 660 3. Simpson A GB, S tev en s JR, Luk e š J. The ev olution a nd dive rsity of k in etopl a stid fla g ella tes . Trend s 661 Para si t o l. 2006 Ap r ;2 2 (4):168 –74. 662 4. Simpson A GB, Luke š J , Roger A J. Th e Evo l u tionary Hi sto r y of Kine topl a stid s and Th eir 663 Kinetopla s t s . M o l Bi ol E vol. 2 002 Dec 1;1 9(12):20 71 –83. 664 5. Calla h an HA, Li t a k er RW, Noga E J. Molec ul ar Tax onomy of t he Su bord er Bod oni n a (Or d er 665 Kinetopla s tida ), I ncludi ng t he Impo rtan t Fi s h P ara s i t e , Ichthyobodo necator . J E uk aryot 666 M i crobiol. 2002 Ma r;49 (2):1 19 –28. 667 6. D e s c hamp s P , L ara E, Ma ra nde W , Lo pez-G arci a P , E kelund F, Mor eira D. Phy loge n omic Ana lysi s 668 of Kinet opla s ti ds Suppor t s That Tryp ano so matid s Ar o s e f rom wi t hi n Bodoni ds . Mo l Bi ol Ev ol. 669 2011 Jan 1;28(1 ):53– 8 . 670 7. Br oo ker BE. F ine st r uc ture o f Bodo sal tan s and Bo do ca uda tu s (Zooma stigo phora : Pr o tozoa ) an d 671 t he ir a f finiti e s with the Trypan o somati da e. Bull Br Mus Nat Hi st Z ool. 197 1;22 (3): 89–1 02. 672 8. Vick er m an K, P re st on T. C omp ara tive cell biology of the k ine topla s tid fl ag ella te s. B iol 673 Kinetopla s tida . 1 979;2:35 –130 . 674 9. D o leze l D, Jirků M, M asl ov DA , Luk e s J . P hyloge ny of the b odonid fla gella t e s (Ki n e t o pla st i da) 675 base d on small- subu nit r RNA gen e s eque nces . Int J S yst E vol Mic robiol. 2 000 Se p t 1;50(5) :1943–676 51. 677 10. Von D er H e yd en S , Ch ao EE, Vic kerman K , Cav alie r-Smi th T. Ribo somal R N A P hy logen y of 678 Bodonid a nd Diplon emi d Fla gella t e s an d the Ev olution o f Eug lenoz oa. J E ukaryot Mic r o biol. 20 04 679 J ul y ; 5 1 ( 4) : 4 02 – 1 6. 680 11. Jack s o n AP , O t to TD, As let t M, A r m stron g SD, B r ing aud F , Sc hlac h t A , e t al. Kin eto plas tid 681 Phylog enomic s Re veal s t he E volution a r y Innova tion s A s s ociat ed with th e Origin s of P a r a siti sm . 682 Cur r Bi ol. 2016 Jan ;26(2 ):161 –72 . 683 12. Szokoli F , Ca s telli M , S aba neyev a E, Schr a ll hammer M, Kren e k S, D o a k TG, e t a l. Di sen tangli ng 684 t he Tax onomy o f Rick ett sial e s a nd D es c ri ption o f Tw o Novel S ymbion ts (“C a ndida tus Beal eia 685 paramac r o n uclea r i s ” and “ Ca ndida tu s F o kinia cryptic a”) S haring t h e Cyt opla sm of the Cili at e 686 Pr o tis t P arame cium b iau relia . Drak e HL, e dito r . Appl Env iron Mic robiol . 2016 D e c 687 15;82(24) :72 36–47. 688 13. Ca stel li M, P e troni G. An Evolution a r y -F o cu s ed R eview of th e Holo spo rale s 689 (A l phap rot eobac t e r i a ) : Diver s i t y , Ho s t In terac tion s , an d Ta xonomi c Re-r ank ing a s Holos porin eae 690 Subord. Nov. Mi crob Ec ol . 202 5 Mar 14;8 8(1):1 5 . 691 14. Schral lhammer M , P o t e khin A . Epidemio l ogy o f N uc leu s -D w ellin g Hol os por a : In fec tion , 692 Trans m i ssion , Ada p t a t i on, and I nt e r a cti o n w ith Param ec ium. In : Kloc M, edito r. Sy mbiosi s : 693 Cellula r , Mol ecul ar, M edic al and Evolu tio n ary Asp ect s [In t e r n e t]. Cham: S p r in ger I ntern ation al 694 Publi shing; 2020 [c i ted 2025 Sep t 24] . p. 105– 35. (R e sult s a nd Probl ems in Cell Dif fer enti ation ; 695 vol. 6 9). Avai labl e f rom: ht tp s: / /li nk . s prin ge r . com/10 .1007 / 978 -3-030 -51849 -3 _4 696 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 32 15. Schral lhammer M , C a stell i M, Pe troni G. Phy logene t i c rela tio n s hip s am ong en do s ymbi ot i c R-697 body produc er : Bac te r ia pr ovidin g their hos t t h e kil l er tra it. Sys t Appl Mic r o bi ol. 201 8 698 M a y;41(3) :21 3–20 . 699 16. Fujishima M. Inf ection a nd Mai n t e nanc e of Holo s po ra Spec i e s in Pa ram ecium cau da t u m. In : 700 Fujishima M , e dito r. Endo symbi onts in P aramec ium [ Int ern et ]. B erlin, H ei delb erg : S pr i nger 701 Berlin Heid el be rg; 2009 [cit ed 2 025 June 1 0]. p. 20 1 –25. ( Mic r o biol ogy Monograp hs; vol . 12 ). 702 A v aila ble from: h ttp s : / /link. spring er.c om / 1 0.1007/978 -3 -54 0-92677 -1_ 8 703 17. Fujishima M, Kawa i M, Ya mamo to R. Paramecium caudatum ac qui re s he at - s hoc k re si st a nc e i n 704 cili ary move ment by infec tion with th e e nd onucle ar s y mbiotic bac te rium Holospora obtusa . 705 FEMS Mi crobiol L et t. 2005 Feb ;243 ( 1 ) : 10 1–5. 706 18. Midha S , Rigde n DJ , Si ozio s S, Hur st G DD , Jack s on AP . Bodo sal tan s (Kine topl a stida ) i s d e pend ent 707 on a nov el Pa r a cae di bac t e r- lik e endo sym bi ont th at po sse s se s mul t i ple pu tativ e to xin-an t i tox in 708 s y stem s. ISME J. 2021 June ;15 (6) :1680– 9 4. 709 19. Jamet A, Na s s i f X. N e w P la yer s in th e Tox in Fi eld: P ol ymorphic Toxi n System s in Ba c teria . mBio . 710 2015 July ;6(3) :e0028 5 -15. 711 20. Sche cke nbac h F, Wyl ez ic h C, My lnikov A P, Wei te re M, A rndt H . Molec ul ar C ompa r iso ns o f 712 Fre s hw at er and Ma r i ne Is olat e s of the Sa me Morpho s p ecie s of H et ero tr o phic F lag e llat es . Appl 713 Environ Mic robiol. 2 006 Oc t;72(10 ):6 638 –43. 714 21. Shaw J, Yu Y W . Fa st an d robust me tagen o mic sequen ce c ompari so n thro u gh spa r s e cha ini ng 715 with skani . N a t M e th od s . 2 023 Nov;20 ( 1 1):1661–5 . 716 22. Schn eider D, Züh lke D, Po e hlein A, Riedel K, Daniel R. M et agen ome -A ss embled Ge n ome 717 Sequen ce s from Dif f e r e n t W a st ewa te r Tre atmen t St age s i n Ge r ma ny. T hr a s h J C, e ditor . 718 M i crobiol Re sour Ann oun c. 20 21 July 8;1 0(27): e00 504-21 . 719 23. Muñoz-Góme z SA, He s s S, Bu rger G, L ang BF , Susko E, Sl amovit s CH, e t al . A n upda ted phylog eny 720 of th e Alphap rot eoba c teri a rev ea ls tha t t he para si tic Rick ettsi ale s a nd Hol o s p or a l es have 721 independ e nt o r i gins . eL if e. 201 9 Feb 25; 8:e425 35. 722 24. Park s DH, Rinke C , Chuvoc hina M, Chau mei l PA, W oodc r o ft B J, Ev an s P N, et al. R ecov ery of 723 nearly 8,000 me t a g enome -a ssembl ed ge nomes sub stan t ia lly expan d s the t re e o f life. Nat 724 M i crobiol. 2017 S ep t 11;2 (11):1 533 –42. 725 25. Liu L , Wang Y, Ya ng Y, W a ng D, Che ng SH , Z heng C, et al . Ch art i ng the c omplex ity of th e 726 activ ated sludg e microbiome th r oug h a h y brid seq uenc ing s tr a t e gy. Mic r obiome . 2 021 727 D e c;9 (1) :205 . 728 26. Park s DH, Chuvoc hina M , Wai t e DW, Ri n ke C, Ska rshew ski A, Chaume il PA, et a l . A s t a nda rdized 729 bacteri al taxonom y ba sed on ge nome ph ylog eny sub st a n tially revi se s t he t r e e o f l ife. N at 730 Biotechn ol . 2018 N ov ;36 (10):996 – 1004. 731 27. Prob st AJ, L add B , Ja ret t JK, Gelle r - Mc G r ath DE, Si eber CMK, Emer s on J B, e t al. Di ff ere ntia l d epth 732 dis t ribu t ion o f microbia l function an d putativ e s y mbiont s t hr ough sedim ent -ho ste d aqu if er s in 733 t he de ep te rre s t ria l s ub s urf ace . Nat Micr obiol. 2018 Ja n 29;3 (3) :328 –36 . 734 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 33 28. H e ss S, S uthaus A , Melko ni an M. “ Can did a t us F inniell a” (Ric ket t siale s, Alpha prot e oba cte ria) , 735 N ov el End o s ym biont s o f Vir i dirap tor i d A moebofl ag ella te s ( Ce rcozoa, R hiz a r i a ) . V oordouw G, 736 edito r. Appl En viron Mic robi ol . 2016 Jan 15;82 (2) :659–70. 737 29. Jain C , Ro dr i gue z-R LM, Phi llippy AM, Konsta nt i ni dis KT, Alu ru S . High th roughp ut ANI a n alysi s o f 738 90K proka r y otic g enome s r evea l s c lear s pec ie s bounda rie s. Nat Commun. 2018 N ov 739 30;9(1) :5114 . 740 30. G ori s J , Kons tan tinidi s KT, Klappe nbach JA, Co enye T, Van da mme P, Tie dje J M . DNA -DNA 741 hybridizatio n v alue s and t h e ir re lati on shi p to who le -genome s e q uence s im ilari tie s . Int J Sy s t Ev ol 742 M i crobiol. 2007 J an; 57 (P t 1):81 –91 . 743 31. Konsta nt i ni dis KT, T iedje J M. T oward s a Genome -B a s ed Taxon omy for Proka r y ote s. J Ba c teri ol. 744 2005 S ept 15;187 (18 ):6258–6 4 . 745 32. Rodriguez -R LM, K on s t antini di s KT. Bypa ssing Cul tivation To Id enti fy Bact eria l Sp e c ies : Cult ure -746 independ e nt gen omi c app r oa ch e s iden tif y credibl y distinc t c lus t e rs , av oid c ultivati on bias , and 747 pr ov ide true in s ig hts int o mi crobial s pe ci es . Mic r ob e Mag. 2 014 Mar 1 ;9(3 ):1 11–8 . 748 33. Park MJ , Kim YJ, Park M , Yu J , Namirimu T, Roh Y R, et al . Establ i s hm en t of Gen ome Ba sed 749 Cr i te r ia for Cla s s if ic ation o f th e Famil y D esul fovibri ona cea e and P r opo s al of Two Nov el Ge ner a, 750 A l kali de s ul f o vibrio g en . nov. and S alid e s ulfov ibr i o gen . nov. F r o n t Mi crobi ol. 2 02 2 May 751 25;13:73 8205. 752 34. G e r h a r dt K, Rui z - Pe r e z C A, Ro d r ig uez - R LM, Jai n C, Tie dje JM , Col e JR, et al . Fa stA AI: e ffi cient 753 estima t i on of gen o me av era ge a mino ac i d ide nt i t y and phylum -le vel r ela t i on s hi p s using 754 t e trame rs o f univ er sal prote in s . N uc leic Acids R e s. 20 25 Apr 22;53 (8) :g kaf348 . 755 35. H a r me r J, Yurchenk o V, N enarokov a A , L ukeš J , Ginge r ML . Fa r mi ng, slavi ng and e nslave m ent: 756 his t orie s o f end o s ym biose s during kine to pl astid evolu tion . Para si tology . 2018 757 Sept;145 (10 ):1 311–23 . 758 36. O ppe r d oe s F R, Butenko A , Fleg on tov P, Y u rchenk o V , Lukeš J. C ompara tiv e Me tab olism o f 759 Free -liv ing Bodo saltans a nd P ara si tic T r y pano somatid s . J Euk aryo t Microbi ol . 2016 760 Sept;63 (5) :657 –78. 761 37. Bur k i F, Sandi n MM , Jamy M. D iv er s ity a nd ec ology of pro ti s t s r evea led by met ab arc oding. C u r r 762 Biol. 20 21 Oct ;31(19) :R126 7–80 . 763 38. Majda S , Boeni gk J , Be i sser D. Intra s pe ci f ic Varia tion i n P rot i st s: Clue s fo r Micro evo lu tion from 764 Pot er i o s pum e lla l acu st r i s (Ch rys ophyc ea e). Ar c hibald J , edi tor . Gen ome Bi ol Evo l. 2019 Sept 765 1;11(9) :2492 –504. 766 39. Lac hanc e MA, Le e DK, H s i ang T . De linea ti ng yea s t sp eci e s with ge nome av e rage n ucleotid e 767 identi ty: a cal ibra tion o f AN I with h aplo n tic , het ero thallic Me t s c hniko wia spec ie s . Antoni e Va n 768 Leeuwenho ek . 2020 Dec;113 ( 1 2) :2097–1 06. 769 40. D e Albuq uerq ue N R M, Haa g KL. Using averag e nuc le otid e id en tity ( A NI ) to eva lu ate 770 microspo ridia spec i e s bounda rie s ba sed on t h eir ge netic r ela tedn e ss . J Euk aryot Microbio l. 2023 771 M a r ; 70 (2): e12944 . 772 41. Konsta nt i ni dis KT, Ti edje J M. P rok aryo tic ta xonomy and phy logen y in the genomic era: 773 adva ncemen ts an d chal lenge s ah e ad . Curr Opin Micr obio l. 2007 O c t ; 10(5 ):504 –9. 774 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 34 42. Q in Q L , Xie B B, Zhang XY, Ch en XL, Zhou BC , Zhou J, e t al . A P ropo sed Genu s B oun da r y for t h e 775 Pr o karyo te s Ba s ed on G enom ic In sigh ts . J Bac t eriol. 201 4 J une 15 ;196 (12):221 0 –5 . 776 43. Makało wski W , Zha ng J, Bog us k i MS. Co m parative an aly s i s o f 119 6 ort h ologous mou se and 777 human full -lengt h mRNA and p r ot ei n s e q uence s. Genome R e s. 1996 Sep t;6 (9):846 – 57. 778 44. W i bberg D , Rupp O , Blo m J, J elon ek L, Kr öb er M , Ve rwaai jen B, e t a l. Deve lo pmen t of a 779 Rhiz oc toni a s olan i A G 1 -I B Spec i fic G e n e Model E nabl es Compa r a t iv e Genom e A n a lyse s be twee n 780 Phytopa thog enic R . sol ani A G1- IA, AG1 -I B, AG3 and AG8 I sol at e s. Yu J H, e dito r . P LOS ONE. 2015 781 D e c 21;1 0(12 ): e0144 769. 782 45. W i bberg D , Stad l er M, La mb ert C , Bunk B, Spröer C, Rüc ke rt C , e t a l . High q uality g enome 783 s eq uence s o f thi rte en Hypo xyla cea e ( A s c omyc ot a ) s tr ength en t h e phylog ene tic fa mily back bone 784 and ena bl e the d i s c ove ry o f new ta xa . F ungal D i ver s . 2021 Ja n; 106(1 ):7 –28. 785 46. G rou s s i n M, M azel F, Alm EJ. C o- ev olutio n a nd Co - s pe cia tio n o f Ho s t -G ut Bac teri a S ystem s . Cell 786 H os t Microb e. 2020 Jul y;28 (1):12 – 22. 787 47. Russo L , Mill er A D, Too ker J, Bjorn s t a d ON, S hea K. Q u an t i tativ e ev olution ary p at tern s in 788 biparti te ne twork s: Vic arianc e, phy logen etic tr ack ing or dif fuse c o -evolu ti on? M e thod s Ecol Evol. 789 2018 Mar;9(3 ):76 1 –72. 790 48. Schral lhammer M , P o t e khin A . Epidemio l ogy o f N uc leu s -D w ellin g Hol os por a : In fec tion , 791 Trans m i ssion , Ada p t a t i on, and I nt e r a cti o n w ith Param ec ium. In : Kloc M, edito r. Sy mbiosi s : 792 Cellula r , Mol ecul ar, M edic al and Evolu tio n ary Asp ect s [In t e r n e t]. Cham: S p r in ger I ntern ation al 793 Publi shing; 2020 [c i ted 2025 Jun e 10] . p. 105–35 . (R e sult s a nd Probl ems in Ce ll Dif fer enti ation ; 794 vol. 6 9). Avai labl e f rom: ht tp s: / /li nk . s prin ge r . com/10 .1007 / 978 -3-030 -51849 -3 _4 795 49. Frelie r PF , Lo y JK, Kruppe nb ach B . Tr an sm ission o f Nec ro t izing H e pat opanc rea ti t i s in P enae u s 796 vanna mei. J I nv erte br Path ol. 19 93 Jan ;6 1(1):4 4 –8. 797 50. Schu lz F , La gkouv ardos I, Wa sch er F, Ai st lei t ne r K , Kos t a nj š e k R, Horn M. L if e in a n unusual 798 intr a cell ular nic he : a b act erial s y mbiont i nfec t i ng the nuc le us o f am o ebae . I SME J . 201 4 Aug 799 1;8(8) :1634– 44. 800 51. Kořený L , Lukeš J, O b or ní k M . Ev olu tion o f the h aem synthe tic pathw ay in kinet opl a stid 801 flagel la te s: An e ssen tia l pa thwa y tha t is not e s s e n t i al a ft er al l? Int J Pa ra s itol . 2010 802 Feb;40 ( 2 ) :1 49 –56. 803 52. Alves JMP , Vo egtly L, M a tveye v A V, L ara AM, D a Si lva F M, Ser rano M G, e t a l. I den tific a tion a nd 804 Phylog ene tic Analy si s of H eme S ynthe si s Gene s in Trypan os omatid s and Th eir Bac terial 805 Endosymbion ts . Lo pez-G arcia P, ed i t or . P LoS O NE. 20 11 Aug 10 ;6(8 ): e23518 . 806 53. Klein C C , Alve s JM P, S e r ran o M G , Buck GA, Va s c onc el o s ATR, Sag o t M F , e t al. Bi o sy nthe s i s of 807 Vitamin s an d Co fa ct or s in Bact erium -Har bouring Trypano soma tid s Dep end s on th e Symbi ot i c 808 A s socia t i on a s R evea led by G enomic Ana lyse s. Parki n s on J , edi tor . PLoS O N E. 201 3 Nov 809 19;8(11) :e79 786. 810 54. Macau lay I C , Hae rty W, Kumar P, Li Y I, Hu TX, Teng MJ, e t al. G & T - s eq : par allel s e qu encin g of 811 s ing le -ce ll ge nom e s and tr an script ome s. Nat M etho ds. 2015 June ;1 2(6 ):519–22 . 812 55. Mart in M . Cutada pt r emov es adap ter se qu enc e s from h igh- throu ghput sequenc i ng read s. 813 EMBnet .journal . 2011 May 2;17(1 ):10 –2 . 814 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 35 56. Banke vich A, Nur k S, A ntipov D , Gure vic h AA, Dvorkin M, Kulikov AS, e t al. S P Ad e s : A N ew 815 G e nome A s s embl y Algorit hm a nd It s A p p lica t io n s to S ingl e -Ce ll S e quenc ing. J Com pu t Biol. 201 2 816 M a y;19(5) :45 5–77 . 817 57. Kang DD, L i F, Kirt o n E, Tho ma s A, Egan R, An H, et al . Me taB AT 2: a n ada p t i ve bi n ning a lgorithm 818 for r obust and e f f i cien t g enome r ec on s t r u ction f r om metage nom e a s s em bli e s. P e erJ . 819 2019;7:e 7359. 820 58. von Mei jenf eldt FA B, Ar k hip ova K, Camb uy DD, C outinh o F H , D u tilh BE. Ro bu st t a xonom ic 821 clas si fica t io n of unch a rted microbia l seq uenc e s a nd bin s wi t h C AT and BA T . Gen o me Biol . 2019 822 O c t 22;20(1 ):217 . 823 59. Lae tsch D R, Blaxter ML . BlobToo l s: In te r r oga t io n o f geno me a s sem blie s. F10 00R es earch. 2017 824 J uly 31;6 :1287. 825 60. Karlick i M, Anton ow icz S, K arnkowska A . Tiara: de ep le a rning-b a sed c la s s i fi cation sys tem for 826 eukary ot i c sequ enc e s . B ioin forma Oxf E n gl. 20 22 Jan 3 ;38(2 ):344–5 0 . 827 61. W e s t PT, Pro bst A J, Grigo riev I V, Thoma s BC, Ban field JF . Genom e- recon structio n for eu k aryot es 828 fr om comple x natural microbia l comm u n itie s. G e nom e Re s. 201 8 Ap r;28(4 ):569 –8 0. 829 62. Park s DH, Imel for t M, Ske nn erton CT, Hugen holtz P, Tyson GW. Chec kM : a sse s sin g the q u ality o f 830 microbia l geno me s r ec overed f rom i s ola te s, sing l e ce ll s, and me t a genom es . G eno me Re s . 2015 831 J uly;25(7 ):10 43–55 . 832 63. H a e s e - Hil l W, Crouch K , Otto T D. Ann ota tion and v isua liza tion o f para si te , fung i a nd art hr opod 833 genome s w ith Compa nio n. Nuclei c Aci d s Re s. 2 024 July 5;52( W1): W39–44 . 834 64. Seema nn T. P r ok ka: r api d proka ryotic ge no me anno tatio n . Bioin fo rmatic s . 2 014 J uly 835 15;30(14) :20 68–9. 836 65. Manni M, Be r k eley MR , Seppe y M, S imã o F A, Zdobnov EM. BUS CO Upd at e: Nove l and 837 Str e amline d Wo r k flow s along with B road e r and D eep er P hylog e ne tic Cov er age f o r S coring of 838 Eukary ot i c, Pr oka r y otic, a nd Vi ral G enom e s. K elley J, edito r. Mol Bi o l Evol . 2021 Se pt 839 27;38(10) :46 47–54. 840 66. Katoh K, S tandl ey DM. MA F FT multipl e s equ enc e al ignm ent so ftwar e v er s ion 7: i mprove ment s 841 in perf ormanc e a nd us abi lity. M ol Bi ol Ev ol. 2013 Apr;30 (4) :772–80 . 842 67. Minh B Q, Sc hmidt HA, Che rnomor O, Sc h r e mpf D , W o o dham s MD, von Hae sel e r A , et al. I Q-843 TREE 2 : N e w Model s an d Ef fici en t Me tho ds f or Ph ylog e net ic Inf e r e nc e in t h e Gen omic Era. Mol 844 Biol Ev ol. 202 0 May 1;37(5) :1530 –4. 845 68. Kalya anamoo rthy S, Minh BQ, W ong TKF, v on Haes eler A , Je rmii n LS . Model Fi nder : fa st mod el 846 s el ection for a ccurat e phylog ene tic e s tim ate s. Nat Meth od s. 2017 Ju n e;14 (6):587 – 9. 847 69. Letunic I, Bo rk P. I nte rac t iv e Tr ee o f Life (i TO L ) v6: r ec en t upda te s to t he p hyloge n etic tr ee 848 dis p lay a nd a nnot ation t o ol . Nucle ic Ac id s Re s. 2024 July 5;52(W1 ): W78–82 . 849 70. Emms DM , Kelly S. O rt h o Finde r: phy loge ne t i c or thology inf erenc e fo r compa r a t i ve g enomic s. 850 G e nome Biol. 2019 D ec ;20(1 ):238 . 851 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 36 71. Capell a- G ut ié rrez S, S illa -M artí nez J M, G abald ón T. t rimAl: a t o ol for au tomat ed a lignme nt 852 t rimmi ng in l arge -sc ale phylog en etic anal y s e s. Bi o inf or m atic s . 2009 Aug 1;25 ( 1 5) : 1972–3. 853 72. Bor o wiec ML. AM AS: a fa st t oo l for alig n ment manip ula tion a nd compu t in g of su mmary 854 s t ati sti cs . Pe er J. 201 6 Jan 28 ;4:e1660 . 855 73. Ar a maki T, Blanc -Ma thieu R, End o H, Oh kubo K , Kaneh isa M, Go t o S, et a l . Kofam KOALA: KE G G 856 O rth olog a s signm e nt ba s ed on pro file HMM a nd adap tiv e score thre sh old. Val en cia A, edit or. 857 Bioinfo rmatic s . 20 20 Apr 1;36(7 ):22 51 –2. 858 74. Rodriguez -R LM, K on s t antini di s KT. The e nveo mics col lec tion: a toolbox f or s pe cial iz ed analy s e s 859 of microbia l ge n ome s and met ag enome s [In tern et ]. Pe erJ Prep rint s; 2016 [c it ed 2025 Apr 1]. 860 A v aila ble from: h ttp s : / /pe erj.com / p rep ri nts/1900 v1 861 862 863 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 37 Figure Legends 864 Figure 1 – SSU rDNA phylo gen y of Bodo spp. A. Pa irwise ANI of SSU rDNA sequences from 865 the seven single cell genomes B. ML phylogeny of SSU rDNA sequences from Bodo spp. and 866 related genera. The tree scale is substitutions per site. The seven single cell genomes are 867 shaded yellow, blue and pink to distin guish the three clades they form. B. s altans lake 868 Konstanz is shaded purple. The tw o columns on the right side of the tree show the pairwise 869 ANI between each taxon and the SSU rDNA from B. s alt ans lake Konstanz and B. edax . 870 871 Figure 2 – Genomic diversity of Bodo spp. A. ML phylogeny of 488 single copy ort hologous 872 proteins. Tree scale is substitutions per site B . Heatmap showing pairwise average AAI 873 values for all Bodo genomes C. Heatmap showing pairwi se AN I (lower triangle) and AF 874 (upper triangle) D. N umber and proport ion of shared ort hogroups for each Bodo spp. E. Bar 875 plot showing the proportion of prot ein annotat ions in each genome that have a Pfam 876 annot at ion F. PCA plot generated fro m the tally of Pfam domain p resent in each genome. In 877 all plots the th ree B odo spp. from sin gle cell genomes are shaded yellow, blue and pink to 878 distinguish the three clades they form while B. saltans lake Konst anz i s shaded purple. 879 880 Figure 3 – Genomic diversity of Holosporales endosymbionts. A. ML phyl ogeny of 24 single 881 copy orthologous proteins. Tree scale is substitutions per site B. Heatmap showing pairwise 882 ANI (lower triangle) and AF (upper triangle) for the Holosporlaes associated with B odo single 883 cells C. Number and proportion of shared orthogroups for each Holosporales species. In all 884 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 38 plot the Holosporales from the seven single cell genomes are shaded yellow, blue and pink 885 to distinguish the three clades they form while “ Ca. B. vickermanii” is shaded purple. 886 887 Figure 4 – Functional diversity among Holosporales endosymbionts. Heatmap in blue 888 showing t he completeness values f or KEGG modules found in each Bodo Holosporales 889 genome and Ca. B. vickermanii genome. Only modules that are at least 50% complete in one 890 or more genome are shown. Heat map in red shows t he number of prot eins belonging to 891 each of the bacterial secretion systems found in the genomes. The block t o the right o f the 892 heatmap shows the tot al number of proteins belonging to each secret ion system, as li sted 893 on the KEGG Brite database. 894 895 Figure 5 – Genomic uniqueness in Genus Bodo. S pecies tree infe rred by STAG. Tree scale is 896 substitutions per site. The support values are the number o f individual gen e trees that 897 contain that bipartit ion. The bar plots show t he tot al number of orthogroups and t he 898 number of unique orthogroups in each genome protein set . Pie charts are shaded t o show 899 the proport ion of orthogroups unique t o the three genera Bodo, Tr y pano s oma and 900 Lei shmania . 901 902 Figure 6 – Functional comparison of Kinetoplastida. Heat map showing t he completeness 903 values for KEGG modules in each genome. Only modules that are at least 60% complete in 904 one or more genomes are shown. 905 906 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint 39 Additional files 907 Additional File 1 - Bodo and related spp. SSU rDNA alignment. FASTA f ile cont aining t he 908 untrimm ed SSU rDNA alignment used to build t he tree in Figure 1. 909 Additional File 2 - AAI summary statistics for Bodo spp. Excel sheet listing the summary 910 metrics out put from the aai.rb script. Column “no. proteins used” shows t hat number of 911 pairwise comparisons used to calculate the average AAI. 912 Additional File 3 – Bodo skani AF full matrix. TSV file showing the AF calcul ated by skani in 913 both directions for each pairwise comparison. 914 Additional File 4 - Pfam counts and PCA loadings. Excell workbook wit h two sheet s. The 915 first sheet shows the Pfam frequencies in each B odo genome. The second sheet shows t he 916 loadings given to each Pfam for PC1-8. 917 Additional File 5 - Bodo SC raw genome assemblies statistics per bin. Excell sheet 918 tabulating statistics for all the raw single-cell genome assemblies, split by MetaBat2 bins. 919 Additional File 6 - FastAAI distance matrix. Excell sheet containing the A AI values f or the 920 Holsoporales heatmap in Figure 3A. 921 Additional File 7 – Holosporlaes skani AF full matrix. TSV file showing the AF calculated by 922 skani in both directions for each pairwise comparison. 923 Additional File 8 - Holosporales all KEGG completeness. Excell sheet showing the 924 completeness values for KEGG modul es ident if ied in each Holosporales protein set. 925 Additional File 9 - Toxin antitoxins in Bodo Holosporlaes B7 F10 figure. P DF do c um e n t 926 showing a: A. schematic of the putative toxin/antitoxin loci in B7 and F10 B. alignment 927 between the large toxin pro tein and alternative toxin, showing region of homology C. 928 SingnalP plots of N termini of the p utative multi-domain toxins from B7 and F10. 929 Additional File 10 - Kinetoplastids all KEGG completeness. Excell sheet showing the 930 completeness values for KEGG modul es ident if ied in each Bodo and Kinet oplast ea protein 931 set. 932 Additional File 11 - SSU rDNA sequences u sed in ph ylogeny . Excell sheet with N CBI 933 accession identifiers for all SSU rDNA sequences used to construct the phylogeny in Figure 1. 934 Additional File 12- all NCBI protein sets . Excell workbook with two sheets listing the N CBI 935 accession identifiers for all pro tein sets used in this study. Kinetoplastea sets are on sheet 1 936 and Holosporales sets are on sheet two. 937 Additional File 13 - TypeVI secretion system protein IDs . Excell sheet listing the Type VI 938 secretion system protein id entifiers used as queries from Ca. B. vickermanii, and the protein 939 identifiers f rom the Holosporales genomes that were signif icant hits. 940 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint A B Bootstrap 100 80-99 50-80 B. saltans AY490226.1 Bodo H10 Bodo G10 Bodo B2 Bodo A8 Bodo A10 B. edax AY028451.1 B. saltans JC02 AY490227.1 B. saltans JF693632.1 B. uncinatus AF208884.1 B. saltans strain Petersburg AF208887.1 B. saltans strain HFCC309 DQ207572.1 Bodo B7 Bodo F10 B. saltans AY490224.1 B. saltans isolate NG MF962814.1 B. saltans JC18 AY490229.1 B. saltans PP5 AY490223.1 B. saltans JC03 AY490228.1 B. saltans AY490222.1 B. saltans AY490232.1 B. saltans AY490230. B. saltans AY490231.1 B. sp. TGS2 AB585965.1 B. saltans strain Konstanz AF208889.1 B. sp. TGKH8 LC000678.1 B. saltans AY028452.1 B. saltans strain HFCC12 DQ207569.1 B. saltans AY490234.1 B. saltans strain SCCAP BS364 AY998648.1 B. sp. ATCC 50149 AY028449.1 B. saltans strain HFCC311 DQ207574.1 B. saltans strain HFCC14 DQ207571.1 B. saltans MH614643.1 B. saltans strain IOW92 KX431511.1 B. saltans AY490233.1 B. saltans strain HFCC310 DQ207573.1 B. saltans strain HFCC323 DQ207575.1 B. saltans strain HFCC13 DQ207570.1 B. caudatus AY028450.1 B. curvifilus AY425015.1 B. sp. isolate COLPROT774 MW355419.1 B. caudatus AY490218.1 B. caudatus strain SCCAP BC330 AY998649.1 B. caudatus AY490215.1 B. sorokini AF208888.1 B. sorokini strain ATCC 50641 AY425018.1 N. designis strain SCCAP BD54 AY998650.1 N. designis strain SCCAP BD55 AY998651.1 N. designis strain SCCAP BD56 AY998652.1 N. designis strain SCCAP BD57 AY998653.1 B. saliens AF174379.1 B. designis AY425016.1 B. designis AF209856.1 N. designis strain SCCAP BD52 AY998646.1 B. rostratus AY425017.1 N. designis strain SCCAP BD50 AY998643.1 N. designis strain SCCAP BD23 AY998644.1 N. designis strain SCCAP BD51 AY998645.1 B. designis AY490235.1 N. designis strain SCCAP BD53 AY998647.1 B. designis strain DH AF464896.1 B. celer AY490221.1 D. trypaniformis strain SCCAP DIM74 AY998641.1 R. nasuta strain SCCAP RH3 AY998642.1 89.6 91.0 82.5 83.4 76.4 76.0 81.6 81.7 94.5 95.0 95.3 0.0 82.5 77.6 83.8 83.6 97.6 95.9 82.3 82.7 96.7 89.5 94.6 95.9 99.9 95.4 82.8 83.6 97.8 94.7 81.7 81.7 80.7 81.4 99.7 88.5 81.4 82.3 97.6 96.1 96.2 96.9 82.3 82.4 80.5 79.6 89.5 89.0 98.4 95.0 94.2 94.8 79.0 79.1 95.1 100 95.3 99.6 95.6 98.2 95.5 96.8 84.6 84.4 97.6 96.0 83.3 84.1 93.5 94.8 94.8 96.0 94.7 95.8 82.4 83.1 82.0 82.7 81.3 81.5 87.0 87.6 0.0 95.3 96.7 94.6 77.5 77.4 85.9 85.9 95.7 98.0 82.0 82.9 95.2 99.6 94.9 99.0 99.9 95.4 98.7 94.6 95.3 99.6 95.2 99.6 78.4 78.5 93.9 94.4 94.7 99.4 98.4 95.0 95.3 99.9 82.8 83.8 94.7 98.0 96.2 96.8 93.4 96.0 87.4 88.0 82.0 82.3 94.2 94.9 B. saltans strain Konstanz B. edax Tree scale: 0.1 Percent identity 75 77 80 82 85 87 90 92 95 97 100 99.82 98.05 98.05 97.68 98.12 97.94 98.05 98.05 97.63 97.99 97.94 99.84 99.93 99.93 99.84 99.93 99.93 100 99.92 99.92 95 96 97 98 99 100 Bodo F10 Bodo B7 Bodo A10 Bodo A8 Bodo G10 Bodo B2 Bodo H10 100 100 100 100 100 100 100 100 F10 B7 A10 A8 G10 B2 H10 18s rDNA pairwise nucleotide identity % Figure 1 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint A C D E F B Bootstrap 100 Bodo H10 Bodo G10 Bodo B2 Bodo A8 Bodo A10 Bodo B7 Bodo F10 Bodo saltans Perkinsela 51.3 62.2 62.2 98.8 100 51.5 62.6 62.6 85.0 85.0 97.9 98.7 100 51.4 62.3 62.1 100 51.6 62.5 62.5 84.9 85.0 100 51.7 100 51.5 98.0 100 51.5 62.4 62.5 84.8 84.8 97.7 100 100 Tree scale: 0.1 Average pairwise amino acid identity (AAI) % H10 G10 B2 A8 A10 B7 F10 B. saltans 50 55 60 65 70 75 80 85 90 95 100 0.00 0.00 0.00 0.00 0.00 0.00 0.00 98.39 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 99.24 84.17 84.05 84.16 84.24 84.14 84.09 98.42 98.23 99.02 0 10 20 30 40 50 60 70 80 90 100 B.sal. B.sal. F10 B7 A10 A8 B2 H10 G10 0.00 0.00 81.56 0.00 0.00 0.00 0.00 0.00 0.00 83.12 0.00 0.00 0.00 15.63 15.46 0.00 0.00 0.00 15.30 15.15 80.14 0.00 0.00 0.00 14.70 14.53 76.72 82.20 0 10 20 30 40 50 60 70 80 90 100 F10 B7 A10 A8 B2 H10 G10 Aligned Fractiom (AF) % Average nucleotide Identity (ANI) % B. saltans B7/F10 A8/A10 B2/ G10/H10 806 (5.8) 571(4.1) 199(1.4) 570(4.1)257(1.9) 11(0.8) 898 (6.5) 77(0.6) 181(1.3) 138(1.0) 280(2.0) 1076(7.8) 530(3.8) 467(3.4) 7663(55.5) Shared orthogroups H10 G10 B2 A10 A8 F10 B7 B. sal 0.00 0.50 1.00 Proportion proteins with Pfam annotations Other proteinsProteins with Pfam annotation B. saltans A8A10 B7 F10 B2 G10H10−0.25 0.00 0.25 0.50 −0.9 −0.6 −0.3 0.0 PC1 (26.38%) PC2 (20.4%) Principle components analysis base on Pfam frequencies Figure 2 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint B C A A10 A8 B2 MAG GCA_018061565 B7 F10 Ca. B. vickermanii MAG GCA_002422845 Ca. C. indipagum Ca. C. primus Ca. C. primus 43202 Ca. N. abundans Ca. N. abundans 44043 Ca. G. agglomerans Ca. H. endosymbioticus H. curviuscula H. elegans E1 H. undulata HU1 H. obtusa F1 Ca. H. penaei Holosp. bac. Namur Ca. B. paramacronuclearis 11III1 Ca. B. paramacronuclearis 15I1 Ca. N. amoebiphila Ca. P. acanthamoebae C. varicaedens Ca. O.thessalonicensis L13 Ca. O. acanthamoebae Ca. P.symbiosus C. Finniella inopinata MAG GCA_002791835 MAG GCA_003542655 G10 H10 MAG GCA_024236035 K. pumila T. marinus eoAcanthamoeba sp. UWC8 A10 A8 B2 MAG GCA_018061565 B7 F10 Ca. B. vickermanii MAG GCA_002422845 Ca. C. indipagum Ca. C. primus Ca. C. primus 43202 Ca. N. abundans Ca. N. abundans 44043 Ca. G. agglomerans Ca. H. endosymbioticus H. curviuscula H. elegans E1 H. undulata HU1 H. obtusa F1 Ca. H. penaei Holosp. bac. Namur Ca. B. paramacronuclearis 11III1 Ca. B. paramacronuclearis 15I1 Ca. N. amoebiphila Ca. P. acanthamoebae C. varicaedens Ca. O.thessalonicensis L13 Ca. O. acanthamoebae Ca. P. symbiosus C. Finniella inopinata MAG GCA_002791835 MAG GCA_003542655 G10 H10 MAG GCA_024236035 K. pumila T. marinus eoAcanthamoeba sp. UWC8 Tree scale: 1 Bootstrap 100 80-99 50-80 AAI >= 55 % , AAI = 95% Pairwise Average Amino Acid Identity (AAI) (%) 34 40 46 52 58 64 70 76 82 88 95 0.00 0.00 0.00 0.00 0.00 0.00 0.00 97.55 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 99.95 0.00 85.84 86.31 0.00 85.42 86.02 99.98 100 100 0 10 20 30 40 50 60 70 80 90 100 Ca. B.v Ca. B.v F10 B7 A10 A8 B2 H10 G10 0.00 0.00 85.33 0.00 0.00 0.00 0.00 0.00 0.00 93.45 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 28.09 27.73 27.89 0.00 0.00 0.00 28.35 29.52 20.79 90.89 0 10 20 30 40 50 60 70 80 90 100 F10 B7 A10 A8 B2 H10 G10 Ca. B.v B7/F10 A8/A10 B2/G10/H10 117 (4.9) 184 (7.7) 156 (6.5) 299 (12.5) 302 (12.7) 1 (0.0) 633 (26.2) 6 (0.3) 19(0.8) 6 (0.3) 3 (0.1) 27 (1.1) 10 (0.4) 49 (2.1) 555 (23.2) Aligned Fractiom (AF) % Average nucleotide Identity (ANI) % Shared orthogroups Figure 3 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint Methionine degradation (M00035) Lysine biosynth., acetyl-DAP pathway, aspartate => lysine (M00525) Lysine biosynth., DAP aminotransferase pathway, aspartate => lysine (M00527) Lysine biosynth., DAP dehydrogenase pathway, aspartate => lysine (M00526) Lysine biosynth., succinyl-DAP pathway, aspartate => lysine (M00016) Citrate cycle (TCA cycle, Krebs cycle) (M00009) Citrate cycle, first carbon oxidation, oxaloacetate => 2-oxoglutarate (M00010) Citrate cycle, second carbon oxidation, 2-oxoglutarate => oxaloacetate (M00011) Gluconeogenesis, oxaloacetate => fructose-6P (M00003) Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate (M00001) Glycolysis, core module involving three-carbon compounds (M00002) Pentose phosphate pathway (Pentose phosphate cycle) (M00004) Pentose phosphate pathway, non-oxidative phase, fructose 6P => ribose 5P (M00007) PRPP biosynth., ribose 5P => PRPP (M00005) Pyruvate oxidation, pyruvate => acetyl-CoA (M00307) Glycogen biosynth., glucose-1P => glycogen/starch (M00854) Glyoxylate cycle (M00012) Methylcitrate cycle (M00982) UDP-N-acetyl-D-glucosamine biosynth., prokaryotes, glucose => UDP-GlcNAc (M00909) Undecaprenyl phosphate (UP) α -L-Ara4N biosynth., UDP-GlcA => UP α -L-Ara4N (M00761) Cytochrome bc1 complex respiratory unit (M00151) Cytochrome bd ubiquinol oxidase (M00153) Cytochrome c oxidase, prokaryotes (M00155) F-type ATPase, prokaryotes and chloroplasts (M00157) NADH (M00144) Succinate dehydrogenase, prokaryotes (M00149) CAM (Crassulacean acid metabolism), dark (M00168) CAM (Crassulacean acid metabolism), light (M00169) Reductive citrate cycle (Arnon-Buchanan cycle) (M00173) Reductive pentose phosphate cycle (Calvin cycle) (M00165) CMP-KDO biosynth. (M00063) KDO2-lipid A biosynth., Raetz pathway, LpxL-LpxM type (M00060) KDO2-lipid A biosynth., Raetz pathway, non-LpxL-LpxM type (M00866) Fatty acid biosynth., elongation (M00083) Fatty acid biosynth., initiation (M00082) Phosphatidylcholine (PC) biosynth., PE => PC (M00091) Phosphatidylethanolamine (PE) biosynth., PA => PS => PE (M00093) C1-unit interconversion, eukaryotes (M00141) Heme biosynth., animals and fungi, glycine => heme (M00868) Heme biosynth., bacteria, glutamyl-tRNA => coproporphyrin III => heme (M00926) Heme biosynth., plants and bacteria, glutamate => heme (M00121) Lipoic acid biosynth., eukaryotes, octanoyl-ACP => dihydrolipoyl-H (M00882) Lipoic acid biosynth., octanoyl-CoA => dihydrolipoyl-E2 (M00884) Lipoic acid biosynth., plants and bacteria, octanoyl-ACP => dihydrolipoyl-E2/H (M00881) Pimeloyl-ACP biosynth., BioC-BioH pathway, malonyl-ACP => pimeloyl-ACP (M00572) Siroheme biosynth., glutamyl-tRNA => siroheme (M00846) Adenine ribonucleotide biosynth., IMP => ADP,ATP (M00049) Adenine ribonucleotide degradation, AMP => Urate (M00958) Deoxyribonucleotide biosynth., ADP/GDP/CDP/UDP => dATP/dGTP/dCTP/dUTP (M00053) Guanine ribonucleotide biosynth., IMP => GDP,GTP (M00050) Pyrimidine deoxyribonucleotide biosynth., UDP => dTTP (M00938) Pyrimidine ribonucleotide biosynth., UMP => UDP/UTP,CDP/CTP (M00052) B2 H10 G10 A8 A10 B7 F10 Ca. B. v. 10 8 0 15 0 10 7 0 13 0 10 7 12 15 0 9 7 9 10 2 10 7 8 10 2 11 8 2 12 2 10 7 7 8 2 11 7 4 8 2 SEC-SRP Type III - Flagellar export apparatus Type IV - Conjugal transfer pilus assembly protein Type VI TAT 14 12 21 16 4 Amino acid metabolism Carbohydrate metabolism Energy metabolism Lipid metabolism Metabolism of cofactors and vitamins Nucleotide metabolism Secretion systems KEGG Module Completeness (%) 0 10 20 30 40 50 60 70 80 90 100 Figure 4 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint Leishmania Orhogroups in genus = 8,237 Orthogroups unique = 475 Trypanosoma Orthogroups in genus = 11,183 Orthogroups unique = 3,877 Bodo Orthogroups in genus = 13,816 Orthogroups unique = 7,321 Perkinsela sp. Bodo saltans Bodo B7 Bodo F10 Bodo A10 Bodo A8 Bodo H10 Bodo B2 Bodo G10 Trypanosoma vivax Trypanosoma brucei gambiense Trypanosoma brucei equiperdum Trypanosoma equiperdum Trypanosoma brucei brucei Trypanosoma grayi Trypanosoma melophagium Trypanosoma theileri Trypanosoma conorhini Trypanosoma rangeli Trypanosoma cruzi marinkellei Trypanosoma cruzi cruzi Trypanosoma cruzi Brazil Angomonas deanei Strigomonas culicis Phytomonas sp. EM1 Phytomonas sp. Hart1 Leptomonas pyrrhocoris Leptomonas seymouri Novymonas esmeraldas Porcisia hertigi Leishmania martiniquensis Leishmania sp. Namibia Leishmania sp. Ghana Leishmania orientalis Leishmania tarentolae Leishmania mexicana Leishmani major Leishmania donovani Leishmania infantum Leishmania naiffi Leishmania lindenbergi Leishmania braziliensis Leishmania utingensis Leishmania shawi Leishmania panamensis 0 5000 10,000 0 200 400 600 Support 0.8 - 1 0.5 - 0.79 < 0.5 Tree scale : 0.5 Pie charts show No. orthogroups in each genus Orthogroups shared with other genera Orthogroups unique to genus Orthogroups in genome Orthogroups unique to genome Figure 5 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint L. donovani L. infantum L. major anacixe m .L sisneilizarb .L sisne manap .L anahG .ps .L aibi maN .ps .L silatneiro .L L. martiniquensis P. hertigi T. brucei brucei esneib mag iecurb .T izurc .T Bodo B2 01G odoB Bodo H10 8A odoB Bodo A10 Bodo F10 Bodo B7 B. saltans KEGG module completeness (%) 0 10 20 30 40 50 60 70 80 90 100 Arginine biosynthesis, ornithine => arginine (M00844) Proline biosynthesis, glutamate => proline (M00015) Proline degradation, proline => glutamate (M00970) Proline metabolism (M00972) Tryptophan metabolism, tryptophan => kynurenine => 2-aminomuconate (M00038) Tyrosine degradation, tyrosine => homogentisate (M00044) Leucine degradation, leucine => acetoacetate + acetyl-CoA (M00036) Cysteine biosynthesis, serine => cysteine (M00021) Methionine degradation (M00035) Methionine salvage pathway (M00034) Histidine degradation, histidine => N-formiminoglutamate => glutamate (M00045) Lysine degradation, lysine => saccharopine => acetoacetyl-CoA (M00032) Glutathione biosynthesis, glutamate => glutathione (M00118) GABA biosynthesis, eukaryotes, putrescine => GABA (M00135) Polyamine biosynthesis, arginine => agmatine => putrescine => spermidine (M00133) Polyamine biosynthesis, arginine => ornithine => putrescine (M00134) Glycine cleavage system (M00621) Threonine biosynthesis, aspartate => homoserine => threonine (M00018) C10-C20 isoprenoid biosynthesis, non-plant eukaryotes (M00367) C5 isoprenoid biosynthesis, mevalonate pathway (M00095) Citrate cycle (TCA cycle, Krebs cycle) (M00009) Citrate cycle, first carbon oxidation, oxaloacetate => 2-oxoglutarate (M00010) Citrate cycle, second carbon oxidation, 2-oxoglutarate => oxaloacetate (M00011) Gluconeogenesis, oxaloacetate => fructose-6P (M00003) Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate (M00001) Glycolysis, core module involving three-carbon compounds (M00002) Pentose phosphate pathway (Pentose phosphate cycle) (M00004) Pentose phosphate pathway, non-oxidative phase, fructose 6P => ribose 5P (M00007) Pentose phosphate pathway, oxidative phase, glucose 6P => ribulose 5P (M00006) PRPP biosynthesis, ribose 5P => PRPP (M00005) Pyruvate oxidation, pyruvate => acetyl-CoA (M00307) Galactose degradation, Leloir pathway, galactose => alpha-D-glucose-1P (M00632) Glyoxylate cycle (M00012) Inositol phosphate metabolism, PI=> PIP2 => Ins(1,4,5)P3 => Ins(1,3,4,5)P4 (M00130) Malonate semialdehyde pathway, propanoyl-CoA => acetyl-CoA (M00013) Nucleotide sugar biosynthesis, glucose => UDP-glucose (M00549) Propanoyl-CoA metabolism, propanoyl-CoA => succinyl-CoA (M00741) UDP-N-acetyl-D-glucosamine biosynthesis, eukaryotes, glucose => UDP-GlcNAc (M00892) UDP-N-acetyl-D-glucosamine biosynthesis, prokaryotes, glucose => UDP-GlcNAc (M00909) Cytochrome bc1 complex respiratory unit (M00151) V-type ATPase, eukaryotes (M00160) C4-dicarboxylic acid cycle, NAD - malic enzyme type (M00171) C4-dicarboxylic acid cycle, NADP - malic enzyme type (M00172) C4-dicarboxylic acid cycle, phosphoenolpyruvate carboxykinase type (M00170) CAM (Crassulacean acid metabolism), dark (M00168) CAM (Crassulacean acid metabolism), light (M00169) Reductive pentose phosphate cycle (Calvin cycle) (M00165) N-glycan precursor biosynthesis (M00055) beta-Oxidation (M00087) beta-Oxidation, acyl-CoA synthesis (M00086) Fatty acid biosynthesis, elongation (M00083) Fatty acid elongation in endoplasmic reticulum (M00415) Acylglycerol degradation (M00098) Ceramide biosynthesis (M00094) Ketone body biosynthesis, acetyl-CoA => acetoacetate/3-hydroxybutyrate/acetone (M00088) Phosphatidylcholine (PC) biosynthesis, PE => PC (M00091) Phosphatidylethanolamine (PE) biosynthesis, ethanolamine => PE (M00092) Phosphatidylethanolamine (PE) biosynthesis, PA => PS => PE (M00093) Sphingosine biosynthesis (M00099) Triacylglycerol biosynthesis (M00089) Cholesterol biosynthesis, FPP => cholesterol (M00101) Ergocalciferol biosynthesis, FPP => ergosterol/ergocalciferol (M00102) C1-unit interconversion, eukaryotes (M00141) C1-unit interconversion, prokaryotes (M00140) Coenzyme A biosynthesis, pantothenate => CoA (M00120) Lipoic acid biosynthesis, plants and bacteria, octanoyl-ACP => dihydrolipoyl-E2/H (M00881) Molybdenum cofactor biosynthesis, GTP => molybdenum cofactor (M00880) NAD biosynthesis, aspartate => quinolinate => NAD (M00115) NAD biosynthesis, tryptophan => quinolinate => NAD (M00912) Riboflavin biosynthesis, plants and bacteria, GTP => riboflavin/FMN/FAD (M00125) Adenine ribonucleotide biosynthesis, IMP => ADP,ATP (M00049) Deoxyribonucleotide biosynthesis, ADP/GDP/CDP/UDP => dATP/dGTP/dCTP/dUTP (M00053) Guanine ribonucleotide biosynthesis, IMP => GDP,GTP (M00050) De novo pyrimidine biosynthesis, glutamine (+ PRPP) => UMP (M00051) Pyrimidine degradation, uracil => beta-alanine, thymine => 3-aminoisobutanoate (M00046) Pyrimidine deoxyribonucleotide biosynthesis, UDP => dTTP (M00938) Pyrimidine ribonucleotide biosynthesis, UMP => UDP/UTP,CDP/CTP (M00052) Amino acid metabolism Carbohydrate metabolism Energy metabolism Biosynthesis of terpenoids and polyketides Glycan metabolism Lipid metabolism Vitamine and cofactor metabolism Nucleotide metabolism Figure 6 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 3, 2025. ; https://doi.org/10.1101/2025.10.03.678719doi: bioRxiv preprint