{'fraction_non_alphanumeric': 0.03896960285890536, 'fraction_numerical': 0.021548148305715685, 'mean_word_length': 5.176037834997373, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 12, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gene fragmentation and RNA editing without borders: eccentric mitochondrial genomes of diplonemids

Abstract Diplonemids are highly abundant heterotrophic marine protists. Previous studies showed that their strikingly bloated mitochondrial genome is unique because of systematic gene fragmentation and manifold RNA editing. Here we report a comparative study of mitochondrial genome architecture, gene structure and RNA editing of six recently isolated, phylogenetically diverse diplonemid species. Mitochondrial gene fragmentation and modes of RNA editing, which include cytidine-to-uridine (C-to-U) and adenosine-to-inosine (A-to-I) substitutions and 3′ uridine additions (U-appendage), are conserved across diplonemids. Yet as we show here, all these features have been pushed to their extremes in the Hemistasiidae lineage. For example, Namystynia karyoxenos has its genes fragmented into more than twice as many modules than other diplonemids, with modules as short as four nucleotides. Furthermore, we detected in this group multiple A-appendage and guanosine-to-adenosine (G-to-A) substitution editing events not observed before in diplonemids and found very rarely elsewhere. With >1,000 sites, C-to-U and A-to-I editing in Namystynia is nearly 10 times more frequent than in other diplonemids. The editing density of 12% in coding regions makes Namystynia’s the most extensively edited transcriptome described so far. Diplonemid mitochondrial genome architecture, gene structure and post-transcriptional processes display such high complexity that they challenge all other currently known systems.


INTRODUCTION
Diplonemids are heterotrophic marine flagellates belonging to the phylum Euglenozoa, which also includes the well-studied parasitic kinetoplastids and free-living euglenids (1,2) ( Figure 1A). Diplonemids have been largely overlooked due to technical limitations, because the SSU rRNA V4 region, typically amplified in the metabarcoding approach, has expanded beyond typical lengths in diplonemids. In a recent survey, which targeted the more conserved V9 region, they have been detected in virtually every sample of seawater (3) and are currently ranked among the most diverse and abundant eukaryotic groups in the world's oceans (4)(5)(6)(7). The ecological role of diplonemids in marine environments has only recently begun being appreciated (6,8,9).
The most conspicuous features of diplonemids are their unique and complex mitochondrial genome architecture  (25). Primary transcripts are end-processed, removing 5 and 3 non-coding regions from the primary transcripts. Certain module transcripts undergo substitution RNA editing and/or appendage RNA editing (nucleotide additions at the module's 3 end). The module transcript that will constitute the transcript's 3 end is poly-adenylated (mRNAs and mtLSU rRNA) or poly-uridylated (mtSSU rRNA) (16,19). Finally, modules are joined together (trans-spliced) yielding mature RNA (mRNA or rRNA). Note that all post-transcriptional processes (gray background) occur in parallel in the diplonemid mitochondrion (19); thus, the arrows do not imply strict sequentiality. and gene expression ( Figure 1B), studied in depth in the type species Diplonema papillatum (reviewed in (15)(16)(17)). For example, the amount of mitochondrial DNA (mtDNA), estimated at 250 Mbp, is so far the highest recorded for an organelle (18). Further, its mtDNA is composed of >80 covalently closed non-catenated 6 and 7 kbplong circular chromosomes. Except for a short unique region called the 'cassette', chromosomes consist mostly of repetitive sequence termed 'constant region', which is essentially identical across chromosomes of a given size class ( Figure 1B). It is the unique cassette that typically encloses a single gene fragment, also called module, the size of which ranges from 40 to 540 bp. Chromosomes including cassettes are transcribed separately, then the non-coding portions are removed leaving behind module-only transcripts that are subsequently joined to their cognate neighboring modules derived from other circles (16,19). In this way, mature transcripts (mRNAs and rRNAs) are assembled via massive trans-splicing, the mechanism of which remains unknown.
Trans-splicing in diplonemid mitochondria differs from that observed in organelles and in the nucleus of nematodes and other eukaryotes including diplonemids (20)(21)(22), because the former process is apparently catalyzed by neither spliceosomes nor Group I or Group II splicing machineries (16,19,(23)(24)(25)(26). With the exception of the mitochondrial small subunit ribosomal RNA (mtSSU rRNA), all D. papillatum genes undergo this assembly process, making the extent of trans-splicing unprecedented.
Moreover, in addition to gene fragmentation compensated by trans-splicing, D. papillatum mitochondrial transcripts are subject to extensive RNA editing of two fundamentally different types: post-transcriptional uridine addition (U-appendage) at 3 ends of certain modules (unique to diplonemids), and deaminations of adenosines to inosine (A-to-I) and cytidines to uridines (C-to-U) at numerous positions within coding regions (16,27).
Studies of three other diplonemid species (D. ambulator, Flectonema neradi and Rhynchopus euleeides) showed little deviation from the features observed in the type species (27,28), except that across these taxa, the size of mitochondrial chromosomes ranges from 2 to 12 kbp (27,28). However, it was reported recently that gene fragmentation, as well as U-additions and A and C substitutions in transcripts are much more frequent in Hemistasia phaeocysticola, the single hemistasiid species examined until now (27,29). While D. papillatum has at most 11 modules per gene (16), the genes of H. phaeocysticola are fragmented twice as much (29).
Does the single examined hemistasiid species represent an exceptional case, or can fragmentation and RNA editing of mitochondrial transcripts reach even higher levels of complexity? To address this question, we examined species that have recently become available in culture, with their morphology, ultrastructure and life cycles described (11)(12)(13). Here, we present the most extensive comparative study of diplonemid mtDNA performed thus far. We show that the degree of mitochondrial RNA editing and gene fragmentation can reach unprecedented complexity, highlighting several questions about the role and evolution of these remarkable features.

Strains, culture conditions and nucleic acids extraction
The six diplonemid species used in this study (Supplementary Table S1) were recently isolated from marine water collected in aquaria, lagoons and sandy beaches of Japan (11)(12)(13). The species were axenically cultivated in a medium containing 3.6% sea salts (Sigma-Aldrich, S9883), supplemented with 1% (v/v) heat-inactivated horse serum (Sigma-Aldrich, H0146) and 0.025 g/l LB broth powder (Sigma, L3522). The medium was filter-sterilized using a 0.22-m filter.
Total DNA from exponentially growing cultures was isolated using MasterPure Complete DNA and RNA Purification Kit (Lucigen, MC85200) specially designed for the isolation of DNA from marine organisms. RNA was extracted from whole cells using TriReagent (MRC, TR118) to prevent the loss of small RNAs corresponding to processing and trans-splicing intermediates. Residual DNA was removed by DNase treatment followed by extraction with a homemade Trizol substitute (30).

Reverse transcription, RT-PCR, 5 and 3 RACE, and poly-A tail site mapping
Reverse transcription was performed with First Strand cDNA Synthesis Kit for subsequent PCR (Roche) or with SuperScript IV Reverse Transcriptase (Thermo). Complementary DNA was amplified with Q5 High-Fidelity DNA Polymerase (New England Biolabs). PCR products were purified using the QIAquick Gel Extraction kit (Qiagen), Wizard SV Gel and PCR Clean-Up system (Promega) or Monarch DNA Gel Extraction Kit (New England BioLabs). Mapping of 3 polyadenylation sites was performed using an oligo-dT primer and a gene-specific primer.
To determine the 5 and 3 ends of modules (5 and 3 RACE), the RNA adapter-oligonucleotide dp124 and the 5 RACE Adapter (from FirstChoice RLM-RACE Kit, Invitrogen) was ligated to the RNA using T4 RNA ligase I (New England Biolabs) and EcRtcB RNA ligase (New England Biolabs), respectively. Detailed protocols are available at https://www.protocols.io/researchers/matus-valach. RT-PCR was performed using specific primers, and amplicons were sequenced at the IRIC Genomics Core Facility (Montreal, Canada) or at Eurofins Genomics (Ebersberg, Germany). Primer and adaptor sequences are listed in Supplementary Table S2.

Genome and transcriptome sequencing and assembly
Both library preparation and sequencing of genomes and transcriptomes were outsourced to the Genome Quebec Innovation Centre (Montreal, Canada). Single DNA-Seq and RNA-Seq library per species were produced due to limited material availability. Illumina genomic paired-end libraries were constructed from total DNA and sequenced in a single Illumina MiSeq lane. DNA reads were assembled using SPAdes v3.11.1 (31) and alternatively, the Tadpole assembler (part of the BBTools suite; https://jgi.doe.gov/dataand-tools/bbtools/).
To avoid the huge variety of trans-splicing and RNA editing intermediates present in diplonemid mitochondria (16,19), which complicate analysis and interpretation, we opted for the enrichment of mature mitochondrial transcripts, i.e. the polyadenylated (poly-A) RNA fraction, which was isolated from total RNA to construct strandspecific RNA-Seq libraries with an average insert size of ∼200 nt. Libraries were sequenced on an Illumina HiSeq platform. For de novo assembly, Trinity v2.2.0 software was used with default parameters (32). Read counts and lengths for both DNA and RNA sequencing are listed in Supplementary Table S1. The raw sequencing data are available at NCBI (https://www.ncbi.nlm.nih.gov/) as BioProject PR-JNA525750.
Two of the examined species (D. japonicum and N. karyoxenos) contain endosymbiotic bacteria (12,13). For the purpose of this study focusing on mitochondrial sequences, it was not necessary to estimate relative abundance of bacterial sequences in the datasets. However, we did perform RNA-Seq read mapping to DNA contigs and found only negligible differences between mapping rates to sequences of endosymbiont-bearing and -lacking species, which suggested that possessing an endosymbiont did not introduce any significant bias to our strategy.

Identification of transcripts and annotation of genomic modules
Candidate contigs originating from the mitochondrial transcriptome were identified by BLASTx searches (33) using protein sequences of previously identified mature mitochondrial mRNAs from D. papillatum, D. ambulator, H. phaeocysticola, F. neradi and R. euleeides as queries, and the transcriptome assemblies of the new diplonemids as queried databases. Genomic modules for each species were annotated based on BLASTn searches using predicted transcript sequences as queries. To validate module assignments, modules were aligned with the respective transcript using the built-in aligner of the Geneious 10.1.3 software (34) and visually inspected. To infer protein-coding ORFs, the nucleotide sequences were conceptually translated using NCBI's genetic code Table 4 (TGA = Trp). Identified modules are cataloged in Supplementary Table S3.

Completion of mitochondrial transcripts
RNA-Seq reads were mapped to reference sequences with Bowtie2 (35). Especially in the case of nad genes, terminal modules were missing. They were recovered by mapping RNA-Seq reads onto partial transcript contigs, and subsequently by extending the sequences via RT-PCR up to the polyA tail using oligo-dT and gene-specific primers (see above). To screen the contigs across the investigated species for the highly divergent nad genes (the previously designated y genes (36), we employed HMMER 3.1b2, a most sensitive method based on profile hidden Markov models (37).

Chromosome classification
The module-containing contigs were first extended by read mapping using the mapper software implemented in the Nucleic Acids Research, 2020, Vol. 48, No. 5 2697 Geneious 10.1.3 package (34) and then compared with each other by BLASTn. Cassettes were identified by the same criteria as described previously (27). Briefly, cassettes are unique sequences surrounding modules, and are flanked by the constant regions of a chromosome, which we define as sequences with >90% identity over >100 bp adjacent to cassettes ( Figure 1B). Contigs bearing the same cassetteflanking regions were assigned to the same chromosome class. For multiple sequence alignments, we used MAFFT v7.388 (38). Chromosome classes were ordered from the highest to the lowest member count and named A, B etc. To calculate mean read coverage of cassettes or modules, DNA-Seq reads were first mapped onto mitochondrial contigs by the Geneious 10.1.3 software. Aligning reads were merged with BBMerge (rem k = 62 extend2 = 50 ecct) and deduplicated with Dedupe (k = 31 ac = f) from BBTools. The resulting reads were mapped with Bowtie2 onto the reference sequence, and mean coverage was calculated using the Pileup tool from BBTools. For N. karyoxenos mitochondrial chromosome sequences, which are highly polymorphic, DNA-Seq reads were mapped with the Geneious 10.1.3 software without subsequent Bowtie2 mapping. The final list of identified chromosomes has been compiled in Supplementary Table S4.

In silico identification of RNA editing and DNA polymorphic sites
RNA editing clusters and longer insertions were identified by comparing the genomic contigs and mature transcript sequences by BLASTn. To distinguish between RNA editing sites and genomic polymorphisms, RNA-Seq and DNA-Seq reads were mapped onto mitochondrial transcripts and genomic contigs, respectively, with the built-in aligner of the Geneious 10.1.3 software. DNA polymorphic sites were identified as those exhibiting two (or sometimes more) nucleotides in the mapped reads, while in the case of RNA editing sites, the consensus nucleotide in genomic reads differed from that in RNA-Seq reads. A genomic position with >10% reads displaying difference from the reference was considered a polymorphic site. For RNA editing, a position was annotated as an editing site if at least 50% reads carried a base change. Note that a vast majority of sites was edited to >90%. RNA editing and DNA polymorphic sites of each transcript are detailed in Supplementary Table S5.

Phylogenomic analysis
We used all 15 assigned mitochondrially encoded protein sequences (Atp6, Cob, Cox1/2/3 and Nad1/2/3/4/4L/5/6/7/8/9) from 11 diplonemids and the corresponding homologs from other discobans, namely Trypanosoma brucei, Bodo saltans and Perkinsela sp. (Kinetoplastida); Euglena gracilis (Euglenida), Acrasis kona, Naegleria gruberi and Stachyamoeba lipophora (Heterolobosea); Tsukubamonas globosa (Tsukubamonadida); and Andalucia godoyi, Reclinomonas americana (ATCC 50394) and Ophirina amphinema (Jakobida). Sequences were downloaded from the NCBI GenBank Protein database except those of B. saltans. We retrieved the latter through tblastn searches (using T. brucei mitochondrial proteins as queries) from transcripts that we assembled by Trinity v2.2.0 using RNA-Seq data deposited in the NBCI Bioproject PRJEB3146. Multiple sequence alignments (MSAs) of proteins were generated with MAFFT v7.38 (38) using the E-INS-i algorithm and default parameters. Protein alignments were stripped of hyper-variable sites (20% gap threshold) with trimAl v1.4.22 (39). Subsequently, protein sequences were concatenated for each species, with the final MSA containing 22 taxa and 5,063 positions. Phylogenetic inferences were performed by a Bayesian approach using posterior probabilities as support values (PhyloBayes v4.1 (40), MrBayes v3.2.6 (41)) and by maximum likelihood with bootstrapping (IQ-TREE v1.6.10 (42) and RAxML v8.2.11 (43)). Bayesian methods were executed in two independent chains and the first 25% cycles were discarded as burn-in. For Phylobayes, we chose the substitution model CAT-GTR and site rate variation modeled as a Dirichlet process (ratecat option); the chains were stopped after they converged (i.e. maxdiff below 0.1 at ∼750 cycles corresponding to ∼25,000 generations). For MrBayes, we chose the GTR model with six discrete categories of gamma rate variation and 200,000 MCMC generations. For ML computations, we chose the substitution matrix LG for amino acid frequencies, which was determined as the best model by Model Finder (44). For IQ-TREE, we used default parameters with the option to calculate 1,000 ultrafast bootstrap replicates. For RAxML, additional parameters were: 50 categories for rate heterogeneity (CAT option), the algorithm 'rapid bootstrap analysis' and 100 distinct alternative runs on distinct starting trees for bootstrap support values. To evaluate the reliability of the inferred tree, we further analyzed the gene-and site-concordance factors (gCF and sCF, respectively) for each branch, as implemented in IQ-TREE v1.7 (45), with default parameters and the option to merge models across loci.

Gene repertoire
We examined mitochondrion-encoded genes from four Diplonemidae species (Diplonema japonicum strain YPF1604, Rhynchopus humris YPF1608, Lacrimia lanifica YPF1601 and Sulcionema specki YPF1618) and two Hemistasiidae species (Artemidia motanka YPF1610 and Namystynia karyoxenos YPF1621). In all six species, we identified the same set of genes described earlier in four Diplonemidae species (27), namely genes encoding ATP synthase subunit 6 (atp6), cytochrome b (cob), three cytochrome c oxidase subunits (cox1, cox2 and cox3), 10 (36), as well as small and large subunit mitoribosomal RNAs (rns and rnl) (Figure 2; Supplementary  Tables S3 and S4). In the two hemistasiids, we failed to detect nad6 [y5], but this was presumably due to the gene's high divergence and not to its genuine absence. Moreover, Lacrimia, Sulcionema and Namystynia also encoded y4, a gene first discovered in D. papillatum; in R. humris, we found a candidate corresponding to module 2 from D. papillatum (y4-m2), but not m1. With homologs of the Nucleic Acids Research, 2020, Vol. 48, No. 5 2699 latter gene at hand, we revisited data from a previous study (27), which enabled us to detect the two y4 modules in R. euleeides, unrecognized previously because of extensive overlaps with cox3-m1 and nad5-m11 in this species. In Lacrimia, we found an additional gene, named y7 (single module), which potentially codes for a protein of 67 amino acid residues. No tRNAs were found; as in other euglenozoan species, they are apparently not mitochondrially encoded, but reside on nuclear DNA and are imported to the mitochondria. The inferred mitochondrion-encoded proteins of all species analyzed here and those studied previously (27) displayed an exceptionally low level of sequence conservation, which made detection of most genes challenging. Cox1 was the most conserved protein across the diplonemids (32.7% identity across 11 species), while Nad3 with a mere 2.5% sequence identity was on the other end of the spectrum (Supplementary Figure S1).

Module numbers and sizes
The four Diplonemidae species investigated here build their 17-19 identified mitochondrial genes from essentially the same number of modules as does the type species D. papillatum. The sole exception is Sulcionema rns, which is encoded by two modules instead of one in all other 'classical' diplonemids. In the two Hemistasiidae species, the total number of modules is doubled (Table 1 and Figure 2), while module sizes are halved ( Figure 3). In fact, a given module observed in Diplonemidae is typically split into two to four modules in Hemistasiidae, since gene breakpoints are typically conserved across diplonemids (for sequences that can be confidently aligned, occasional shifts are less than 6 bp).
About 4% of modules in hemistasiids are shorter than 20 bp--even as short as 3 bp--and referred to in the following as mini-modules ( Figure 4A and Supplementary Figure S2). It should be noted that mini-modules cannot be unambiguously distinguished from appendage RNA editing (see also below) by inspection of DNA-RNA sequence differences alone. Still, two lines of evidence support mini-modules. First, we confirmed by 3 RACE and subsequent sequencing of Artemidia cox3 the existence of an mRNA trans-splicing intermediate containing the putative cox3-m9 mini-module. The 3 RACE RT-PCR sampled two amplicon populations: a minor one with 5 Us appended to the 3 terminus of the upstream module 8 and a major one with the triplet CAG, corresponding to the mini-module 9, joined to the aforementioned U-tract ( Figure 4B). The triplet thus appeared to have been added as a whole, i.e. trans-spliced, rather than as a succession of unrecognized editing events. However, to completely rule out the latter alternative, a more extensive sampling of 3 -end RNA processing and editing intermediates by RNA-Seq would be necessary. Second, in the available RNA-Seq data, we detected RNA processing and trans-splicing intermediates, which contained in addition to the diminutive module its flanking sequence, thus indicating where in the genome it resided. This way, we could trace back the genomic source of four and eight such gene pieces in RNA-Seq reads of Artemidia and Namystynia, respectively. (Note, however, that for six additional short segments of 2-6 nt in five transcripts of Namystynia this was not pos-sible and the corresponding regions have been marked as unresolved [ Figure 2F].) For the putative cox3-m9 minimodule, we could detect two RNA-Seq reads containing the CAG triplet joined upstream to the cognate cox3-m8 (with the appended U-tract) and downstream to cox1-m21 ( Figure 4C). This indicated that cox3-m9 and cox1-m21 of Artemidia were actually juxtaposed on the same chromosome U01. The observation of the intermediates containing mini-modules with their flanking sequences led us to hypothesize a possible assembly scenario for mini-modules, where a larger precursor acts as a mini-module carrier (Figure 4D).

Classes of mitochondrial chromosomes
Knowing the module sequences from transcriptome data allowed us to identify the corresponding regions in genomic contigs, while the repetitive sequences adjacent to the unique module-flanking regions allowed us to delimit cassettes (see Figure 1). Further, cassette-flanking repetitive sequences were presumed to be part of the constant regions of chromosomes, according to the classification scheme of mitochondrial chromosomes in D. papillatum (16), and thus allowed categorization of chromosomes into multiple classes (Table 1;  For Namystynia chromosomes, in contrast to the other species, we could not employ the criterium of recurring cassette-flanking sequences (representing the constant regions), because sequences around modules frequently consisted of unequally spaced tandem and dispersed repeated homooligomeric motifs that could not be unambiguously aligned. However, numerous chromosomes shared motif 1 (5 -GGGCCAAAAA-3 ) upstream and motif 2 (5 -TTTTGGGCC-3 ) downstream of the cassettes. Consequently, all chromosomes bearing these motifs were classified as class X, which is much more diverse than the classes from the other diplonemids. Finally, in every species, a handful of chromosomes did not fit into a defined class and therefore were grouped into the category 'unclassified'. Total counts of classified and unclassified chromosomes for each species are summarized in Table 1 and Supplementary  Table S4.
Since the sequence repeats prevented the assembly of whole chromosomes from the available short reads (except for two cases in Sulcionema; see below), chromosome sizes remain unknown. Nevertheless, the assembled genomic contigs indicate that the overall chromosome architecture conformed to that previously observed in other diplonemids (24,27,28), namely cassette sizes varied from ∼0.2 to ∼2 kbp with median length ∼330 (±50) bp. Two types of deviations were observed. First, the median size of D. japonicum cassettes was ∼1 kbp; as we detail below, this was due to the unusually high number of modules per cassette. Second, Sulcionema and Artemidia chromosomes contained cassettes with sizes well above the 2 kbp mark (from 3.2 to 12.3 kbp). Based on the complete assembly of three Sulcionema class D chromosomes (Supplementary  Table S4), we calculated that these long cassettes covered   n.d., not determined. a atp6-m1 on two different chromosomes (counted once). b nad2-m2 on two different chromosomes (counted once). c nad3-m5 on three, and nad4-m3 and nad4-m5 on two different chromosomes (counted once). d Chromosomes assigned to a class based on different criteria than in other species.
83-90% of the circles, the complete opposite of the situation in other analyzed diplonemid chromosomes, where a cassette represents only 5-10% of the chromosome length (24,27,28).

Module content and arrangement
Diplonema papillatum has 81 distinct mitochondrial chromosomes, 76 of which carry a single cassette that in turn contains a single module (mono-module/mono-cassette organization). The remaining chromosomes contain one cassette each that encloses two modules (three instances; multimodule/mono-cassette organization) or cassettes without any identified module (two instances). In the other diplonemids examined previously and here, additional arrangements coexist, notably three or more (up to 11) modules per cassette and also two cassettes per chromosome (multi-module/multi-cassette organization) (27). In these Nucleic Acids Research, 2020, Vol. 48  latter instances, modules are either separated, overlapping or nested. The six species analyzed here differed considerably in their total number of distinct chromosomes, ranging from 30 in D. japonicum to 163 in Namystynia (Table 1 and Supplementary Table S4). These differences were not only due to a different number of modules in a given species, but also to the fact that some chromosomes encoded multiple modules. For example, D. japonicum contained 23 multi-module chromosomes, the highest number in this category among all species analyzed (27), but only 7 mono-module chromosomes (Table 1). In contrast, among the 86 chromosomes of Lacrimia, only four contained multiple modules. The highest number of modules detected in a single chromosome was 11 in Sulcionema; this species also contained by far the largest number of apparently module-less cassettes (all 16 from its A class chromosomes; Supplementary Table S4).

Conserved types of mitochondrial RNA editing and their distribution
In all species studied here, we identified mitochondrial transcripts that underwent multiple events of C-to-U and Ato-I substitution editing and U-appendage editing (Figures  2 and 6; Supplementary Table S5). These types of editing had also been described previously in four Diplonemidae species and H. phaeocysticola (16,27). In the type species, the presence of inosines in transcripts has been demonstrated experimentally indicating that post-transcriptional A-to-G DNA-RNA differences arose by deamination, a process that most certainly applies to C-to-U changes as well (16). Substitution editing sites occurred in clusters in similar gene regions across species, although individual sites did not necessarily coincide. Earlier reported substitution editing clusters (e. g. nad2-m4, nad3-m2, nad4-m1, nad6-m1,  nad7-m3 and m5, nad9-m1 and rns) were present in the species studied here as well, although with some exceptions, such as the complete absence of substitution editing sites in nad3, nad7 and nad9 of Lacrimia ( Figure 2). We also identified new sites of C-to-U and A-to-I substitutions and U appendage. These included one and two new C-to-U editing sites in Lacrimia cob-m5 and cox3-m1, respectively. Further, Sulcionema possessed a novel A-to-I substitution site in cox3-m3, and in R. humris, we discovered two new editing sites at the junction of m1 and m2 of cox3 (Supplementary  Table S5).
All the above editing types were much more frequent in the hemistasiids and affected every single transcript (Figures 2 and 6; Supplementary Table S5). In addition to the editing clusters documented previously (16,27), we identified several novel instances, mostly located at the ends of modules (which complicated the recognition of the corresponding modules). Although more numerous, substitution editing sites in Artemidia (493 A-to-I and 620 C-to-U sites in >100 editing clusters) were amassed in half the number of clusters compared to Namystynia (458 A-to-I and 588 C-to-U sites in >210 editing clusters; Supplementary Table S5). Interestingly, certain editing sites in Namystynia coincided with one of the >300 genomic polymorphisms (mostly single-nucleotide polymorphisms, SNPs) dispersed across modules (Supplementary Table S5 and Supplementary Figure 3A).

Novel types of RNA editing
By inspecting DNA-RNA differences, we detected two new types of editing not documented before in diplonemids. First, the two hemistasiids carried G-to-A substitutions ( Figure 7A; Supplementary Figure S3A,B), notably 14 unambiguous sites in six different transcripts of Namystynia and one such site in nad5 of Artemidia (Figures 2 and 6; Supplementary Table S5); a second site may exist in Artemidia nad1; however, the corresponding A in RNA could have also originated by A-appendage (see below).
The second new type of editing was detected in Sulcionema and Namystynia. The nad4 transcript of the former contained between m4 and m5 not only a non-encoded U but also an additional A (Figure 2), which we confirmed by 3 RACE (Supplementary Figure S3C). Such Aappendage editing appeared far more frequent in Namysty- The majority of >300 RNA-Seq mate 1 reads extending to or spanning the nad5-m12/m13 junction (not shown) map exactly to the sequence of the transcript containing nad5-m12, the RNA editing-appended tract and nad5-m13. A minor population shown here represents nad5-m12 modules (4.6%) containing a 3 terminal A-tract, which we interpreted as RNA editing intermediates. The apparent G-appendage in N. karyoxenos might thus originate from an A-appendage followed by A-to-I deamination.
nia, in which we spotted one to four such sites in eight different transcripts, summing up to a total of 17 editing positions and 41 post-transcriptionally added As (Figure 2 and Supplementary Table S5). The A-appendage site of Sulcionema nad4 also occurred in the same transcript of Namystynia (m13/m14 junction). Curiously, at all sites except the one in nad4L, A-appendage in Namystynia coincided with U-appendage, occasionally in an interspersed fashion, as for example between nad4-m1 and m2 (5 -UUAAUUUUUUUUU-3 ).
In Namystynia, we also observed six cases of apparent G-additions between modules, and again intermixed with other post-transcriptionally added nucleotides, for example 5 -GAUUU-3 between nad5-m12 and m13 (Supplementary Table S5). The Gs could have arisen in two ways, by (i) genuine G-appendage editing or (ii) A-appendage followed by A-to-I deamination. While editing by rare G-insertions had been observed in an amoebozoan and a heterolobosean (46,47), we considered the second alternative more likely because it would not imply an additional machinery required to specifically add G residues. Further, RNA-Seq read mapping to the junction of nad5-m12 and nad5-m13, where such a G-addition was noted, revealed a minor population of RNA-Seq reads that displayed 3 terminal A-tracts, which we interpret as not-yet edited intermediates ( Figure 7B). It will be interesting to validate this hypothesis by biochemical assays once Namystynia becomes more amenable to experimental work.

Phylogenetic relationships among diplonemids
In molecular phylogenies based on 18S rRNA (11,12,14) ( Figure 6), the genera Diplonema, Rhynchopus, Lacrimia, Flectonema and Sulcionema formed Diplonemidae, which are also referred to as 'classical' diplonemids, to the exclusion of Hemistasiidae, Eupelagonemidae (formerly DSPD I (48)), and the lineage currently described by its acronym 'DSPD II'. To build a more robust phylogeny, we used here the concatenated protein sequences inferred from 15 different mitochondrial transcripts, identified in this and previous studies. The resulting tree ( Figure 6, right tree; Supplementary Figure S4A) resolved the relationships between diplonemids with high confidence but differs in topology from the nuclear 18S rRNA-based trees including the same species ( Figure 6, left tree). First, in the mitochondrial phylogeny, Lacrimia grouped together with Diplonema, while it is placed at the base of the Diplonema + Rhynchopus clade in the 18S rRNA trees. Second, Namystynia, and not Artemidia, was the sister taxon of Hemistasia. Finally, the most significant deviation was the position of Sulcionema, as it branched together with Diplonemidae in the nuclear trees, but with Hemistasiidae in the mitochondrial tree.
To examine possible reasons for this incongruence, we calculated gene-and site-concordance factors for each branch in the tree (45) (Supplementary Figure S4B). Compared to the concatenated dataset, single-locus phylogenies showed comparably low support for the positions of Lacrimia and D. papillatum; however, the majority of informative sites in the concatenated dataset supported the tree topology. The conflicting positions within Hemistasiidae were mainly due to the limited data currently avail-able for H. phaeocysticola (i.e. four proteins instead of 15). More importantly, most single-protein phylogenies placed Sulcionema prior to the divergence of hemistasiids, i.e. the topology shown here (Supplementary Figure S4A), but the site-wise support for this topology was almost identical with that of the two other mutually exclusive topologies, in which Sulcionema formed a sister group to either all diplonemids or the Diplonemidae clade. Further taxon sampling of basal diplonemids should allow to resolve the described inconsistencies.

DISCUSSION
In the past, most studies were performed on the type species D. papillatum which, incidentally, has recently become genetically tractable (49). Together with three other classical diplonemids that had formerly been examined at the molecular level (11)(12)(13)(14), these species represent a tiny and ecologically restricted fraction of the highly diverse group (4,14). Here, we have considerably expanded the knowledge about the diplonemid mitochondrial genomes by examining six recently isolated species from both the Diplonemidae and Hemistasiidae clades, thus covering a substantial part of diplonemid diversity (11)(12)(13).

Diplonemids are record holders in mitochondrial genome content and organization
It is worth noting that these protists, which were neglected until recently, carry the largest amount of mtDNA documented in an organelle, which in D. papillatum even exceeds that of nuclear DNA (18). In comparison to human mtDNA, which typically constitutes only about 1% of total cellular DNA and consists of a single, circularmapping 16.5 kbp molecule encoding complete proteincoding genes (13 in human versus 16 in diplonemids), mtDNA in diplonemids is unprecedented in its magnitude, while its gene expression mode adds a supplementary layer of complexity. Most of the diplonemid mtDNA is noncoding, with the extensive constant regions apparently carrying only the origin of replication and transcription initiation signals. The baroque organization of the mitochondrial genome and transcriptome in diplonemids is partially met by the well-studied case of sister trypanosomatids, which tells us that sustaining this extravagancy must require an enormously complex cellular machinery. Moreover, observing such unusual features in the free-living diplonemids challenges the common view that extreme oddities are synonymous with a parasitic lifestyle.

The enigmatic y4 gene
All diplonemids analyzed in this and previous studies (27) encode the same set of mitochondrial genes ( Figure 6). The only gene with patchy distribution, encountered in half of the species, is y4 encoding a hypothetical protein, which is poorly conserved at the sequence level and for which no homolog was found outside diplonemids. The Y4 protein of D. papillatum was recently detected by mass spectrometry in a respirasome supercomplex and, therefore, might represent a novel diplonemid-specific subunit of one of the respiratory chain complexes (36). Alternatively, Y4 might specify a Nucleic Acids Research, 2020, Vol. 48, No. 5 2705 highly derived mitochondrion-encoded mitoribosomal protein. For example, the kinetoplastid RPS12 (encoding the mitoribosomal uS12m) and MURF5 are renowned for their extreme sequence divergence, and the latter has been uncovered as RPS3 (uS3m) only by structure determination of the Trypanosoma brucei mitoribosome (50). No homologs encoding uS12m and uS3m were detected in the D. papillatum nuclear genome (our unpublished observations), which suggests that Y4 may be an extremely divergent mitoribosomal protein. Structure determination of the Diplonema respirasome and mitoribosome will be the ultimate test of these hypotheses.

Mitochondrial gene fragmentation at new heights in hemistasiids
An earlier study of four genes indicated high fragmentation in Hemistasia mtDNA (29). Our more systematic investigation of complete mitochondrial transcriptomes from two other hemistasiids generalized this finding, uncovering putative mini-modules as short as 3 bp. The vast majority of mini-modules is embedded in another module, which indicates that increasing gene fragmentation facilitates double use of coding sequence for distinct genes. Importantly, reuse can be even multiple: in Namystynia, cox1-m15 and cox1-m6 are both embedded in cox1-m13, while in Artemidia, nad2-m1 and nad5-m20 extensively overlap with nad5-m13, with a 9-bp region contributing to all three gene pieces (Figure 5 and Supplementary Table S6).
Ultra-short (1-30 nt) coding sequences are found in nuclear and mitochondrial genomes of numerous organisms (51)(52)(53). These micro-exons are joined to their neighbors--by the spliceosome or the Group I or Group II splicing machineries--via cis-splicing, thus relying on a physical connection between exons for proper joining. The hemistasiid case suggests that mini-module joining might proceed through an intermediate where a larger precursor acts as a 'carrier' of the mini-module ( Figure 4C and D), ensuring the correct trans-splicing of a sequence that alone is presumably too short to ensure specificity. The actual mechanism of module transcript match-making still remains an intriguing puzzle.

RNA editing at an unprecedented level
Among diplonemids, the largest number of mitochondrial editing sites was counted in Namystynia, notably over 1,000 A-to-I and C-to-U substitutions, 14 G-to-A changes, and 94 U+A-tracts that sum up to >600 nt added to modules (Supplementary Table S5). As in previously studied diplonemids (16), RNA editing had an overall restorative effect on coding sequences, allowing production of functional proteins from a priori defective gene pieces.
Physarum polycephalum is also one of the few species known to employ more than one mode of mitochondrial RNA editing: co-transcriptional nucleotide insertions and occasional deletions (57), and post-transcriptional C-to-U substitutions (47,58), while diplonemids feature substitution (C-to-U, A-to-I and G-to-A) as well as U-and Aappendage editing.

New types of RNA editing
In hemistasiids, we detected two types of editing novel for diplonemid mitochondria, which involve G-to-A substitutions and A-appendage to internal modules. G-to-A editing has been only rarely reported in mitochondria, e.g. in dinoflagellates such as Hematodinium (59), while such events are extremely uncommon in the nucleus (60,61). Attesting to the importance of this type of editing in hemistasiids, the G-to-A substitution site in nad5-m6 of Artemidia is also conserved in Namystynia. The editing event contributes to the replacement of a Ser by an Asp codon that corresponds to the function-critical residue at position 179 in mammalian Nad5, an amino acid involved in the proton relay of complex I (62).
The molecular mechanism of G-to-A editing remains a matter of speculation; while C and U interconversion can proceed by transamination (U-to-C) and deamination (Cto-U)--since the two bases differ only in the absence or presence of an amino group--G and A differ in two groups, and no single chemical reaction is known to interconvert these two bases.
More concrete notions exist about A-appendage editing, which is a crucial step in the maturation of the dinoflagellate cox3 transcript (63,64). The reaction is presumably catalyzed by the poly(A) polymerase that otherwise adds poly-A tails to mitochondrial transcripts. The corresponding enzymes have been characterized in mammals and trypanosomes (65). In the latter, the 3 tails are actually a mix of A+U residues, generated in two steps. Prior to editing, which in trypanosomes involves U-insertions and Udeletions (66-68), a 20-25 residue-long 3 A-tail is added. Once editing is completed, this tail is elongated to a 200-300 nt-long A+U heteropolymer, earmarking the transcript for translation and allowing its association with the mitochondrial ribosome (69). Both short and long tails are synthesized by the kinetoplast poly(A) polymerase 1 (70), which forms a complex with two pentatricopeptide proteins called kinetoplast polyadenylation/uridylation factors (KPAFs) 1 and 2 (71). RNA-editing terminal uridylyl transferase 1, which forms a complex with 3 exonuclease (72), is involved in the formation of long mRNA 3 tails (73,74), and possibly also in uridylation of rRNAs and gRNAs.
We expect a similar protein complex to operate in diplonemid mitochondria. However, in diplonemids, the addition of As--frequently together with Us--takes place in two fundamentally different contexts: not only at the 3 end of terminal modules, thus generating mRNA tails, but also of internal modules, as we documented here for Namystynia (Figures 2 and 7A; Supplementary Table S5). It will be interesting to examine in diplonemids whether a single protein complex is responsible for generating both Atails and A-appendages or whether distinct specialized complexes have evolved for the two purposes.
Mitochondrial genes of Sulcionema--primitively simple or reduced upon divergence?
The mitochondrial system of Sulcionema appears in several aspects less complex than that of the other diplonemids. Furthermore, this species has the shortest branch in both mitochondrial ( Figure 6A and Supplementary Figure 4A) and nuclear phylogenies (11). Some of these features might have been reduced upon divergence, while others might be primitive.
Of particular interest is the absence of the posttranscriptionally added Us between modules m4 and m5 of Sulcionema cox1 (Supplementary Figure S5). Besides being the first editing site identified in a diplonemid (23), this Uappendage had apparently an important evolutionary impact on the Cox1 protein structure. In all species except diplonemids, the protein region corresponding to junction m4/m5 (loop 1) is positively charged and, in the folded protein, interacts with a downstream loop 2 composed of small and hydrophobic residues. The inverse situation applies to all diplonemids that feature cox1 U-appendage editing. Here, the U-tract added at the m4/m5 junction and its environs specify a hydrophobic patch, whereas the downstream loop contains a polar Arg residue (28). Interestingly, the Sulcionema Cox1 protein has exactly the same hydropathy pattern in loops 1 and 2 as the other diplonemids, only the Us at the m4/m5 junction are genome-encoded as part of m4. More extensive taxon sampling will be necessary to untangle the order of the two evolutionary events, loop 1/loop 2-polarity switching and U-appendage editing.
Other less complex features in Sulcionema include the lack of nested modules (Supplementary Table S6) and the low RNA editing frequency (Supplementary Table S7) compared to the other diplonemids. Furthermore, the range of copy numbers across its module-bearing chromosomes is only ∼13 (∼30 when including its 16 module-less chromosomes), but 50-150 in the other Diplonemidae, and even ∼600 in the hemistasiid Artemidia (up to ∼300 of its B class chromosomes alone) (Supplementary Table S4) (27). It was noted that unequal copy numbers of chromosomes in multipartite genomes may cause chromosome loss during random mtDNA segregation to daughter cells (75). One solution to this problem is to over-amplify mtDNA, which in D. papillatum represents >50% of total cellular DNA (18,75). It would thus be interesting to see whether the more even chromosome copy number distribution in Sulcionema correlates with a lower mtDNA to nuclear DNA ratio.

Gene fragmentation and complexity in the most diverse diplonemids
Given the very high estimate of diplonemid species in the ocean (5,7), it can be safely predicted that species with even more complex mitochondrial genomes and transcriptomes will eventually be discovered, especially among hemistasiids and eupelagonemids. For example, when reanalyzing the published genome sequences from 10 single-cell eupelagonemids (4), we detected in the data from 'cell 13' three potential mitochondrial modules encoding highly conserved regions of cox1 and nad7. One candidate module corresponds exactly to the hemistasiid nad7-m2, the second to the upstream half of the hemistasiid nad7-m3, and the third is a homolog of cox1-m11 from Hemistasia (Supplementary Figure S5). This indicates that the mitochondrial genomes of Eupelagonemidae species are similar to those of the Hemistasiidae clade with respect to module sizes and RNA editing (Figure 2 and Table 1; Supplementary Table  S5).

CONCLUSIONS AND OUTLOOK
Since diplonemids are a highly successful group as to their geographic distribution and habitat diversity, their extravagantly complex mitochondrial system has apparently little if any impact on their fitness. We see in it an excellent example of constructive neutral evolution (76)(77)(78)(79), which postulates that a stepwise increase in the complexity of a given cellular machinery can occur with no associated selective consequence. Regardless of whether the flexible gene module structure allows sequence reuse for unrelated genes and/or de novo generation of 'improved' gene pieces, the variety of splicing and editing events makes the mitochondrial genome of diplonemids a laboratory for new inventions. Indeed, in a diplonemid cell, it takes 'the whole village' to decode the handful of mitochondrial genes. Our next challenge is to identify the individual players involved in decoding and to establish how the complex gene expression is coordinated.