Introduction

The kinetoplastids (Kinetoplastea; Euglenozoa) are a major group of heterotrophic protists that are of substantial medical and ecological importance, and that have bizarre basic cell biology: kinetoplast, a mitochondrial genome with massive DNA content, composed largely of thousands of small circular plasmids which encode guide RNAs required for editing of mitochondrial transcripts), paraxonemal rods, intronless nuclear genes transcribed in polycistronic clusters, a lack of transcriptional regulation of nuclear gene expression and glycosomes (modified peroxisomes where glycolysis occurs), among other features. (6; 49; 48; 52; 16; 32). The best studied kinetoplastids are the trypanosomatids, which include agents causing diseases in humans, live-stock and plants (51; 20). All other kinetoplastids are historically known as ‘bodonids’. Analyses of sequence data divide kinetoplastids into Prokinetoplastina and Metakinetoplastina (9; 46; 37; 11; 54). Metakinetoplastina includes Trypanosomatida plus three additional major taxa: Eubodonida, Parabodonida and Neobodonida (37; 26). 18S rDNA phylogenies divide neobodonids into subclades ‘1A’, ‘1B’ (Rhynchomonadidae), ‘1C’, ‘Nd6’, ‘1D’, and ‘1E’ (Allobodonidae) (55; 56; 18).

To date, phylogenies 18S rDNA sequences and 1-2 heat shock proteins have not resolved the deep branches within Metakinetoplastina. Unrooted analyses of Metakinetoplastina frequently place eubodonids adjacent to trypanosomatids, generally with poor support (47; 18). Rooted analyses occasionally place some/all neobodonids as the deepest branch within metakinetoplastids (47; 46; 18), indicating that eubodonids, then parabodonids, are successive sisters to trypanosomatids (Figure 1-A). Others place the root elsewhere: within trypanosomatids, or with trypanosomatids sister to all other metakinetoplastids (55; 56; 37) (Figure 1-C). While 18S rDNA datasets have by far the best taxon sampling, an extremely long basal branch separates metakinetoplastids from any outgroup (47; 18; 54). Predictably, the root position is highly unstable and any particular position poorly supported (18). This pattern is less extreme in protein datasets, however taxon sampling is much poorer and the majority of the neobodonid subgroups are missing (46; 13; 58).

Fig. 1:

Resolving the deep-level relationships within Metakinetoplastina is crucial to understand the evolution of the multiple bizarre traits found in kinetoplastids and for a rational high-level systematics of the group.

Results and Discussion

We sequenced and de novo assembled the transcriptomes of 11 neobodonid isolates. Novijibodo darinka and Allobodo yubaba, were described previously (40), whereas the remaining 9 isolates are undescribed Figure S1, Figure S2. We combined the predicted proteins of these transcriptomes with publicly available datasets representing Metakinetoplastina, Prokinetoplastina and discoban ougroups and performed de novo orthogroup inference. We selected single-copy orthogroups (OGs), aligned them individually and concatenated the alignments. The final phylogenomic matrix contained 144,376 sites from 444 de novo inferred OGs and 49 species, with an average OG occupancy of 90%, and an average of 86.9% of sites present per species Figure S9.

Maximum Likelihood (ML) and Bayesian analyses recovered identical phylogenies with maximum support at almost all nodes Figure 2; the exceptions were the two internal nodes within Prokinetoplastina, which were poorly resolved in the Bayesian analysis, while the node splitting 1C neobodonids and Rhynchomonadidae had poor support in the ML analysis Figure 2.

Fig. 2:

The monophyly of outgroup lineages (Diplonemida, Euglenida, Heterolobosea) and Trypanosomatida is consistent with previous studies (e.g. (41; 30; 58)), and supports the overall reliability of dataset construction and phylogenomics analysis workflows of our study. Our recovery of Prokinetoplastina and Metakine-toplastina as monophyletic sister taxa confirms previous analyses based on single to several dozens of genes (58; 9; 46; 11; 54).

Previous analyses specifically aiming to resolve deep kinetoplastid phylogeny have included up to two neobodonodis (13; 58). Our taxon sampling covers most of the 18S rDNA-based subgroups of neobodonids for the first time, and arguably represents the first real test of neobodonid monophyly using a phylogenomic approach. Our analyses are in fact consistent with previous studies (13; 58), but the much-improved taxon sampling of neobodonids, and larger number of OGs included (444 versus 64 (13) or 43 (58)) allowed a credible test of neobodonid monophyly, as opposed to inferring their placement under an implicit assumption that they represented a clade. Crucially, we find strong evidence that neobodonids are a paraphyletic assemblage, from which other metakinetoplastids descend as a monophyletic group.

The rooting of Metakinetoplastina splits them into two sister clades that both include ‘neobodonids’. One clade contains Trypanosomatida, Eubodonida, Parabodonida, Allobodonidae (’1E’) and neobodonid clade ’1D’; the second one contains the other three neobodonid clades that were included in the analysis: ‘1C’, ‘Nd6’, and ’1B’ (Rhynchomonadidae). The monophyly of each of ’1E’ (Allobodonidae), ’1D’, and ’1C’ is recovered with full support.

Within the first clade, Eubodonida then Parabodonida are recovered as the closest sister groups to Trypanosomatida, which is consistent with most relevant analyses of hsp90 sequences and with previous phylogenomic analyses (47; 46; 13; 58). The placement of Eubodonida (i.e. Bodo saltans) as the closest sister to trypanosomatids has been widely assumed for decades (4; 32), but actually was not clearly resolved by 18S rDNA analyses (47; 37; 18). Our analysis confirms this relationship, however, bolstering previous studies (47; 46; 13; 58).

The placement of Allobodonidae as sister to the Parabodonida, Eubodonida and Trypanosomatida clade is in harmony with the consistent bipartition between Allobodonidae and all other neobodonids in 18S rDNA trees (40; 18; 55; 56). The placement of 1D separate from Nd6, 1C and Rhynchomonadidae in our tree is neither supported nor contraindicated by previous analyses. No prior analyses of protein coding genes included 1D while analyses of 18S rDNA sequences show various sets of poorly supported relationships amongst 1A, Rhynchomonadidae, 1C, 1D and Nd6 (40; 18; 55; 56). The biology of 1D kinetoplastids remains poorly studied. The first isolate from 1D (referred to as ‘Cryptaulaxoides-like’) reportedly had a spiral groove running around the cell (55); Cryptaulaxoides (= Cryptaulax Skuja 1948) refers to elongate, rapidly swimming cells similar to Rhynchobodo spp. However, subsequent isolates were listed as ‘Neobodo sp.’ (56) and the 1D isolates studied here (BEAP0071-“10D”, BEAP0117-“13G”) superficially resemble Neobodo designis (Figure S1, Figure S2).

In the second major clade, the single representative of ’Nd6’ (KUC3-KIN1) falls sister to ’1B’ and ’1C’, although the node grouping ’1B’ and ’1C’ has poor support in the ML analysis. The poor support may be because both Nd6 and 1B are represented by a single species; sampling of more representatives is required. The grouping of Nd6, 1C and 1B is consistent with previous phylogenomic analyses placing the 1B and 1C clades together with strong support to the exclusion of other metakinetoplastids (58). Nd6 placement as a sister to clades 1B and 1C disagrees with some of the formerly produced 18S rDNA phylogenies that showed Nd6 branching with 1C (56) (yet the support in both cases is weak). In other 18S rDNA phylogenies, Nd6 is not recovered as the sister to 1C at all (18).

Our sampling of significant kinetoplastid clades is not exhaustive. The lack of 1A neobodonids is notable, since this clade includes Rhynchobodo. Phylogenetic analyses of 1 or 2 heat shock proteins have placed a Rhynchobodo isolate separately from other included neobodonids, often as the deepest branch within Metakinetoplastina, albeit statistical support was usually weak, and taxon sampling of both neobodonids and close outgroups was poor (46). If this position is correct, however, it implies that neobodonids are even more serially paraphyletic than our current results already show. A possible deep-branching position for Rhynchobodo is intriguing given that the recently documented free-living prokinetoplastids (54), are similar to Rhynchobodo: they are polykinetoplastidic cells with tubular extrusomes, large rostrum, a spiral groove or twist to the cell, and tend to swim through fluid rather than along surfaces (5), (39). It might be speculated that the common ancestor of all kinetoplastids was similar to free-living prokinetoplastids and Rhynchobodo. However, this might imply multiple independent origins of N. designis-like morphotypes: all representatives from 1C, Nd6 and 1D examined in this study are N. designis-like Figure S1, Figure S2. Because of this, including neobodonid clade 1A in future analyses would be especially valuable.

Phylogeny, sequences and cultures available from this study will be valuable to address questions on the ancestral states of metakinetoplastid kDNA organisation, glycolysis compartmentalisation, cytoskeleton organisation (as discussed by (47; 32; 2; 40; 59)), metabolic pathways, gene content and gene structure (addressed in (22; 8; 28), but with data from only one to three bodonids), and possible pre-adaptations of ancestral kinetoplastid entities to parasitism.

Materials and Methods

Sampling, cultivation and light microscopy

Isolates BEAP0154-“D44”, BEAP0117-“13G” and BEAP0071-“10D” were obtained from plankton samples of subsurface water near Blanes Bay, Mediterranean Sea, Spain (GPS coordinates 41.66667 N, 2.8 E) collected during the year 2022. Isolates BEAP0230-“RhMon” and BEAP0214-“G65133” were obtained from subsurface water of the Caribbean Sea, roughly 10 km south of the island of Grenada (GPS coordinates 11.893732 N, 61.699506 W), in February 2023. We enriched the plankton samples with the nutrient component of RS medium (dilutions from 1:5 to 1:500), and obtained stable mixed protist cultures from the enrichment samples, which were cultivated in RS medium (https://mediadive.dsmz.de/medium/P4, https://mediadive.dsmz.de/medium/P5) with nutrient (NM) to non-nutrient (NNM) component ratios between 1:50 and 1:500, at 17^◦C (for BEAP0154-“D44”, BEAP0117-“13G” and BEAP0071-“10D”) and at 23^◦C (BEAP0230-“RhMon” and BEAP0214-“G65133”). The isolates we used in these studies were monoclonal cultures obtained from initial stable mixed cultures by performing 3 to 5 sequential rounds of dilution to extinction. Isolates BEAP0117-“13G”, BEAP0154-“D44”, BEAP0230-“RhMon” and BEAP0214-“G65133” were all initially and routinely cultivated with a 12h:12h light:dark cycle, while isolate BEAP0071-“10D” was cultivated in darkness. We note that, if necessary, all of them are also capable of growth in darkness.

Isolate “OB23” was established from a seawater sample collected at a depth of 5 m at the entrance of Osaka Bay (GPS coordinates 34.324444 N, 135.120833 E), Japan, in July 2023. The sample was enriched in Daigo’s IMK medium (Fujifilm Wako Pure Chemical Corporation). Single cells were then isolated from the enrichment culture using a glass pipette, followed by several sequential rounds of serial dilution to extinction. The established culture was maintained at 20^◦C under dark conditions in Daigo’s IMK medium.

Freshwater isolates “BAB10-KIN1”, “INU10-KIN1”, and “KUC3-KIN1” were isolated from water samples collected in October 2023 in the Masurian Lakeland, Poland. “BAB10-KIN1” originated from a sample of the hypolimnion of Lake Babiety (53.719288 N, 21.129903 E), while “INU10-KIN1” and “KUC3-KIN1” were isolated from samples of the epilimnion of Lake Inulec (53.805704 N, 21.491791 E) and Lake Kuc (53.815844 N, 21.398909 E), respectively. To enrich for kinetoplastids, 10 mL of each sample were incubated with an autoclaved rice grain for 7 days. Individual kinetoplastid cells were manually picked, together with associated bacteria, using a glass microcapillary and transferred to separate wells of 12-well microplates containing NM component diluted 1:100 in Bold’s Basal Medium with threefold nitrogen and vitamins (3N-BBM (Bischoff, 1963); https://mediadive.dsmz.de/medium/C15). The isolates were routinely maintained in 3N-BBM supplemented with 1:100 dilution of NM at +18^◦C in the dark.

The isolation of Novijibodo darinka-“Gemkin” and Allobodo yubaba BEAP0186 “16Ckin” was described previously, in (40).

All isolates in this study were grown as xenic cultures (i.e., with several uncharacterized bacterial species). Light microscopy data from selected isolates used in this study (BEAP0230-“RhMon”, BEAP0214-“G65133”, BEAP0154-“D44”, BEAP0117-“13G”, BEAP0071-“10D”, “BAB10-KIN1”, “INU10-KIN1” and “KUC3-KIN1”) was obtained with a Zeiss Axio Observer.Z1 inverted microscope with a Differential Interference Contrast (DIC) Plan-Apochromat 63x Oil immersion objective (NA=1.4). The imaging setup was equipped with a HDCamC13440-20CU Hamamatsu camera coupled with ZEN image acquisition software. Imaging data was further processed in Fiji (44). Transmission Electron Microscopy was performed at Electron Cryomicroscopy Unit at the Scientific and Technical Centers of the University of Barcelona. Cells were collected by centrifugation (2000 x g for 5 min), resuspended in 20% BSA in artificial sea water and cryo-immobilized using a Leica HPM100 high-pressure freezer (Leica Microsystems, Vienna, Austria). Samples were freeze-substituted with 2% osmium tetroxide and 0.1% uranyl acetate dissolved in the pure acetone at -90^◦C for 72 hours in an EM AFS2 (Leica Microsystems, Vienna, Austria), followed by gradual warming up to RT and infiltration of Epon-812 resin and polymerization in Epon-812 at 60^◦C for 48 hours. Ultrathin Sections of 60 nm were obtained with a UC6 ultramicrotome (Leica Microsystems, Vienna, Austria) and placed on Formvar-coated copper grids. Sample Sections were stained with 2% uranyl acetate for 30 min, lead citrate for 5 min and examined in a TEM Jeol JEM 1010 (Gatan, Japan) equipped with a tungsten cathode. Images were acquired at 80 kV with a 1k x 1k CCD Megaview III camera. Complete imaging data, including movies, are available via Figshare at the link https://doi.org/10.6084/m9.figshare.30214630.

RNA isolation and transcriptome sequencing

The Novijibodo darinka isolate “Gemkin” was grown to high abundance in a 400 ml parafilm-sealed vented culture flask, initially containing 20 ml of cell culture, 100 ml of CR10 media ((17), with CR number indicating salinity in ppt), and 12 autoclave-sterilised wheat grains. The culture was grown for 12 days in a 21^◦C dark incubator before being pelleted via centrifugation (1500 x g for 12 min, then 6500 x g for 3 min). The pelleted cells were lysed using 4 ml TRIzol^TM (ThermoFisher), and stored at -20^◦C for 2 weeks. Total RNA was isolated using the standard TRIzol^TM protocol, except that the RNA was precipitated at -20^◦C overnight. Total RNA was washed twice using the RNeasy PowerClean Pro cleanup kit (QIAGEN).

The RNA was stored at -80^◦C before being sent to Génome Québec for poly-A selection, cDNA library preparation, and paired-end sequencing (100 bp) using Illumina NovaSeq.

Isolate OB23 was cultivated as described above in eight 50-mL flasks (Violamo) containing Daigo’s IMK medium. After two weeks of cultivation, the cells were pelleted by centrifugation at 2,000 x g for 10 min. RNA was extracted from the pelleted cells using TRIzol reagent and Phasemaker Tubes (Thermo Fisher Scientific) according to the manufacturer’s instructions. Total RNA was then used for library construction using the MGIEasy Fast RNA Library Prep Set (MGI Tech). The resulting libraries were circularized using the MGIEasy Dual Barcode Circularization Kit (MGI Tech) and converted into DNA nanoballs (DNBs) using the DNBSEQ-G400RS High-throughput Sequencing Kit and High-throughput Pair-End Sequencing Primer Kit (App-D) (MGI Tech). Finally, 150-bp paired-end sequencing was performed on a DNBSEQ-G400 platform (MGI Tech).

Total RNA from Allobodo yubaba (isolate BEAP0186-“16Ckin”), and from the remaining isolates used in this study (BEAP0230-“RhMon”, BEAP0214-“G65133”, BEAP0154-“D44”, BEAP0117-“13G”, BEAP0071-“10D” as well as “BAB10-KIN1”, “INU10-KIN1” and “KUC3-KIN1”) was extracted following the procedure from (42). All of these isolates used in transcriptomics studies were strictly monoclonal (i.e., they were initiated from a single eukaryotic cell and any co-isolated environmental bacteria), except for “INU10-KIN1”, for which the culture originated from a few manually picked cells . For RNA isolation, cultures were grown in RS medium (https://mediadive.dsmz.de/medium/P4, https://mediadive.dsmz.de/medium/P5), of nutrient (NM) to non-nutrient (NNM) component ratios between 1:20 and 1:500 for marine isolates (BEAP0230-“RhMon”, BEAP0214-“G65133”, BEAP0154-“D44”, BEAP0117-“13G”, BEAP0071-“10D”) or 1:100 dilutions of NM component prepared in distilled water in 3N-BBM (Bold’s Basal Medium) (Bischoff, 1963) for freshwater isolates (“BAB10-KIN1”, “INU10-KIN1” and “KUC3-KIN1”), in 75 cm² rectangular canted-neck cell culture flasks with vented caps (353136, Corning Life Sciences) and pelleted by centrifugation (13000 x g for 20 min, at +4^◦C). RNA isolation was performed with an RNAqueous kit (AM1912, ThermoFisher Scientific), followed by genomic DNA digestion using a TURBO DNA-free kit (AM1907, ThermoFisher Scientific). The kits were used according to the manufacturer’s protocols with the modifications described in (42). Double the normal rounds of polyA selection were applied to the total RNA before cDNA synthesis (to reduce the amount of prokaryotic RNA in the sample, which is common in non-axenic cultures). The integrity of total RNA was checked using Bioanalyzer 2100 RNA Pico chips (Agilent Technologies) at the CRG Genomics Core Facility in Barcelona, and total RNA was stored at -80^◦C before being sent to the same facility for poly-A selection, library preparation, and sequencing. Stranded libraries were prepared using the NEBNext Ultra II Directional RNA Library Prep Kit. Finally, the RNA was sent for paired-end sequencing (150 bp) on an Illumina NextSeq 2000.

Transcriptome assembly, decontamination, protein prediction and redundancy reduction

The procedure described below and in Figure S3 was performed for all the transcriptomes obtained in this study (isolates “Gemkin”, “OB23”, BEAP0186-“16Ckin”, BEAP0230-“RhMon”, BEAP0214-“G65133”, BEAP0154-“D44”, BEAP0117-“13G”, BEAP0071-“10D”, “BAB10-KIN1”, “INU10-KIN1” and “KUC3-KIN1”).

Read trimming and transcriptome assembly

Raw reads were trimmed using fastp(v0.20.1) (12); depending on the sequencing technology used, settings varied (see exact command list at https://github.com/beaplab/kineto_phylo/blob/main/MASTER_commands_V2.sh). Assembly was performed using RNAspades v3.13.1 (7), in stranded mode (with the exception of “OB23”, which was assembled in non-stranded mode), with the other settings left as their defaults.

Identification of putative contaminants and creating the database of contaminant genomes

Potential contaminants were primarily identified by searching for 18S ribosomal RNA gene sequences of the non-target eukaryote species and for prokaryotic 16S ribosomal sequences in the assemblies.

18S and 16S of contaminants were retrieved by blasting 18S and 16S sequences of Mus musculus (NR 003278) and Escherichia coli (NR 114042) against the assembled contigs. Matches were manually filtered out, and target contigs producing reasonable match length, E-value and % identity were selected (see table on FigShare at https://doi.org/10.6084/m9.figshare.30214630 for the exact list of contigs and contaminants). These contigs that putatively contained 18S or 16S of contaminants were used as queries to perform blastn against the NCBI nt database (a custom Perl script was used as a wrapper for performing blastn of multiple queries vs. multiple targets, available on Github https://github.com/beaplab/kineto_phylo/blob/main/Suppl_v1_scripts/web_blast.pl). A single top hit was selected for each contig. NCBI taxonomy IDs of the targets from the top hits were retrieved.

In addition to the taxonomy IDs matching 16S- and 18S-containing contigs identified by the BLAST approach, we also retrieved taxonomy IDs of 18S sequences assembled by phyloFlash v3.3b3 (19) from the raw transcriptomic reads (“Full-length SSU rRNA sequences assembled by SPAdes, matched to SILVA database with Vsearch” list of PhyloFlash output).

Genomes corresponding to all the taxonomy IDs of potential contaminants were downloaded using the Entrez v14.6 and NCBI datasets v16.22.1 tools. If there were multiple genomes available for a single taxonomy ID, up to 5 genomes were selected for download, and RefSeq genomes were given preference.

Removal of contamination

For decontamination, assembled contigs were used as blastn queries against the database of putative contaminant genomes obtained in the previous step (a custom Perl script was used as a wrapper for performing blastn of multiple queries vs. multiple targets; it is available on Github https://github.com/beaplab/kineto_phylo/blob/main/Suppl_v1_scripts/run_BLAST.pl).

The distribution of % identity, E-value and match length was visualised for blast hits of 60,000 randomly selected contigs from all assemblies vs. putative contaminant genomes Figure S5 using a custom R script available on Github https://github.com/beaplab/kineto_phylo/blob/main/Suppl_v1_scripts/blast_res_viz.Rmd.

The minimum threshold for blastn hit match length was set to 150 for all assembly-contaminant genome pairs, based on visual inspection of the aforementioned distribution. Minimum thresholds of % identity were selected individually for each pair of assembly and contaminant genome, based on an automated algorithm designed to select a value separating identical or nearly identical matches resulting from contamination versus true matches resulting from sequence homology (see the same custom R script for details).

Contigs that matched contaminant genomes with % identity higher than the automatically identified cutoff and with the length of match exceeding 150 bp, were removed from the corresponding assemblies (using seqkit v0.10.0 (45)).

Removal of Cross-contamination

To remove potential cross-contamination among transcriptomes which were sequenced together on the same flow cell (BEAP0230-“RhMon”+BEAP0214-“G65133”, BEAP0186-“16Ckin”+BEAP0154-“D44”+BEAP0117-“13G”+BEAP0071-“10D” and “BAB10-KIN1”+“INU10-KIN1”+“KUC3-KIN1” were sequenced together), the WinstonCleaner https://github.com/kolecko007/WinstonCleaner tool was run with default parameter values. WinstonCleaner removes cross-contamination by mapping trimmed reads from the species sequenced together against the assembled contigs for each species in order to estimate coverage, and when a pair of contigs from different assemblies is found to have a high sequence similarity, the contig with lower read coverage is considered a contaminant and is removed from the corresponding assembly. The cutoffs of minimum % identity for the overlapping contigs to be considered cross-contaminants were automatically identified for each assembly pair and additionally verified by manually inspecting the distribution of % identity of all contigs in assembly pairs.

Protein prediction and assembly completeness

Protein prediction was done with the TransDecoder version v5.0.1 tool run in stranded mode (with the exception of “Gemkin” and “OB23”, for which it was run in non-stranded mode) https://github.com/TransDecoder/TransDecoder to produce predicted CDSs and proteins for each decontaminated assembly.

Obtaining and processing data from other species

Data for other species used in the phylogenomic inference was retrieved from various sources (see table on FigShare at https://doi.org/10.6084/m9.figshare.30214630). Raw reads for Cruzella marina (ERR13669942) and Azumiobodo hoyamushi (SRR10586159) were downloaded from NCBI and assembled with RNAspades, in stranded mode, with other settings default. Assembled contigs of Neobodo designis and Trypanoplasma/Cryptobia borreli were obtained from EukProt (43). Assemblies of these four kinetoplastid species (Neobodo designis, Trypanoplasma/Cryptobia borreli, Cruzella marina, Azumiobodo hoyamushi) were cleaned using the procedure to remove contamination described above. Assembled contigs of the outgroup species Pharyngomonas kirbyi were also obtained from EukProt (43).

For the data from Tsukubamonas globosa and Ploeotia vitrea, the same procedures were applied as for Pharyngomonas kirbyi. However, these two species were excluded from the final set of taxa used in phylogenomic inference, because of low assembly completeness.

Proteins were predicted as described above for Cruzella marina, Azumiobodo hoyamushi, Neobodo designis, Trypanoplasma/Cryptobia borreli and Pharyngomonas kirbyi (further details are provided in Figure).

Proteomes/CDS from the rest of the species were downloaded from the sources listed in table on FigShare at https://doi.org/10.6084/m9.figshare.30214630, CDSs were translated to amino acid sequences, if necessary, and used in the subsequent procedures in an unaltered state (without additional decontamination or quality control) (see Figure).

Redundancy reduction

A redundancy-reduced proteome set was generated from the initial set of proteomes using CD-HIT v4.8.1 (31; 15). The redundancy-reduced set (“NR90”) retained proteins with < 90% amino acid global sequence identity for each proteome (the default identity metric is calculated by CD-HIT as the number of aligned identical amino acids divided by the full length of the shorter sequence (31; 15)).

Orthology inference and verification

The final dataset used for orthology inference contained 11 new predicted proteomes produced in this study, and 38 more predicted proteomes: 28 from kinetoplastid species and 10 from other Discoba lineages (as outgroups), resulting in 49 proteomes in total.

Orthogroups (OGs) were inferred using OrthoFinder v2.5.5 (14), using the “diamond ultra sens” option for the sequence search algorithm and the “msa” option for gene tree inference, with other options left as their defaults.

Initially, 72407 OGs were identified in the set of 49 proteomes. The subset of OGs used for single-OG tree reconstruction and further manual curation was selected based on the species represented and on the number of proteins per species in each OG. Within each OG, we estimated the mean and the standard deviation of the number of proteins present per species, and we counted the number of species with at least one protein. Based on the distributions of these three parameters across all 72407 OGs (see Figure S8), we identified a set of thresholds which would theoretically allow us to produce a subset OGs containing the highest number of single-copy orthologs (i.e., one copy per species) in the largest number of taxa. The summary statistics on the OGs were calculated using a custom R script available on Github https://github.com/beaplab/kineto_phylo/blob/main/Suppl_v1_scripts/Orthofinder_stats_UNIVERSAL_V2.Rmd.

Our thresholds resulted in the selection of OGs meeting all three of the following criteria: (i) mean number of proteins per single species between 0.9 and 1, inclusive, (ii) standard deviation of the number of proteins per single species less than or equal to 0.5, and (iii) number of species with at least one protein present greater than 35 (out of 49 species in total).

The subset obtained by applying these criteria contained 461 OGs. In order to produce a set of single-copy orthologous sequences, first, proteins within each OG were aligned with MAFFT v7.490 (24), using the L-INS-i (“–localpair”) strategy, with the iteration limit set to 16. Resulting alignments were trimmed using trimAl v1.4.rev22 (10) with the “–gappyout” option, and maximum likelihood single-OG trees were inferred using IQ-TREE v1.6.12 (38) with node supports estimated by 1000 ultrafast bootstrap replicates, and other options as their defaults.

Each single-OG tree was manually examined to mark inparalogs, contaminant sequences, and other artefacts for removal from the corresponding OG. As a result of single-gene tree checks, 17 OGs were completely removed from further analysis, 377 OGs were edited by removing some of the proteins, and 67 OGs were left unedited.

Protein sequences from the set of 461 OGs, alignments, single-OG trees, as well as the lists of sequences and OGs that were removed by manual inspection, can be found at table S3 and “orthogroup” directory on FigShare at https://doi.org/10.6084/m9.figshare.30214630.

Ortholog alignment and trimming

After paralog and artefact removal, proteins of each OG were re-aligned with MAFFT v7.490 (24), using the L-INS-i (“–localpair”) strategy, with the iteration limit set to 1000. Resulting alignments were trimmed using trimAl v1.4.rev22 (10) with the “–gappyout” option. Trimmed alignments of all OGs were concatenated using PhyKit v2.0.1 (50).

Phylogenomic tree inference

Bayesian tree inference was performed with PhyloBayes-mpi v1.9 (29). Five independent chains were run with the CAT+GTR+G4 model, sampling once every generation. The run was stopped at the 16000th generation, at which point the five chains did not converge, however the only differences between them were in the inter-relationships among the three prokinetoplastids. A consensus tree topology was obtained with a burn-in of 4000 generations (25% of the total number of generations per chain). The mean discrepancy across all bipartitions of all three chains was 0.00631579.

A maximum likelihood analysis was run using IQ-TREE v2.3.6 (36). The LG4M+R10 model was determined by the IQ-TREE ModelFinder (23) as the best-fitting, according to BIC, as well as AIC and AICc. The robustness of the most likely ML tree with the LG4M+R10 model was inferred using ultrafast bootstrapping (1000 replicates) (21) along with standard bootstrapping (200 replicates).

Trees were visualised in tvBOT (57) and edited for presentation with Inkscape.

Data availability

Raw transcriptome sequencing reads for BEAP0186-“16Ckin”, BEAP0230-“RhMon”, BEAP0214-“G65133”, BEAP0154-“D44”, BEAP0117-“13G”, BEAP0071-“10D”, “BAB10-KIN1”, “INU10-KIN1” and “KUC3-KIN1” and “Gemkin” are available via GenBank under the BioProject PRJNA1332963. Raw transcriptome sequencing reads for “OB23” are available via GenBank under the BioProject PRJDB40754. In addition, 18S ribosomal RNA gene sequences obtained from transcriptomic assemblies of BEAP0230-“RhMon”, BEAP0214-“G65133”, BEAP0154-“D44”, BEAP0117-“13G”, and BEAP0071-“10D”, “BAB10-KIN1”, “INU10-KIN1” and “KUC3-KIN1” were deposited in GenBank under IDs PZ229093-PZ229100. The 18S ribosomal RNA gene sequence from “OB23” was deposited in GenBank under IDs LC924540.

A BEAP0186-“16Ckin” cell culture (with mixed bacteria) is available in the Roscoff Culture Collection (RCC) with accession number RCC11336. The following cell cultures (with mixed bacteria) are available in the Culture Collection of Algae and Protozoa (CCAP): BEAP0230-“RhMon” (accession number CCAP 1982/5), BEAP0214-“G65133” (CCAP 1982/4), BEAP0154-“D44” (CCAP 1982/3), BEAP0117-“13G” (CCAP 1982/2), and BEAP0071-“10D” (CCAP 1982/1). All commands that were run during data analysis are stored in the file available on GitHub https://github.com/beaplab/kineto_phylo/blob/main/MASTER_commands_V2.sh. Assembled transcriptome contigs and predicted proteins for all newly sequenced isolates in this study are available on FigShare at https://doi.org/10.6084/m9.figshare.30214630.

All files generated in the process that cannot be directly obtained by running the commands from the aforementioned files (for example, the manually curated sets of orthologous genes) are available either through the GitHub repository https://github.com/beaplab/kineto_phylo or Figshare at https://doi.org/10.6084/m9.figshare.30214630.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 949745), from the Departament de Recerca i Universitats de la Generalitat de Catalunya (exp. 2021 SGR 00751), from grant PID2023-152955NA-I00 funded by MICIU/AEI/10.13039/501100011033 and by ERDF/EU (to D.J.R.), from the Natural Sciences and Engineering Research Council of Canada (NSERC) discovery grant 298366-2019 (to A.G.B.S.), from Grants-in-Aid for Scientific Research (No. 21H05057, No. 24K21929, and No. 25K02334) from the Japan Society for the Promotion of Science (JSPS), from the project No. 2022/47/P/NZ8/02074 (to DT) cofunded by the National Science Centre and the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 945339. RK is also supported by Institute for Fermentation, Osaka (2024-G-1-011); some of the environmental samples were collected thanks to funding from the National Science Centre (Poland), grant number 2020/37/B/NZ8/01456 (to Anna Karnkowska).

Computational analyses were enabled by the supercomputing infrastructure from the Galician Supercomputing Center (CESGA; funded by the Spanish Ministry of Science and Innovation, the Galician Government and the European Regional Development Fund; ERDF).

We are deeply grateful to Clare Morrall for arranging access to sampling and facilities in Grenada. We thank Clara Cardelús, operating the Blanes Bay Microbial Observatory (BBMO), and Josep M. Gasol and Ramon Massana for providing access to monthly water samples from Blanes Bay. We thank Margarita Skamnelou and À lex Gàlvez i Morante from the Biology and Ecology of Abundant Protists lab for isolating, identifying and establishing initial cultures of isolates of BEAP0154-“D44”, BEAP0230-“RhMon” and BEAP0214-“G65133”. We thank Dr. Keigo Yamamoto (Research Institute of Environment, Agriculture and Fisheries, Osaka Prefecture, Osaka, Japan) for arranging access to sampling and facilities in Osaka Bay. We thank Anna Karnkowska for providing environmental samples from which isolates “BAB10-KIN1”, “INU10-KIN1”, and “KUC3-KIN1” were isolated. We also thank Elizabeth Weston (Dalhousie University) for training JP in RNA extraction and maintaining cultures.

Additional files

Supplemental information.

Additional information

Funding

EC | Horizon 2020 Framework Programme (H2020)

https://doi.org/10.3030/949745

• Daryna Zavadska

• Daniel J Richter

Government of Catalonia | Departament de Recerca i Universitats, Generalitat de Catalunya (DURSI) (2021 SGR 00751)

• Daniel J Richter

• Daryna Zavadska

MICIU/AEI/10.13039/501100011033 (PID2023-152955NA-I00)

• Daniel J Richter

• Daryna Zavadska

ERDF/EU (PID2023-152955NA-I00)

• Daniel J Richter

• Daryna Zavadska

Natural Sciences and Engineering Research Council of Canada (NSERC) (298366-2019)

• Alastair Simpson

• Julia A Packer

Japan Society for the Promotion of Science (21H05057)

• Ryoma Kamikawa

MEXT | Japan Society for the Promotion of Science (JSPS) (24K21929)

• Ryoma Kamikawa

MEXT | Japan Society for the Promotion of Science (JSPS) (25K02334)

• Ryoma Kamikawa

National Science Centre Poland (2022/47/P/NZ8/02074)

• Daria Tashyreva