The complete mitochondrial DNA of the carnivorous sponge Lycopodina hypogea is putatively complemented by microDNAs

Sampling, DNA extraction, library preparation and shallow shotgun sequencing

An individual of L. hypogea was taken from an aquarium at the Station Marine d’Endoume (UMR-CNRS), Marseille, France (see Vacelet et al., 2022), immediately preserved in RNAlater^® (Thermo Fisher Scientific, Waltham, MA, USA), and kept refrigerated at −20 °C in the laboratory (and under room temperature during transit). In Brazil, the sample was deposited at the Museu Nacional (Universidade Federal do Rio de Janeiro, UFRJ), under the voucher number MNRJ 22102, and at the Laboratório de Genética Marinha (Universidade do Estado do Rio de Janeiro, UERJ), a fragment was used for total genomic DNA purification (BioSample ID SAMN40943203), following the protocol published by Salgado et al. (2007) using CTAB with minor modifications. Briefly, the sample was incubated in lysis buffer (CTAB 2%, NaCl 1.4 M, EDTA 20 mM, Tris-HCl 100 mM, pH 8.0, 2-Mercaptoethanol 0.2%, Proteinase K 50 $\mathrm{\mu }\mathrm{g}$/mL), following extraction with chloroform, precipitation with isopropyl alcohol, and resuspension in RNAse A (20 $\mathrm{\mu }\mathrm{g}$/mL) solution. The purified DNA was quantified through a fluorometric method using Qubit (Thermo Fisher Scientific, Waltham, MA, USA) and stored at –20 °C. The sample was sent for library preparation and sequencing to the Rush University Genomics and Microbiome Core Facility (GMCF, Chicago, IL, USA) through a third-party service. Library preparation was conducted using the Illumina DNA Prep kit (Illumina, San Diego, CA, USA; previously known as Nextera DNA Flex Library Prep) and 100 ng of the sample, and 2 $\times$ 150 bp paired-end sequencing was performed in a Illumina NovaSeq6000 system according to the workflow indicated by the manufacturer, using roughly 0.4% of a SP flow cell. Shallow shotgun metagenomic sequencing (SSMS) data from L. hypogea is deposited under BioProject ID PRJNA1099585.

Quality control, mapping, assembly, annotation and analysis

Low-quality reads (SLIDINGWINDOW:4:20), sequencing adapters (ILLUMINACLIP:2:30:10), and short reads (MINLEN:30) were removed from the raw data using Trimmomatic v0.39 (Bolger, Lohse & Usadel, 2014). The BBTools package v37.62 (Brian Bushnell, available at https://sourceforge.net/projects/bbmap/) was used to trim 5 bp at the start and the end of each read using the ‘bbduk’ script, to avoid contamination with any remaining adapter sequence. The reads were then mapped against 23 Heteroscleromorpha (Demospongiae) mitogenomes (Table 1) using the ‘bbmap’ script, with default settings. Unverified and unannotated mitogenomes available at Genbank were disregarded. The mapped reads were de novo assembled using the software Megahit v1.2.9 (Li et al., 2015) with default parameters. Protein coding genes were predicted using ORF Finder (https://www.ncbi.nlm.nih.gov/orffinder/) translated with minimally derived sponge mitochondria genetic code (the same as mold, protozoan, and coelenterate mitochondrial code) and alternative initiation codons other than “ATG”, and ontology inferred using BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Transfer RNA genes were predicted using tRNAscan-SE v2.0 (Chan et al., 2021), also assuming the genetic code as before. The boundaries of protein coding and ribosomal RNA genes were determined by BLAST and multiple sequence alignments (MSA) with genes retrieved from poecilosclerid mitogenomes using MAFFT v7.505 (Katoh, Rozewicki & Yamada, 2019). The assembled mitogenome was annotated manually and verified through the MITOS2 WebServer (Donath et al., 2019) using 89 reference sequences from opisthokonts and the mold/protozoan mitochondrial genetic code. Composition skew values were computed according to Lobry (1996), as follows: AT skew = $(A-T)/(A+T)$; and GC skew = $(G-C)/(G+C)$. The base and amino acid compositions and codon usage were obtained using SMS v2 (Stothard, 2000). Detection of terminal direct repeats (TDRs), terminal inverted repeats (TIRs), and miniature inverted-repeat transposable elements (MITEs) longer than 10 bp in IGR sequences was conducted using Generic Repeat Finder (GRF; Shi & Liang, 2019). Genes putatively missing from L. hypogea mitochondrial genome were individually mapped against the FASTQ reads and assembled into contigs using the same protocol above. Positive results were searched among the mitochondrial intergenic regions (IGRs) and the metagenome of L. hypogea, which was de novo assembled with Megahit using filtered reads. Since all positive hits were amongst the highly coverage metagenomic contigs, all contigs over 100 $\times$ depth had their homology checked through BLAST searches against the NCBI NT database. In addition, possible ORFs in IGRs and metagenomic contigs (all sequences between STOP codons over 30 nt) were extracted using EMBOSS’s ‘getorf’ (Rice, Longden & Bleasby, 2000) and their orthology were determined using eggNOG-Mapper (Huerta-Cepas et al., 2017). Metagenome analyses were conducted on a Galaxy server (The Galaxy Community, 2022; https://usegalaxy.eu/). The gene map of the mitochondrial genome of L. hypogea was generated using CGView v2.0.3 (Stothard & Wishart, 2005).

Phylogenetic analyses

The phylogenetic analyses took into account all 23 mitochondrial genomes from heteroscleromorphan sponges used for mapping (Table 1). In addition, the mitogenomes of the haplosclerids Amphimedon compressa (NC_010201), Callyspongia plicifera (NC_010206), and Xestospongia muta (NC_010211) were used as outgroups, according to current phylogenetic hypotheses (e.g., Lavrov, Wang & Kelly, 2008; Thacker et al., 2013). All protein coding genes were used as amino acids for the analyses, while the rRNA genes (12S and 16S) were kept as nucleotides. Sequences were aligned individually for each gene using a global alignment algorithm implemented in the software MAFFT v7.505, with default parameters, and then compiled into a single matrix. The software Gblock v0.91b (Castresana, 2000) was used to eliminate ambiguously aligned positions. Phylogenetic reconstructions were performed using the Maximum Likelihood (ML) method implemented in RAxML-NG v1.1.0 (Kozlov et al., 2019). The best-scoring tree was selected among 20 independent runs from random and parsimony starting trees, using, for the protein partition, the mtZOA+F0+G10 model, and for the DNA partition, the GTR+F0+G10 model. Branch lengths were estimated independently for each partition (unlinked). Branch support estimated through 1,000 replicates of bootstrap (or until convergence was achieved) using the Transfer Bootstrap Expectation support metric (Lemoine et al., 2018). Since we were able to recover the complete nuclear rDNA gene cluster of L. hypogea (see below) from the de novo metagenome assembly, we conducted a phylogenetic analysis using 18S sequences from representative species as above since 28S sequences are not available for most of them.

Ethical standards and sample collection permit

Despite the overall deficit of ethical guidelines governing the use of invertebrates in science (see Drinkwater, Robinson & Hart, 2019), all care was taken to minimize animal suffering during the procedures. Sample collection permit to Jean Vacelet, Arrêtê N°107, Prêfecture des Bouches–du-Rhône, Marseille, France.

Genome structure and composition

The complete mitochondrial genome of Lycopodina hypogea was recovered split into three circular molecules: mitochondrial chromosome 1 (chrM1, GenBank accession PP657140), 31,099 bp long; and putative mitochondrial chromosomes 2 (putative chrM2, PP657141), 291 bp, and 3 (putative chrM3, PP657142), 140 bp (Fig. 1; Table 2). As mentioned before, the mtDNA sequence reported by Aquatic Symbiosis Genomics Project (OZ017794) is identical to chrM1 reported here, which includes 13 proteins, 19 tRNAs, and two rRNAs (rnl/16S and rns/12S) genes. As usual for sponges, and in contrast to bilaterians, the mitogenome of L. hypogea lacks an organized control region and includes the atp9 gene. Regarding the transfer RNA genes in the poecilosclerid mitogenomes, L. hypogea chrM1 lacks the trnI(gau), trnL(uaa), trnM(cau), trnR(ucu), and trnT(ugu) genes, found in most poecilosclerids; and the trnY(aua) gene, found in the mitogenome of Negombata magnifica. It also lacks the protein coding gene nd4l. The mitochondrial gene order in chrM1 of L. hypogea differs vastly from other poecilosclerids, including the distancing of the rnl and rns genes. The base composition of L. hypogea mitochondrial genome was: $G=23.6\mathrm{\%}$, $A=25.7\mathrm{\%}$, $T=36.1\mathrm{\%}$, and $C=14.6\mathrm{\%}$ (Table 3). The characteristic intergenic regions in the mitochondrial genome of sponges are highly accentuated in L. hypogea, almost all genes present some sort of IGR between them, ranging from 1–1,791 bp. The L. hypogea mitogenome’s AT-skew was negative (–0.1677) while its GC-skew was positive (0.2362), comprising a higher abundance of Ts and Gs than As and Cs. The chrM1, as assembled from mapped reads, differs from a 31,240 bp long contig recovered from the metagenomic data (see below) solely by a 141 bp long insert that matches exactly an adjacent region, which is likely an assembly artifact, since no read was found directly on the FASTQ files (using the ‘grep’ command) comprising this apparent duplication event. The protein gene nd4l and the tRNA gene trnI(gau), both missing on chrM1, were recovered into distinct circular contigs. These putative mitochondrial chromosomes, putative chrM2 and putative chrM3, differ vastly from other contigs recovered from the de novo metagenomic assembly, comprising short, circular structures. Putative chrM2 comprises solely the nd4l gene, without any IGR between the stop and back to the start codon of the gene, while putative chrM3 comprises trnI(gau) and a short (67 bp) non-coding region. We recovered raw reads spanning several sequences the opposite side of these circular chromosomes in which both start and stop codon of nd4l, in putative chrM2, and reads presenting the full putative chrM3 sequence were found (Supplemental Material S1). The mitochondrial genome of L. hypogea comes from a specimen (voucher MNRJ 22102) belonging to a population cultured in aquariums from 2015 to 2018, originally sampled from a Mediterranean submarine cave (3PP Cave, La Ciotat, France). The ASGP’s sample (BioSample SAMEA9463981), whose mtDNA (OZ017794) is identical to chrM1 reported here, was collected in the same locality, but directly preserved from the field for DNA purification. Thus, short-term culturing of L. hypogea poses no effect over mtDNA organization and overall structure, despite evidences for relaxation of selective constraints, and higher substitution rates in the mitogenome of animals (e.g., Björnerfeldt, Webster & Vilà, 2006; Moray, Lanfear & Bromham, 2014).

In the possibility the unusual features in L. hypogea mtDNA, in particular the putative small circular chromosomes, could be artifactual, we carefully examined alignments of mapped reads, and redo the analyses several times, with different software and parameters in order to corroborate our results, to no change whatsoever (TS de Paula, 2023–2024, personal observation). We also performed PCR amplifications from the genomic DNA of the sample using specific primers designed here (Supplemental Material S1). These were conducted in order to validate arrangements and organization of the mtDNA of L. hypogea as presented, most of which were corroborated to some level given the limitations of the method. These amplifications were conducted prior the L. hypogea genome assembly reported by Aquatic Symbiosis Genomics Project (GCA_963969325) became available, which does not report small extra chromosomes for the mitochondrial genome. Through BLAST searches (Supplemental Material S2), we were able to find partial nd4l sequences in nuclear chromosome 1 (74 bp, 92% identity). In nuclear chromosome 2, we found a 324 bp long, single copy ORF (OZ017780:13,3345,923..13,346,246) encoding the first 235 bp of nd4l gene in putative chrM2 of L. hypogea, but highly divergent from the 236^th position till its end; upstream that ORF, we also found a 82 bp copy of the 3’ end of nd4l gene (Supplemental Material S2). This nuclear copy is likely a pseudogene, and seemingly presents a recombination within it. Given the depth of coverage of putative chrM2 in L. hypogea metagenome (396.7 $\times$; see below), we believe nd4l was not loss to the nucleus. Reads from ASGP’s Illumina run (ERR12668765), mapped against putative chrM2 sequence, also assembled into a circular small contig (average depth = 333.1 $\times$; data not shown). This contig differs from the nuclear copy reported above, and only slightly differs from putative chrM2 by presenting a 29 bp duplication from an adjacent region that leads to a frameshift in ORF comprising nd4l gene, which is likely artifactual since such duplicated sequence is not found among raw reads. We were also able to find several hits for putative chrM3 in the nuclear genome of L. hypogea, all comprising partial sequences spread across several chromosomes (Supplemental Material S2). Thus, these nuclear mitochondrial DNA segments (NUMTs), although presenting similarities to putative chrM2 and putative chrM3, do not appear to be functional, and likely comprise pseudogenes.

Phylogenetic analyses

All 13 protein coding genes and both rRNA genes from chrM1 were included in the final alignment. The nd4l gene was disregarded in these analyses since it was found later and on a different, putative chromosome. The percentage of positions remaining after removing ambiguously aligned regions differs greatly among genes. For the proteins, ATP9 and ATP6 contributed the most, relative to their initial positions (100%), while only 56.9% and 37.8% of the positions were kept for the ND6 and ATP8 genes. For the rRNA genes, a high proportion of ambiguously aligned positions were found, which in addition to long indels, resulted in only 49.5% and 47.9% of initial positions surviving. After removing ambiguously aligned regions, the final alignment had 6,877 positions split between the protein and DNA partitions. The protein partition comprised 3,830 amino acids of which 1.25% were gaps and 45.07% were invariant sites, while the DNA partition had 3,047 nucleotides, 5.85% of those were gaps and 44.27% were invariant sites. The phylogenetic reconstructions took into account 131 free parameters (model + branch lengths) and bootstrap support was achieved after convergence at 250 replicates. The best-scoring tree (Fig. 2) shows the L. hypogea mitogenome sequence closely related to other poecilosclerid sponges, despite its disparate features. The reconstruction also recovers a highly divergent mitogenome sequence for L. hypogea, as depicted by its long branch, particularly in contrast to its ordinal sister species and its own divergence on the 18S tree. No inner node support was achieved for relationships within Poecilosclerida. The phylogenetic reconstruction using 18S data was similar to the mtDNA tree, but for the recovering of Clionaida closer to Poecilosclerida. All other heteroscleromorphan relationships are in agreement to previous surveys (Lavrov, Wang & Kelly, 2008; Hajdu et al., 2013; Redmond et al., 2013; Thacker et al., 2013).

Protein coding genes and codon usage

The thirteen protein coding genes in the chrM1 of Lycopodina hypogea mitogenome covered 38.4% of its total length, encoding a total of 3,969 amino acids (Table 3). The protein genes included six NADH dehydrogenases (ND1–ND6), three cytochrome c oxidases (COX1–COX3), three ATP synthases (ATP6, ATP8 and ATP9) and one cytochrome b (CYTB), ranging in size from 216 bp (ATP8) to 1,839 bp (ND5). As usual for sponges, all thirteen genes lie on the heavy (H) strand. Ten out of the thirteen proteins used the typical initiation codon AUG, while ATP9, COX1 and ND6 used GUG. All protein genes were terminated with the stop codon UAA, but ATP8 and COX1, terminating with UAG. All 13 protein coding genes followed the AT-skew (negative) and GC-skew (positive) pattern from the overall mitogenome. The codon degeneracy pattern for the protein genes encoded in L. hypogea mitogenome (Fig. 3A) shows the amino acids arginine (R), serine (S), and leucine (L) encoded by six synonymous codons, and the remaining amino acids encoded by four (A, G, P, T, and V) or two (C, D, E, F, H, K, N, Q, W, Y) codons, but methionine (M) and isoleucine (I), with one and three codons, respectively. Overall, both two- and four-fold degenerate codons are biased toward the use of codons with U, and against the use of C, in the third position. Amino acid usage across all protein genes revealed leucine (L; 13.5%), glycine (G; 10.3%), and isoleucine (I; 9.0%) as the most used ones, and glutamine (Q; 2.2%), histidine (H; 1.6%), and cysteine (C; 1.4%) as the least used amino acids (Fig. 3B). The loss of nd4l from chrM1 in L. hypogea would be unprecedented in the mitogenome of sponges if not for its occurrence the putative chrM2, which also lies in the heavy strand and uses AUG and UAG as start and terminating codons. The nucleotide composition in the 3^rd codon position on nd4l gene (nd4l: $G_{3rd}=44.3\mathrm{\%}$, $A_{3rd}=20.6\mathrm{\%}$, $T_{3rd}=19.6\mathrm{\%}$, and $C_{3rd}=15.5\mathrm{\%}$) differs greatly from overall protein genes in chrM1 (overal: $G_{3rd}=21.3\pm 4.2\mathrm{\%}$, $A_{3rd}=28.4\pm 2.6\mathrm{\%}$, $T_{3rd}=38.0\pm 6.3\mathrm{\%}$, and $C_{3rd}=12.2\pm 3.7\mathrm{\%}$). The loss of some protein coding genes across non-bilaterian mitogenomes is relatively common (Lavrov & Pett, 2016). Usually, the loss of mitochondrially encoded genes are followed by their transfer to the nucleus (i.e., gene loss to the nucleus), as it seems to be the case of atp9 in Amphimedon queenslandica (Erpenbeck et al., 2007). Conversely, atp8 was thought to have been lost in the mitogenome of glass sponges by many authors (e.g., Haen et al., 2007; Rosengarten et al., 2008; Haen, Pett & Lavrov, 2014) until it was later found in Oopsacas minuta (Jourda et al., 2015). The extreme dissimilarity of atp8 sequences among glass sponges (with an average kimura-2-parameter distance of 0.686 $\pm$ 0.065) was likely one reason it had not been easily found by previous authors. However, the loss of nd4l from the main mitochondrial chromosome, chrM1, to a putative small extrachromosomal circular DNA, chrM2, is until now unknown for any organism, which remains to be properly proved by experimental means.

Transfer RNAs and ribosomal RNAs

There were only 19 tRNAs in the main mitochondrial chromosome, chrM1, of L. hypogea, in contrast to the usual 24–25 found in other poecilosclerids. All 19 tRNAs are encoded on the H-strand (Fig. 1) and present a typical clover structure (Supplemental Material S3). The tRNA genes comprise 4.4% of the length of L. hypogea mitochondrial genome, ranging 71–73 bp in size. The base composition of tRNAs genes was $G=24.2\mathrm{\%}$, $A=28.1\mathrm{\%}$, $T=32.1\mathrm{\%}$, and $C=15.5\mathrm{\%}$ (Table 3), and their AT and GC skew values of tRNA genes matched the mitogenome overall, but a slight prevalence of Ts over As. Among the 19 tRNAs, the amino acid serine (S) presented two genes with distinct anti-codon sites (GCU and UGA), while the amino acids tyrosine (Y), methionine (M), leucine (L), and threonine (T) presented no mitochondrial tRNA gene. The only amino acids which tRNAs were prevalently or equivalently used in contrast to other anti-codons were glutamate (Q), lysine (K), glycine (Q), arginine (R), and tryptophan (W). For all other amino acids, the anti-codon sites from mitochondrial tRNAs presented a lower prevalence across codon usage. The $tRN{A}_{GAU}^{Ile}$ was found encoded elsewhere, in the putative chrM3 by itself. Several NUMT sequences from trnI(gua) were found in the nuclear genome of L. hypogea, in different degrees of mismatch and coverage (see Supplemental Material S2). We were unable to find significant hits to mitochondrial trnM(cau) and trnT(ugu) genes from other poecilosclerids in the nuclear genome of L. hypogea, which may be i) too divergent to be found by BLAST searches or ii) putatively lost from the sponge genome. However, we did found hits to trnL(uaa), trnR(ucu), and trnY(aua) genes, which may be lost to the nucleus, with promising hits on chromosomes 4 and 14, 13, and 6 and 7, respectively. The genes for the large subunit (rnl; 16S rRNA) and small subunit (rns; 12S rRNA) of the ribosomal RNA are both located on the H-strand but appear on opposing sides of L. hypogea mitogenome. The two rRNA genes are, respectively, 2,526 and 1,197 bp long, and present a similar base composition (Table 3). Contrary to the overall mitogenome, the AT-skew value of both genes were positive (0.0303 and 0.0315), indicating slightly more adenine in both, while the GC-skew values remained about the same, but about 40% higher in the rns gene (0.3655), indicating a higher concentration of guanine nucleotides in that rRNA (Table 3).

Gene rearrangements and IGRs

Compared with the previously known and relatively homogeneous gene order in the mitochondrial genome of Poecilosclerida, the Lycopodina hypogea mitochondrial chromosome 1 (chrM1) underwent profound and obvious rearrangements (Fig. 4). Given the sequencing depth of L. hypogea mitogenome across all sites (Fig. 1B), it is unlikely that its gene order, as recovered through de novo assembly, is artifactual. Seemingly, the position of all genes had changed, and the L. hypogea mitogenome is packed with intergenic regions (IGRs). The IGRs comprise 45.2% of the length of L. hypogea mitochondrial genome, and fall in five categories given their size: I, very short, ranging 1–20 bp ( $N=11$); II, short, ranging 66–109 bp ( $N=5$); III, intermediate, ranging 198–426 bp ( $N=7$); IV, long, ranging 592–1,108 bp ( $N=6$); and V, very long, ranging 1,496–1,791 bp ( $N=4$). The base composition of all IGRs was, on average, $G=24.5\pm 2.7\mathrm{\%}$, $A=25.2\pm 4.3\mathrm{\%}$, $T=35.3\pm 3.7\mathrm{\%}$, and $C=15.1\pm 3.1\mathrm{\%}$ (Table 3), and their AT- and GC-skew values matched the overall mitogenome. BLAST searches on most IGRs returned non significant results, but four: IGR 2527–3123, between the rnl and trnI(cau) genes; IGR 7741–8332, between trnS(uga) and trnC(gca); IGR 14911–15019, between cytb and trnA(ugc); and again in IGR 30083–30843, between nd1 and trnY(gua). All these positive BLAST results comprise hits with 42–74 bp and 79–92% similarity with terrestrial animals, three artropods, Lagria hirta (Coleoptera), Chrysoteuchia culmella (Lepidoptera), and Timema monikensis (Neoptera), and a mammal, Bos mutus (Bovidae), but with no hits to known functional regions.

The large size of mitochondrial chromosome 1 of Lycopodina hypogea is due solely to the prevalence and size of its IGRs. Although no MITE sequence was detected, all long and very long IGRs presented at least one imperfect terminal direct repeat (TDR), up to 67, and most presented 5–38 imperfect terminal inverted repeats (TIR) longer than 10 bp. The length of IGR was proved to be a significant factor for the observed number of TDRs ( $R^2=0.9032$; $F=121.2$; $p<0.0001$) and TIRs ( $R^2=0.8326$; $F=64.6$; $p<0.0001$). Large mitogenomes (over 30 Kbp) with extensive IGRs (comprising 30–65% of their sequences) are highly unusual in metazoans, with cases found among freshwater sponges (Lavrov et al., 2012), placozoans (Signorovitch, Buss & Dellaporta, 2007), arcid bivalves (Kong et al., 2020), and amphipods (Romanova et al., 2021). For freshwater sponges, repetitive hairpin-forming elements were abundantly found in their mitogenomic IGRs. For the amphipod mitogenome, conversely, vestiges of rRNA could be found within their IGRs. Both examples are associated with two processes that are argued to be involved in the extension of IGRs in animal mitogenomes: i) duplication of genetic regions and subsequent loss of redundant copies (Moritz, Dowling & Brown, 1987; San Mauro et al., 2006; Lavrov, Boore & Brown, 2002; Schirtzinger et al., 2012; Shi, Miao & Kong, 2014); and ii) proliferation of repetitive elements within the mitogenome (Erpenbeck et al., 2009; Lavrov, 2010; Lavrov et al., 2012). It remains to be proved which (if any) of these processes are responsible for the extensive IGRS in chrM1 of L. hypogea.

Metagenomic analyses

Although a complete metagenome analysis of L. hypogea is beyond the scope of this study, the analysis of deeply sequenced reads was able to reveal that some common mitochondrial genes, the protein gene nd4l and tRNA gene trnI(gua), missing from the main mitochondrial genome, chrM1, were recovered each in a separate contiguous, circular sequence (Fig. 5). After careful examination of these contigs, which were artifactually assembled due to looped sequences, we concluded they comprise putative mitochondrial chromosomes 2 and 3 (see above). The analysis of these contigs also revealed that most of these sequences comprise genes from mobile elements, either transponsable elements, or genes from prophages or plasmids (Material S4). Analysis of ORFs extracted from sequences belonging to mitochondrial IGRs did not return significant results, meaning they present no homology to proteins curated in the database. This contrast with the deeply sequenced contigs from sources other than the mitochondrial or rDNA of L. hypogea, in which almost all comprised at least one positively identified ORF. In addition, no read mapped onto chrM1 was found among reads mapped onto these other contigs. This suggests mitochondrial IGRs of L. hypogea have no direct homology to these prevalent mobile elements from distinct sources in particular. Transponsable elements are not only the most abundant and ubiquitous genetic sequences in nature (Aziz, Breitbart & Edwards, 2010), they are also thought to provide elevated rates of change and adaptation in bacterial populations (e.g., Brazelton & Baross, 2009; Vigil-Stenman et al., 2017).

Small extrachromosomal circular DNAs

We are relatively secured that putative mitochondrial chromosomes 2 and 3 do not comprise single, linear molecules. And although we did not conclusively exclude the possibility they comprise tandem repeats within longer molecules, the more suitable explanation for their presence in L. hypogea metagenome is in the form of small and circular molecules (Material S1). While the nature of these molecules is de facto unclear, pending further analyses such as Southern blotting and DNA cloning, indirect evidences lean towards one hypothesis over the other. Firstly, given how deeply sequenced putative chrM2 and chrM3 were among reads (Fig. 1B), the diversity around the ends of their sequence, which would result from sequences flanking the tandem repeat, was null, which is highly improbable, unless they cover several dozens of repeating units. In addition, the genomic DNA of L. hypogea shows segments of the expected size range of putative chrM2 and chrM3. Assuming these gDNA segments comprise the sequences from the assembled contigs (since it would be uncanny for them to share the expected sizes but not their sources), to comprise linear fragments from longer molecules bearing tandem repeats, they would have to have been naturally cut in the exact same or on a single place, while other gDNA molecules remained intact, with high molecular weight. Thus, it seems more likely they are organized into small circular molecules. But in the absence of further evidences, we opt to regard them as putative mitochondrial chromosomes.

The organization of putative chrM2 and chrM3 into small circular molecules is congruent with the growing body of knowledge about extrachromosomal circular DNA (eccDNA) (Paulsen et al., 2018; Ain et al., 2020; Liao et al., 2020; Zhao et al., 2022). These small (<400 bp) eccDNAs, called microDNAs, or small polydisperse circular DNAs (spcDNA), as initially denoted (Cohen, Regev & Lavi, 1997), are associated, in humans, to genomic instability in malignant tumors, cellular degeneration, and senescence processes (Ain et al., 2020; Liao et al., 2020), but are also widely dispersed in healthy tissues (Dillon et al., 2015; Mϕller et al., 2018). Many mechanisms are seemingly involved in the biogenesis of microDNAs. In the nucleus, microDNA are formed mainly by DNA repair pathways, especially after double-strand breaks and during replication. However, since eccDNA sequences may comprise a substantial proportion of transposable elements, such as LINE-1 retrotransposons, and long terminal repeats (LTRs) (Dillon et al., 2015; Mϕller et al., 2016, 2018, 2020), some authors suggest regions associated to DNA:RNA hybrids (R-loops) are structurally prone to eccDNA biogenesis (Dillon et al., 2015). R-loops are also major contributors to genome stability (Brickner, Garzon & Cimprich, 2022), specially for mitochondrial DNA replication, since R-loops are important intermediates for the D-loop formation in mtDNA control region (Holt, 2019; Brickner, Garzon & Cimprich, 2022). Thus, it is a striking coincidence to find both high levels of transposable elements and putatively microDNAs of mitochondrial origin in the metagenome of L. hypogea.

A strong opposing argument can be drawn from the fact that both putative chrM2 and chrM3 lack canonical signaling sequences. Although they could unknowingly to us exist within their sequences, in the negative case, to be viable mitochondrial extrachromosomes, these microDNAs must be replicated and transcribed independent of such signalling mechanism. It was already shown that microDNAs are transcribed in vivo even without a canonical promoter sequence (Paulsen et al., 2019), thus they theoretically could be expressed within L. hypogea mitochondria. The lack of an organized control region in the mitochondrial DNA of sponges suggests its replication mechanism may differ from their bilaterian relatives (Lavrov & Pett, 2016). Given the primitive features of sponge mtDNAs, it is reasonable to assume they may also share similarities to the replication process observed in other opistokhonts, such as yeast. Although yeasts mtDNA may present rep/ori elements containing a consensus transcriptional promoter that may control replication through RNA polymerase primed initiation, other mechanisms must exist since several yeast species and strains lack rep/ori segments (Ling, Hori & Shibata, 2007; Chen & Clark-Walker, 2018). Thus, there is strong evidence suggesting yeast mtDNA could be replicated via a recombination-mediated initiation mechanism, in which mtDNA undergo homologous (catalysed or uncatalyzed) recombination that primes a rolling-circle replication, a mechanism that resembles replication and packaging of phage DNA. If sponges present similar mechanisms, mtDNA replication may bypass the requirement of signalling sequences in either normal mitochondrial chromosomes or putative microDNAs.

We also found a third deeply sequenced, circular contig within the metagenome of L. hypogea, but at this time with no significant orthologs hit. Given that we poorly understand the role (if any) of microDNAs in the evolution of sponge mitogenomes, this contig could comprise another mitochondrial chromosome, (see n.s. contig at Fig. 5A). After further analysis, we determined this microDNA presents 133 bp, containing small ORFs with no orthologs found in the searched databases. Thus, we opt to disregard it as an additional mitochondrial chromosome of L. hypogea until its origin/function is revealed.

Nuclear ribosomal DNA cluster

As part of the efforts to find the missing mitochondrial genes among metagenomic contigs, we were able to recover the full nuclear ribosomal DNA (nrDNA) cluster of L. hypogea, which is 9,046 bp long and includes, in order: the intergenic spacer (IGS); 5′ external transcribed spacer (5′ ETS); 18S ribosomal RNA gene (18S); internal transcribed spacer 1 (ITS1); 5.8S ribosomal RNA gene (5.8S); internal transcribed spacer 2 (ITS1); 28S ribosomal RNA gene (28S), internal transcribed spacer 3 (ITS3), 5S ribosomal RNA gene (5S), and 3′ external transcribed spacer (3′ ETS) (Genbank accession PP658212; Fig. 6). This is noteworthy since, to the best of our knowledge, this is the first time a full rDNA cluster is described for sponges that shows the 5S rRNA gene clustered in the same cassette with 18S, 5.8S and 28S rRNA genes, including a third internal transcribed spacer (ITS3). The boundaries between 3′ ETS and IGS, and IGS and 5′ ETS, are unclear, but reads mapped onto the IGS and ETS regions evidence they are contiguous, what would be expected for a multi-copy repeat of clusters. Here they seem delimited by major hairpin formations on their secondary structures (data not shown), what must be reviewed after analysis of pre-rRNA sequences. However, we choose to tentatively distinguish transcribed from the non-transcribed, intergenic spacers since i) this arrangement could be useful for designing primers to amplify the 5S rRNA gene and the 5′ end of the 18S rRNA gene, and ii) comparisons with rRNA sequences extracted from transcriptomic data might help elucidate the structure of these regions. Evidence for intragenomic variations (IGVs) can be found among reads, in particular in the IGS region, which presents variation in number of repeats (easily noticeable by a sharp drop in coverage of reads mapped onto it). For instance, a poly- $(\mathrm{T}{)}_{11-17}$ repeat at 175–190 bp of the rDNA, and the poly- $(\mathrm{A}{)}_{25-31}$ repeat at 490–515 bp, or the motif $(\mathrm{G}\mathrm{C}\mathrm{A}\mathrm{A}{)}_{18-24}$ at 674–745 bp. Other repeats in the IGS regions, such as the motifs $(\mathrm{C}\mathrm{C}\mathrm{A}\mathrm{A}{)}_7$ at 631–658 bp, and the $(\mathrm{C}\mathrm{T}\mathrm{T}\mathrm{G}{)}_{12}$ at 878–925 bp, or poly-nucleotide and terminal direct repeats are also present, but seemingly without variation among reads (Fig. 6). These repeating sequences are features commonly found in the IGS region of almost all species, from plants (e.g. Hu et al., 2019), to yeasts (e.g., James et al., 2009), and animals (e.g., Dyomin et al., 2019; Hori, Shimamoto & Kobayashi, 2021). We were able to find the partial 28S sequences spread in L. hypogea nuclear chromosome 9, with a cluster of three repeats ranging 8,660–9,148 bp each (OZ017785:2,841,318–2,868,141, minus strand). With exception of variation within the IGS region, those repeats are in agreement with the rDNA cluster presented here.