High molecular weight DNA was extracted from single individuals of the nemertodermatid N. westbladi stored in either ethanol, RNAlater, or RNA Shield using two different methods: the salting-out protocol (Guimaraes, 2018) and the QIAamp Micro DNA kit (Qiagen). The Qubit dsDNA HS kit, a 2% agarose gel, and a Femto Pulse system were used to ensure the extraction met the minimum requirements for DNA yield and fragment size (the majority of gDNA over 20 kb). These worms were isolated in dishes of seawater before their long-term preservation.

Library preparation and sequencing followed the PacBio Ultra-Low DNA Input protocol with small modifications. This protocol has been heavily tested at the Welcome Sanger Institute by Dr. Laumer, who kindly shared his experience and recommendations with us (Laumer, pers. comm.). Briefly, DNA was sheared to 10 kb using Megaruptor 3 instead of Covaris g-TUBE (Covaris). After removing single-strand overhangs and repairing the fragment ends (SMRTbell Express Template Prep Kit 2.0; PacBio), DNA fragments were ligated to the amplification adapter and PCR amplified in two independent reactions (Reaction Mix 5A and 5B) for 15 cycles each (SMRTbell® gDNA Sample Amplification Kit; PacBio). Amplified DNA was purified using ProNex Beads (Promega), pooled in a single sample, damage-repaired for the second time, and ligated to the hairpin adapters. Size selection of the prepared SMRTbell library was done using a 35% dilution of AMPure PB beads (PacBio), which removed all fragments shorter than 3 kb, instead of the BluePippin system (Sage Science). Finally, the library was sequenced in one SMRT cell on the Sequel IIe platform.

The ‘Trim gDNA Amplification Adapters’ pipeline from SMRT Link v11 was used to remove the adapter sequences. Three genome assembly strategies were attempted and compared: the IPA HiFi Genome Assembler included in SMRT Link v11 (PacBio), Hifiasm v.0.7 (Cheng et al., 2021)⁠, and Flye v.2.8.3 (Kolmogorov et al., 2019)⁠. Default parameters were used in the three approaches, but with the “--meta” option activated in Flye. Based on genome length and completeness (measured with BUSCO and the metazoa odb10 database; Supplementary Table S1), the Flye assembly was selected for downstream analyses, which included two additional scaffolding approaches. First, the two N. westbladi transcriptomes were mapped to the genome using HISAT2 v.2.0.5 (Kim et al., 2019)⁠ and fed to P_RNA_SCAFFOLDER (Zhu et al., 2018)⁠. Second, the genome of S. roscoffensis was used as a reference to map the assembled genome with RagTag v.2.0.1 (Alonge et al., 2022)⁠. Unfortunately, none of these attempts improved the genome contiguity any further. Finally, assembly redundancy was removed using the kmerDedup pipeline (https://github.com/xiekunwhy/kmerDedup). First, the Kmer spectra (k = 21) from the trimmed reads and the assembly were extracted and then merged with Jellyfish v.2.3.0 (Marçais and Kingsford, 2011). The list of kmers was converted to a fasta file and mapped to the assembly with Bowtie2 v.2.5.1 (Langmead and Salzberg, 2012), using the very-sensitive mode and the end-to-end approach, with a minimum alignment score of “L, −0.6, −0.2”. Then, kmerDedup was run with different minimum kmer coverage (mcv) and maximum duplication percentage (mpr), but the best combination according to BUSCO was to set both to 30%.

The raw assembly was decontaminated following the BlobTools2 pipeline (Challis et al., 2020)⁠. Coverage data was calculated by mapping the filtered HiFi reads to the assembled genome using Minimap2 with default parameters (Li, 2021)⁠, genome completeness inferred with BUSCO v.5.2.2 (Seppey et al., 2019) and the Metazoa odb10 database, and taxonomic information was identified through BLAST searches of the contigs versus the UniProt database (Release 2022_05) using diamond v.0.9.26.127 (Buchfink et al., 2021). Only the contigs identified as “Metazoa” were kept at this stage. Additionally, a BLAST search and a custom database were used to remove mitochondrial contigs. Finally, Minimap2 was used to map the reads back to the decontaminated genome to separate the nemertodermatid reads. The k-mer approaches GenomeScope v.2.0 (assuming a diploid genome and a 35x coverage) and SmudgePlot (Ranallo-Benavidez et al., 2020) were used to calculate the genome heterozygosity and ploidy before and after the decontamination step with a Kmer length of 21. To identify the contaminant contigs, the diamond output was used to extract the Taxid information of the hits, which is associated with a unique taxonomic category on the NCBI database.

RepeatMasker v.4.1.2-p1 (Smit et al., 2015)⁠ was used to soft mask the repeats in the decontaminated genome with the rmblast engine, for which a custom repeat database was generated with RepeatModeler v.2.0.1 (Smit and Hubley, 2015) and the -LTRStruct option activated. Afterwards, the genome was annotated with BRAKER2 (Brůna et al., 2021) using transcriptomic and proteomic evidence. The two available transcriptomes for N. westbladi, sequenced from a whole-body extraction, were downloaded and quality filtered in a two-step approach. Adapters removal and the first trimming were performed with Trimmomatic v.0.36 (as implemented in Trinity v2.6.6 with default parameters, (Grabherr et al., 2011)), followed by a more thorough cleaning with PRINSEQ v.0.20.3 (Schmieder and Edwards, 2011): trim all terminal bases with a quality below 30 and filter out reads whose mean quality is below 25, low complexity sequences (minimum entropy 50), and reads shorter than 75 bp. Clean reads were mapped to the soft-masked genome with STAR v.2.7.9 (Dobin et al., 2013) and the options “--sjdbOverhang 100 --genomeSAindexNbases 13 --genomeChrBinNbits 15” and “--chimSegmentMin 40 --twopassMode Basic”. For the proteomes, the gene models from the acoel P. naikaiensis (Arimoto et al., 2019)⁠, the BUSCO Metazoa odb10 database, and a custom set of single-copy orthogroups, inferred from published transcriptomes with OrthoFinder v.2.4.1 (Emms and Kelly, 2019), were concatenated and mapped to the N. westbladi genome using ProtHint v.2.6 (Brůna et al., 2020). The inferred gene models were functionally annotated by pfam_scan v.1.6 (Mistry et al., 2007) and the PFAM 31.0 database.

The quality of the decontaminated genome was assessed using QUAST v.5.2.0 (Gurevich et al., 2013) and the completeness of the genome and the annotation with BUSCO v.5.2.2 using the Metazoa and Eukaryota odb10 databases. Since all the metazoan contigs were kept during the decontamination step, two approaches were followed to ensure they belong to the nemertodermatid genome. First, a distance tree was inferred with FastMe v.2.1.5 (Lefort et al., 2015) based on a distance matrix calculated with Skmer (Sarmashghi et al., 2017), an alignment-free method designed to estimate genomic distances, over the N. westbladi genome and 18 metazoan genomes downloaded from GenBank (Supplementary Table S2). Second, a phylogenetic tree was inferred from these genomes except for three for which the annotated proteome was not available. Briefly, orthogroups were inferred with OrthoFinder v.2.4.1 (Emms and Kelly, 2019) and cleaned from paralogs with PhyloPyPruner v.1.2.3 (Thalén et al., 2021) using the “Largest Subtree” method, collapsing nodes with bootstrap support lower than 60, and pruning branches more than five times longer than the standard deviation of all branch lengths in the tree. Then, orthogroups were aligned with MAFFT v.7.475 using the L-INS-i algorithm (Katoh and Standley, 2013), cleaned from poorly aligned sites with BMGE v.1.12 (Criscuolo and Gribaldo, 2010), tested for stationarity and homogeneity (symmetry tests) with IQ-TREE2 v.2.1.3 (Minh et al., 2020), and concatenated with FASconCAT v.1.05 (Kück and Longo, 2014). Finally, a phylogenetic tree was inferred using coalescence [ASTRAL; (Zhang et al., 2017)] and concatenation by maximum likelihood with a site-specific heterogeneous model (assuming 20 amino acid categories, C20) with IQ-TREE v.1.6.12 (Nguyen et al., 2015).

All the genome metrics, including length, contiguity, number of genes, and completeness, among others, were compared to the acoel genomes from P. naikaiensis (Arimoto et al., 2019)⁠ and S. roscoffensis (Martinez et al., 2022)⁠, which were also tested for contaminants using BlobTools2, following the same pipeline and with the same filtering criteria. The genomes of H. miamia and X. bocki (Gehrke et al., 2019; Schiffer et al., 2023) were not considered because an annotation file with details of protein structure is not available for any of them. Additionally, a second N. westbladi genome sequenced in an Illumina HiSeq2500 platform was also included in the comparisons to estimate the improvement in genome quality with HiFi data relative to a short-read approach. Briefly, DNA was extracted from a pool of 12 individuals, collected in the same location at the same time, the sequencing library was prepared with a Rubicon kit, and the sequencing generated more than 385 million reads. The Illumina reads were assembled with SPAdes v.3.14.1 (Bankev et al., 2012), with four Kmer lengths (21, 33, 55, 75) and error correction activated. Finally, this genome was analysed with the same parameters as the HiFi genome to eliminate contamination contigs (but not redundancy), produce completeness stats, and annotate gene models.

To analyse the evolution of gene content in Acoelomorpha, the annotated genomes of 18 metazoans were compared, including N. westbladi (Nemertodermatida) and P. naikaiensis and S. symsagittifera (Acoela) as representatives of Acoelomorpha, eight protostome genomes, four deuterostomes, and three cnidarians as the outgroup to Bilateria (Supplementary Table S2). Redundancies in the gene models of all genomes were removed with CD-HIT (Fu et al., 2012), clustering all sequences more than 95% identical, and then functionally annotated with pfam_scan v.1.6 (Mistry et al., 2007) and the PFAM 31.0 database. The annotated proteins were clustered using OrthoFinder v.2.4.1 (Emms and Kelly, 2019) and used to calculate the number of genes specific to or shared among the four main clades of interest: Cnidaria, Acoelomorpha, Deuterostomia, and Protostomia. Genes present in at least one cnidarian and one bilaterian were considered to be shared across Metazoa, whereas genes present in at least two of Acoelomorpha, Deuterostomia, and Protostomia were considered to be shared across Bilateria. The proportion of “metazoan” and “bilaterian” genes absent from each of the three bilaterian clades was calculated based on these two datasets.

This analysis was based on the results of Gąsiorowski et al. (Gąsiorowski et al., 2021), who used spatial transcriptomics to identify the genes involved in the development of the ultrafiltration excretory system in several protostomes and one hemichordate species. All the protein sequences annotated in that study were downloaded from GenBank except Hunchback, as they found no evidence of this gene being involved in nephridiogenesis, for a total of three structural proteins: Nephrin, Kirrel, and ZO1; and six transcription factors: Eya, Lhx1/5, Osr, POU3, Sall, and Six1. These genes were annotated in the same genomes used to analyse gene content evolution through BLAST searches with diamond v0.9.26.127 (Buchfink et al., 2021). The correct identification of these genes was later confirmed through phylogenetic analyses with IQ-TREE v.1.6.12 (Nguyen et al., 2015) and manual BLAST searches on the NCBI web server. The identification of the Lhx1/5 and Six1 transcription factors was not always straightforward, as they are thoroughly mixed in the phylogenetic tree with many other gene variants and sometimes different isoform names were proposed in the BLAST searches for the same sequence, and thus they represent a mixture of isoforms of the same gene. All genes were later confirmed to be correctly annotated by mapping them to the same proteins annotated in the N. westbladi transcriptome. A custom R script was written to locate the filtered genes in the GFF files and extract three metrics related to gene architecture: protein length, number of exons per protein, and average exon length per gene. Unfortunately, the GFF annotation file was not available for all these genomes, so not all of them could be included in this analysis (Supplementary Table S2). To ameliorate the misleading effect of highly fragmented genes we filtered out all proteins shorter than half of the average protein length of the respective gene (a total of 10 proteins). To test if the observed differences in the three gene metrics were statistically significant, the Shapiro-Wilk’s method and the Barlett test were used to check if they follow a normal distribution and the homogeneity of their variances, respectively. For each gene, the differences among clades were tested with either an ANOVA or a Kruskal-Wallis test, depending on the result of the normality and homoscedasticity tests. Finally, the Bonferroni correction (ANOVA) and the Dunn test (Kruskal-Wallis) were selected to run pairwise comparisons in all cases identified as statistically different. In all cases, a p-value of 0.05 was set as the significance threshold.

The best extraction was produced from a sample stored in RNAlater using the QIAamp Micro kit, obtaining a fragment size over 20 kb and ca. 20 ng of total DNA, which yielded 990 ng after DNA shearing and whole-genome amplification. About half of this DNA was selected for sequencing. A total of 2,313,071 reads were produced during HiFi sequencing, later reduced to 2,297,478 after quality filtering with an average length of 6.6 kb.

Flye produced the best assembly (Supplementary Table S1), which was 678.9 Mb long and contained 26,880 contigs, and later reduced to 458 Mb and 11,625 contigs when assembly redundancies were removed (Figure 1A). The longest contig was 2 Mb long, with an N50 of 57.1 kb and contained 82.6% of the BUSCO Metazoa odb10. The assembly contained two repeats of 507 and 531 bp with 70,000 and 79,000 copies, respectively, corresponding to 11% of the assembled (raw) genome. BlobTools2 revealed the presence of many contaminants, with only 78.9% of the contigs identified as metazoan (Supplementary Table S3). Thus, the decontaminated assembly was only 407.7 Mb, split into 9,167 contigs with an N50 of 60 kb (Figure 1A; Table 1), but 80.2% of the Metazoa and 88.6% of the Eukaryota BUSCO genes were still present (Supplementary Figure S1). The two databases reported a relatively high proportion of duplicated BUSCO genes, 7.8% and 11% in the Metazoa and Eukaryota databases, respectively (Supplementary Figure S1). The smudgeplot was markedly different before and after the decontamination step, as the inferred ploidy went from triploid to diploid after the decontamination (Supplementary Figure S2). The genome size estimated by GenomeScope was 237.1 Mb, with an average coverage of 39.6, and high heterozygosity (3.99%), although these numbers must be taken cautiously given the poor fit of the model (39.3%; Supplementary Figure S3).

The decontaminated Illumina genome was also relatively complete, with 76.8% of the metazoan BUSCO genes present in the assembly (with barely any duplicates, 0.6%), but much shorter (62.2 Mb) and much more fragmented (49,310 contigs; N50: 4 kb) (Figure 1A). Despite being sequenced from cultured, starved and symbiont-free populations, BlobTools also identified some contaminants in the published genomes of P. naikaiensis and S. roscoffensis. The former went from 656.1 Mb and 12,525 contigs to 581.4 Mb and 7,104 contigs, whereas the latter went from 1103 Mb to 3,460 contigs to 1064.9 Mb and 2,730 contigs (Figure 1A). The N50 of the two genomes rose from 127 to 130 kb in P. naikaiensis, and from 1.04 to 1.08 Mb in S. roscoffensis. Despite the observed differences in genome size and contiguity, the four genomes show very similar completeness results. Around 90% of the Eukaryota BUSCO genes were identified in the decontaminated genomes of all species, with P. naikaiensis presenting the fewest genes (14.9% of missing genes) (Supplementary Figure S1A). Differences among genomes were slightly higher with the Metazoa database, with almost a 10% difference between the most (S. roscoffensis; 18.5% missing genes) and the least (P. naikaiensis; 27%) complete genomes. In N. westbladi, the HiFi genome was almost as complete as S. roscoffensis (19.8% missing genes), whereas the Illumina genome was in an intermediate position (23.1%) (Supplementary Figure S1B).

The number of gene models in the four genomes ranged from 20,303 (P. naikaiensis) to 28,513 (S, roscoffensis), although the differences were reduced when only functionally annotated genes were considered: 10,074 (N. westbladi, HiFi), 13,708 (P. naikaiensis), 14,486 (N. westbladi, Illumina), and 17,717 (S. roscoffensis) (Table 1). The organisation of these genes in the genome somehow reflected the differences observed in genome contiguity. In the N. westbladi genome sequenced with Illumina, the average number of genes per contig was just 0.876, with a single gene in almost 90% of the contigs (Figure 1B), and the contig with the highest number of genes presented 33 gene models (Table 1). In the HiFi-sequenced N. westbladi genome, up to 92 genes were found in a single contig, with an average of 2.4 genes per contig. Similarly, an average of 2.9 genes per contig were annotated in the P. naikaiensis genome, but in this case, the maximum number of genes in one contig was only 37. The S. roscoffensis genome stands out, with a maximum of 280 genes in a single contig and more than 10 genes in almost 40% of the contigs (Figure 1B; Table 1). This trend, however, was not observed in gene architecture. The gene models in P. naikaiensis, S. roscoffensis, and the HiFi genome of N. westbladi were similar, ranging between an average of 3.2–6.3 exons per gene, whereas almost all the genes presented a single exon in the Illumina genome (average 1.5) (Figure 1D). One gene with over 1200 exons was annotated in the N. westbladi genome. This is likely an annotation error and it was not included in the calculation of these metrics, but it is shown in Figure 1D. The intron size was very variable in all genomes, ranging from 6 (P. naikaiensis) to 193,733 (S. roscoffensis) bp. The intron size distribution was similar between N. westbladi and P. naikaiensis, but with generally longer introns in S. roscoffensis (Figure 1C). Nevertheless, the intron size range was similar in the three genomes, but visibly smaller in the N. westbladi Illumina genome.

According to RepeatMasker, the N. westbladi genome is very repetitive, masking up to 59.66% of the genome (Supplementary Table S4). The majority of these repeats are interspersed throughout the genome (57.84%) and more than a fifth (22.61%) were not classified into any known repeat family. Among the classified repeats, the most common ones are retroelements (31.19%), particularly the long terminal repeats (LTR, 19.32%) and long interspersed nuclear elements (LINEs, 11.53%). The Illumina genome presents a sharp contrast, with just 16.40% of the genome masked as repetitive, although LINEs (4.28%) and LTR (3.43%) are still the most abundant repeat elements (Supplementary Table S4).

More than half of the taxonomic groups identified within the set of contaminant contigs were bacteria, including several of the major taxonomic groups: Bacteroidetes, Tectomicrobia, Proteobacteria (including Alpha-, Beta-, Delta/Epsilon-, and Gammaproteobacteria), Planctomycetes, Actinobacteria, Cyanobacteria, and Firmicutes. None of the “Candidate Phyla Radiation” phyla (Brown et al., 2015) were identified. More specifically, there are nine genera that have been reported as statistically more abundant in the microbiome of microscopic animals than in environmental samples (Boscaro et al., 2022) and thus might be part of the Nemertoderma microbiome: Algoriphagus, Alteromonas, Francisella, Photobacterium, Roseobacter, Shewanella, and Vibrio. Other important sources of contamination besides bacteria are algae (Chlorophyta, Rhodophyta, and Streptophyta), and fungi (Ascomycota, Basidiomycota, Microsporidia, Mucoromycota, and Zoopagomycota). These groups accumulate 87% of the taxonomic diversity within the contaminants. In addition, we also found Protista (Amoebozoa, Perkinsozoa, Endomyxa, and Oomycota), and Virus (Uroviricota). A complete description of these results is provided in Supplementary Table S5.

The comparison of 18 animal genomes, representing Acoelomorpha, Cnidaria, Deuterostomia, and Protostomia revealed a high degree of specificity in gene content: 17.2% of all orthogroups present in Cnidaria are exclusive to this phylum (1754 out of 10,172), 23.1% in Acoelomorpha (2,101 out of 9,080), 45.3% in Deuterostomia (7,976 out of 17,593), and 48.7% in Protostomia (9,573 out of 19,669; Figure 2A). Hence, only 33.4% (10,736 out of 32,141) of all orthogroups were annotated in at least two of the four groups. Among these, more than half (53.4%) were present in at least one species of each clade, whereas only 3.8% were present in all bilaterian clades but Cnidaria. A total of 8,418 genes were identified as shared across Metazoa (present in Cnidaria and at least one Bilateria), and 2,318 for Bilateria (present in at least two bilaterian clades). Acoelomorpha had 71.5% of the metazoan genes and 41.5% of the bilaterian ones, contrasting with deuterostomes (91.3% and 83.1%) and protostomes (94.4% and 92.9%) (Figure 2C). The proportion of missing BUSCO genes was below 11% in all four groups (Figure 2B), and so genome completeness does not explain this pattern. Within Acoelomorpha, almost half (43.2%) of the genes were shared between Acoela and Nemertodermatida (Figure 2A), 41.5% were unique to Acoela, and 15.3% to Nemertodermatida.

Nine genes, selected because of their participation in the development of the ultrafiltration excretory system, were investigated: three structural proteins (Nephrin, Kirrel, and ZO1), and six transcription factors (Eya, Lhx1/5, Osr, POU3, Sall, and Six1). All of them were annotated in both protostomes and deuterostomes. In Acoelomorpha, all genes but Osr were annotated, whereas only three out of the nine genes were found in the two cnidarian species (ZO1, Six, and Lhx; Figure 3A). According to GenBank, three more genes (Nephrin, Eya, and POU3) are also present in this phylum (Figure 3A).

The gene architecture (in terms of protein length, number of exons per gene, and average exon length) was compared for the nine genes among four clades: Cnidaria, Acoelomorpha, Deuterostomia, and Protostomia. Almost half of the 27 comparisons returned statistically significant differences among clades, most of them related to acoelomorphs (Figure 3B). Despite the evident variation in protein length, both within and among clades, only three out of the nine genes were considered to be statistically significant: Kirrel, which is significantly longer in acoelomorphs; ZO1, longer in deuterostomes; and Lhx, but in this case the differences were only significant between acoelomorphs (longer) and protostomes (shorter). As for the number of exons per gene, ZO1 and Eya presented fewer exons in acoelomorphs than in both deuterostomes and protostomes. Finally, the last gene with a significantly different number of exons is POU3. This is a relatively short protein, on average shorter than 500 amino acids in all clades, and with very few exons: only one exon in all deuterostomes but Branchiostoma floridae (three), between one and three in protostomes, and between one and four in acoelomorphs. Only the differences between deuterostomes and acoelomorphs were statistically significant. Two remarkable outliers were found when comparing the number of exons per gene. Three chordate ZO1 sequences were divided into more than 80 exons (average 29.5) and one of the POU3 sequences annotated in P. naikaiensis presented 15 exons (average in Acoelomorpha: 2.6). Nonetheless, these proteins were roughly of the same size as the others and their identity to the most similar protein was above 90%.

In an attempt to avoid the misleading effect of errors in the annotation (partial proteins will be generally shorter and with fewer exons), the average exon length was also considered. In this case, six out of the nine proteins were significantly different among clades. The average exon length was significantly longer in acoelomorphs in three genes (Kirrel, Eya, and Lhx), and two in deuterostomes (Sall and Osr, although the latter was only present in deuterostomes and protostomes). The only instance with significantly shorter exon lengths is the protostome’s ZO1 gene. Finally, among the nine comparisons including at least one cnidarian species (three genes, three metrics) no significant differences were found but in the average exon length of Lhx, which is significantly shorter than that of acoelomorphs, as also observed in deuterostomes and protostomes.

The steady development of sequencing technologies is allowing the generation of genomes spanning the diversity of life, which now includes minute organisms. Indeed, thanks to the latest low and ultra-low DNA input protocols, sequencing high-quality genomes from millimetric animals is now possible (Yoshida et al., 2018; Schneider et al., 2021; Lord et al., 2023). In this study, we used the Pacbio Ultra-Low DNA Input protocol to sequence the genome of N. westbladi, reporting the first nemertodermatid genome, sequenced from a single microscopic worm. The estimated genome length is shorter than in any sequenced acoel, and considerably shorter than S. roscoffensis and H. miamia (Gehrke et al., 2019)⁠. Although the P. naikaiensis genome is slightly more contiguous than N. westbladi, all the metrics compared are similar between the two genomes. In contrast, both S. roscoffensis and H. miamia were scaffolded using proximity ligation data, and hence both show much higher contiguity. Beyond the differences in contiguity, annotation metrics are comparable among N. westbladi, P. naikaiensis, and S. roscoffensis. In this case, N. westbladi is more similar to S. roscoffensis than to P. naikaiensis, which shows the lowest genome completeness and number of gene models. In particular, the analysis of gene architecture shows that the number of exons per gene and intron size is also comparable, likely meaning that the annotated proteins are complete or nearly complete, facilitating the study of gene properties, such as intron-exon structure. Likewise, all genomes are similarly repetitive: N. westbladi 59.66%; P. naikaiensis 69.8%; S. roscoffensis 61.14%; and H. miamia 53%, but this is where the difference between the short- and long-read genomes of N. westbladi strikes the most. Although they have similar completeness and number of gene models, the Illumina genome is only 62.2 Mb long and only 16.4% repeats, which is probably explained by the difficulty to assemble repetitive areas of the genome (Tørresen et al., 2019).

It is obvious from the comparisons above that achieving a highly contiguous genome from single-millimetre worms is still challenging. One potential explanation for this is genome size. Although the N. westbladi genome is shorter than the maximum genome size advised by PacBio, due to the contaminants the raw assembly is still almost 700 Mb long. The Ultra-Low DNA Input protocol has insofar been tested in animals whose genome size ranges between 200 and 300 Mb, returning significantly more contiguous genomes than that of N. westbladi (Kingan et al., 2019; Korlach, 2020; Schneider et al., 2021)⁠. The generally lower coverage of the raw assembly of the nemertodermatid genome, due to its larger size, could have resulted in a more fragmented assembly. Yet, sequencing a second HiFi SMRT cell was not feasible due to the low DNA yield. This protocol has been successfully tested on C. elegans and freshwater flatworms, demonstrating its good performance (Laumer, pers. comm.). However, it has never been used for sequencing >300 Mb genomes. Thus, it would be interesting to sequence the genome from a novel species with a similar genome size to confirm low coverage, and not unexpected artefacts, explain the fragmentation observed. Alternatively, one straightforward solution to improve genome contiguity is complementing this approach with ligation data, which has shown great results both in S. roscoffensis and H. miamia (Gehrke et al., 2019; Martinez et al., 2022)⁠. However, this approach would require pooling tens of individuals to obtain the required amount of DNA, which is not feasible for all animals. N. westbladi cannot be cultured in the lab and collecting worms in enough numbers is challenging. Interestingly, the P. naikaiensis genome (the most similar to N. westbladi) was sequenced from a pool of individuals in 52 SMRT Cells (Arimoto et al., 2019)⁠, whereas the N. westbladi genome comes from a single worm and one HiFi SMRT Cell. Altogether, these results highlight the potential of combining this protocol and HiFi reads to generate good-quality genomes from single, microscopic organisms, even for relatively large genomes.

The BlobTools analysis identified a high degree of contamination in the raw assembly of N. westbladi, which is to be expected from a microscopic organism caught in the wild. Although N. westbladi is known to not carry internal symbionts (based on hundreds of observations), a TEM analysis revealed the presence of gram-negative bacteria throughout the epidermal cilia (Lundin, 1998). All bacteria identified during the decontamination are gram-negative, matching this observation. Thus far, DNA extraction was performed from a whole specimen, thus sequencing the gut microbiome, and other contaminants might have been transferred from the DNA suspended in the seawater. A common practice to limit the presence of contaminants in the organism, and that was applied to both P. naikaiensis and S. roscoffensis, but not N. westbladi, is to starve the animals before DNA extraction. Besides, the two acoel genomes were sequenced from juveniles, before they incorporate the symbiotic algae, and rinsed with filtered seawater (Arimoto et al., 2019; Martinez et al., 2022). However, as seen here this is not enough to prevent the presence of contaminants. This was particularly problematic in the case of P. naikaiensis, as almost 4% of the contigs (75 Mb, over 10% of the genome) were identified as bacterial contigs. It is important to notice that a big fraction of the genomes did not have any hit against the Uniprot database (N. westbladi 18.1%; P. naikaiensis 8.4%; S. roscoffensis 1.9%; Supplementary Table S3), showing the importance of sequencing underrepresented groups to improve the reference databases.

The increasing availability of animal genomes has unveiled a remarkable diversity in genome sizes, ranging from 15.3 Mb in the orthonectid Intoshia variabilis to the 43 Gb of the lungfish genome (Slyusarev et al., 2020; Meyer et al., 2021)⁠. It has been observed that miniaturised animals tend to have smaller genomes, which has been noted both in vertebrates and invertebrates (Liu et al., 2012; Decena-Segarra et al., 2020; Xu et al., 2021), but with notable exceptions to this rule, as observed in nematodes and platyhelminths (Consortium, 2019). Genome length in the latter ranges between 700 and 1200 Mb, the same size range as birds, some gastropods, and many freshwater fish, among others (Zhang et al., 2014; Nam et al., 2017; Yuan et al., 2018). Similarly, acoelomorph genomes vary between 407 (the genome of N. westbladi is the shortest reported for Acoelomorpha) and 1059 Mb but contrast with the chromosome-level genome of X. bocki, estimated at 110 Mb (Schiffer et al., 2023). Comparisons of eukaryotic genomes proposed that variations in genome sizes and proportion of repeat elements are correlated (Elliott and Gregory, 2015; Shah et al., 2020), which might also apply within Xenacoelomorpha. Acoelomorph genomes show a much higher repeat content than the small genome of Xenoturbella (Schiffer et al., 2023).

In turn, acoelomorph genomes seem to be characterised by an important reduction of gene content. Indeed, almost 60% of the genes shared between protostomes and deuterostomes are missing in acoelomorphs, which could be explained by the morphological simplicity of these worms compared with other bilaterians, but the evolutionary interpretation depends on the phylogenetic hypothesis. Under the Xenambulacraria hypothesis, their absence must be explained by massive secondary losses. The Nephrozoa hypothesis, on the other hand, suggests that the evolution of the genes exclusively shared by deuterostomes and protostomes occurred in the stem line of Nephrozoa and no ad hoc hypotheses of gene loss are required.

Despite the absence of a specialised excretory system in Xenacoelomorpha, Andrikou et al. (2019) detected active excretion in this group through the digestive tissue and annotated several genes known to participate in the excretory mechanisms of nephrozoan animals. Here, we annotated in the genomes of Acoela and Nemertodermatida seven of the nine genes involved in the development of the nephridia and one more (Sall) in Acoela. Regardless of their phylogenetic position, whether as a sister to Ambulacraria or Nephrozoa, the presence of these genes might be explained by their participation in other important functions. A spatial transcriptomics analysis in the acoel Isodiametra pulchra and the nemertodermatid Meara stichopi located the expression of Nephrin in the brain and the nerve cords (Andrikou et al., 2019), which resembles observations in mammals and Drosophila, the latter through the Nephrin homolog Sns (Putaala et al., 2000; Putaala et al., 2001; Bali et al., 2022). In contrast, no homologs to the Osr gene (named Odd in Drosophila) could be annotated in any of the acoelomorph genomes. A BLAST search over the two Xenoturbella transcriptomes failed to annotate this gene in these species, confirming its absence is a general trait of Xenacoelomorpha. This is noteworthy, as Osr is essential in the formation of the excretory organs: in vertebrates, it participates in the formation of the pronephros, the first stage in kidney formation, and its knock-out results in the absence of kidneys (James et al., 2006); whereas in Drosophila, Odd participates in the embryogenesis of the tubules of Malpigi (Tena et al., 2007). Overall, it seems that the molecular machinery that participates in the functioning of a complex ultrafiltration excretory system is present in acoelomorphs, but they lack the one gene necessary to promote the formation of discrete excretory organs.

This pattern fits well within the Nephrozoa hypothesis. In this scenario, the origin of the molecular machinery associated with the excretory organs would be the result of gene co-option, a common phenomenon in the origin of key innovations, such as the development of the radula and shell evolution in molluscs (Hilgers et al., 2018) or the multiple origins of cnidarian eyes (Picciani et al., 2018). Interestingly, six of the nine genes investigated have been annotated in different cnidarian species, strengthening the idea of these genes pre-dating the appearance of this specialised excretory system (Gąsiorowski et al., 2021). Thus far, Osr has not been annotated in any phylum outside of Nephrozoa, supporting the origin of this gene in the ancestor of this clade. Nevertheless, given the ongoing debate around the phylogenetic position of xenacoelomorphs, the Xenambulacraria hypothesis also needs to be taken into consideration. If Xenacoelomorpha is the sister group of Ambulacraria, additional ad hoc hypotheses have to be invoked: either the Osr gene was independently gained in Protostomia, Ambulacraria, and Chordata or it was lost in Xenacoelomorpha. The Drosophila Odd gene has been shown to activate the formation of kidney tissue in vertebrates (Tena et al., 2007), which suggests a common origin of both genes in protostomes and deuterostomes. Likewise, the function of this gene is not limited to the development of the excretory organs, but it participates in the development of the foregut in vertebrates (Han et al., 2017) and it is known to be expressed in the digestive tract of spiralians and hemichordates (Gąsiorowski et al., 2021). Although its general anatomy is variable, the presence of a sack-like gut is considered a plesiomorphy within Xenacoelomorpha (Gavilán et al., 2019) and the involvement of Osr in its development could be expected. In this light, the reduction of the excretory organs alone would not explain the secondary loss of Osr, as it would need to be completely nonfunctionalized before that. Regardless of which of the two phylogenetic hypotheses is true, the acquisition of Osr and the development of a discrete excretory system seem to be interconnected. The two are a nephrozoan novelty and their origin likely dates back to the most recent common ancestor of protostomes and deuterostomes. More importantly, knock-out experiments demonstrate that there are no other current transcription factors capable of replacing Osr in the early formation of the pronephros (James et al., 2006), which suggests this gene must have assumed this function rather early.

We found statistically significant differences in the gene architecture of all genes but Nephrin and Six, six of them related to the average exon length. Acoelomorpha is responsible for two-thirds of the differences observed, which fits with the co-option of these genes into the development of the excretory system in the ancestor of Nephrozoa. Changes in gene structure are a strong generator of diversity, particularly after gene duplication, as part of the neofunctionalization of proteins (Xu et al., 2012). Alternatively, the differences observed might simply be explained by changes in the selective pressures during the acquisition or the reduction of this system, something that might be supported by the observations in Bryozoa. Within protostomes, Bryozoa, which also lack an excretory system, is responsible for most of the variation observed. Notably, half of the gene metrics that are visibly different in this phylum are shared with acoelomorphs: ZO1 and Lhx length, ZO1 number of exons, and Sall average exon length. However, the variation does not always go in the same direction (e.g., the number of exons in ZO1 increases in Acoelomorpha, but decreases in Bryozoa), likely because the absence of the excretory organs in the two animal groups represents two independent evolutionary events. Some authors have argued that the rapid evolutionary rates observed in Acoelomorpha might be associated with other traits observed in this group, such as chromosomic rearrangements or changes in gene content, misleading comparative analyses and making Xenoturbella a better model for studying the evolution of Xenacoelomorpha (Philippe et al., 2019; Schiffer et al., 2023). Unfortunately, the genomic data of X. bocki is yet not available so we have inferred a gene tree for each of the nine genes analysed and compared the differences in branch lengths among clades to explore this possibility (Supplementary Figure S4). Although branch lengths are indeed significantly longer in acoelomorphs than in any other clade (except in Lhx and Six), they are also longer in deuterostomes compared to protostomes despite the similarities between the two clades. In more detail, protostomes present the shortest branches in the gene trees, while Bryozoa is one of the phyla with the most changes in gene architecture. Hence, the accelerated evolutionary rates of Acoelomorpha do not seem to be the main factor underlying the differences observed in these genes, although it would be interesting to confirm this once all the data from the Xenoturbella genome is publicly available.