<title>Introduction</title> <p>Two major evolutionary events occurred early in animal history and shaped the majority of animals, as we know them today: the origin of multicellularity and the origin of bilateral symmetry. The phylogenetic boundaries of these events are well defined among extant taxa and correspond to the traditional groups Metazoa (multicellular animals) and Bilateria (all animal phyla except Porifera, Placozoa, Cnidaria, and Ctenophora). Multiple genomic changes must have occurred in association with these morphological transitions, and current genome sequencing projects give us the first glimpses into these changes ##REF##18273011##[1]##, ##REF##17615350##[2]##.</p> <p>Surprisingly, the transitions to multicellular and bilaterally symmetrical animals also correlate with multiple changes in mitochondrial genome architecture ##REF##21669754##[3]##, although the main function of mitochondria themselves remained unchanged. In particular, the origin of animal multicellularity is associated with the loss of all ribosomal protein genes from mtDNA, the disappearance of most introns, and a large reduction in the amount of non-coding DNA ##REF##21669754##[3]##. The origin of bilaterality correlates with further compaction of mtDNA, multiple changes in the genetic code and the associated losses of some tRNA genes, along with the appearance of several genetic novelties ##REF##1282405##[4]##. Obviously, the picture presented above is an extrapolation of our knowledge of extant organisms into the ancient past and as such can be affected by artifacts of ancestral state reconstruction ##REF##21238344##[5]##. It is also based on a relatively limited sampling of mitochondrial genomes, especially from non-bilaterian animals, and additional data from Cnidaria, Ctenophora, Porifera, as well as the closely related lineages of eukaryotes (e.g., Choanozoa) are essential to support, expand, or refute it.</p> <p>Class Demospongiae ##UREF##0##[6]## is the largest (&gt;85% of species) and most morphologically diverse group in the phylum Porifera. It contains sponges of various shapes and sizes that occupy both freshwater and marine environments from shallow to abysmal depths and includes such oddities as carnivorous sponges ##UREF##1##[7]##. Within the extant Demospongiae 14 orders are recognized that encompass 88 families, 500 genera and more than 8000 described species ##UREF##2##[8]##, ##UREF##3##[9]##. Although traditionally three subclasses have been distinguished, two of them do not appear to be monophyletic. Instead, recent molecular studies ##REF##15288059##[10]##, ##UREF##4##[11]## provide strong support for five major clades within the Demospongiae: Homoscleromorpha (G0) (Homosclerophorida), Keratosa (G1) (Dictyoceratida+Dendroceratida), Myxospongiae (G2) (Chondrosida, Halisarcida, and Verongida), Marine Haplosclerida (G3), and all the remaining groups (G4) (##FIG##0##Figure 1##). Our knowledge of mtDNA diversity within the demosponges has been rudimentary, with only five sequences representing 3 of the 5 major groups available ##REF##15703239##[12]##–##REF##17053047##[15]##. Previous studies revealed that demosponge mtDNA resembles that of most other animals in its compact organization, lack of introns, and well-conserved gene order, but at the same time contains extra genes, including <italic>atp9, trnI(cau), trnR(ucu),</italic> encodes bacterial-like ribosomal and transfer RNAs, and uses a minimally derived genetic code in protein synthesis ##REF##15703239##[12]##. Furthermore, additional unusual features found in the mitochondrial genomes of <italic>Oscarella carmela</italic>
##REF##17090697##[14]## and <italic>Amphimedon queenslandica</italic>
##REF##17053047##[15]## suggested that more mitochondrial genomic diversity might exist among the demosponges. Here we describe complete mitochondrial sequences from 17 species of demosponges and analyze them with five previously published mitochondrial genomes from this group that were available at the time this study was conducted. Taken together, our sampling covers all recognized order-level diversity within the Demospongiae and provides the first analysis of general evolutionary trends in mitochondrial genome organization for this group. Such a comprehensive approach to the analysis of demosponge mtDNA is needed because, at least in the fossil record, the evolution of demosponges closely mirrors the evolution of all bilaterian animals with the first demosponge fossils appearing in Precambrian deposits and a major radiation occurring in the Lower Cambrian ##REF##15665105##[16]##, ##UREF##5##[17]##.</p>
<title>Methods</title> <title>Genome sequencing and phylogenetic analysis</title> <p>Taxon sampling, DNA extraction, PCR amplification, and sequencing were described in our previous article ##UREF##4##[11]##. Phylogenetic analysis of demosponge relationships was conducted with the PhyloBayes program ##REF##15014145##[38]## as described previously ##UREF##4##[11]##, except that mitochondrial sequences from several taxa (mostly cnidarians) have been added: <italic>Agaricia humilis</italic> NC_008160, <italic>Anacropora matthai</italic> NC_006898, <italic>Aphrocallistes vastus</italic> EU000309, <italic>Branchiostoma floridae</italic> NC_000834, <italic>Capsaspora owczarzaki</italic>, <italic>Colpophyllia natans</italic> NC_008162, <italic>Discosoma</italic> sp. CASIZ 168915 NC_008071, <italic>Hydra oligactis</italic> EU237491, <italic>Montipora cactus</italic> NC_006902, <italic>Mussa angulosa</italic> NC_008163, Placozoan sp. BZ2423 NC_008834, Placozoan sp. BZ49 NC_008833, <italic>Pocillopora damicornis</italic> NC_009797, <italic>Porites porites</italic> NC_008166, <italic>Pseudopterogorgia bipinnata</italic> NC_008157, <italic>Rhodactis</italic> sp. CASIZ 171755 NC_008158, <italic>Seriatopora caliendrum</italic> NC_010245, <italic>Siderastrea radians</italic> NC_008167.</p> <title>Annotation and analysis of coding sequences</title> <p>We used flip (<ext-link ext-link-type="uri" xlink:href="http://megasun.bch.umontreal.ca/ogmp/ogmpid.html">http://megasun.bch.umontreal.ca/ogmp/ogmpid.html</ext-link>) to predict ORFs in assembled sequences; similarity searches in local databases and in GenBank using FASTA ##REF##8004177##[39]## and NCBI BLAST network service ##REF##18073190##[40]##, respectively, to identify them. Protein-coding genes were aligned with their homologues from other species and their 5′ and 3′ ends inspected for alternative start and stop codons. Inferred amino acid sequences of encoded proteins were aligned with ProbCons ##REF##15687296##[41]## using default parameters. Genetic distances between demosponges and four species of cnidarians (<italic>Briareum asbestinum, Metridium senile, Montastraea annularis</italic> and <italic>Ricordea florida</italic>) were calculated with the TREE-PUZZLE program ##REF##11934758##[42]##, using the mtREV matrix, estimated frequencies of amino acids and 8 gamma rate categories. Effective numbers of codons ##REF##2110097##[43]## were calculated with the chips program within the EMBOSS package ##REF##10827456##[44]##.</p> <title>Annotation and analysis of RNA genes</title> <p>Genes for small and large subunit ribosomal RNAs (<italic>rns</italic> and <italic>rnl</italic>, respectively) were identified based on their similarity to homologous genes in other species, and their 5′ and 3′ ends were predicted based on sequence and secondary structure conservation. The secondary structures of selected rRNA genes were manually folded by analogy to published rRNA structures, and drawn with the RnaViz 2 program ##REF##12538259##[45]##.</p> <p>Transfer RNA genes were identified by the tRNAscan-SE program ##REF##9023104##[46]## and aligned manually in MacGDE 2.3 ##REF##7704666##[47]## using their secondary structure as a guide. This alignment was used to calculate sequence conservation at each position and average pairwise identity values for individual tRNAs. For the latter calculation we excluded all tRNAs from <italic>Plakortis angulospiculatus, Amphimedon queenslandica</italic> and all species in G1, which encode incomplete sets of tRNAs in their mtDNA.</p> <title>Intronic sequences</title> <p>We used intron prediction programs RNAweasel ##REF##17280737##[48]## and Rfam ##REF##12520045##[49]## to search for introns in coding sequences. The exact positions of two introns found in <italic>cox1</italic> of <italic>Plakortis angulospiculatus</italic> were adjusted based on <italic>cox1</italic> alignments with homologous sequences from other demosponges. The inferred amino acid sequence of the large ORF found in one of the <italic>P. angulospiculatus</italic> introns was aligned with the sequences of LAGLIDADG ORFs analyzed by Rot <italic>et al.</italic>
##REF##16972986##[23]## and Fukami <italic>et al.</italic>
##REF##17437148##[24]## and used for a phylogenetic analysis. We selected the best model for these ORFs with the ProtTest program ##REF##15647292##[50]## and performed a maximum likelihood search and bootstrap analysis in TREEFINDER ##REF##15222900##[51]##, using the WAG model of sequence evolution, estimated amino acid frequencies and 4 gamma categories.</p> <title>Intergenic regions and repeated sequences</title> <p>Intergenic regions were extracted from each genome with the PEPPER program (<ext-link ext-link-type="uri" xlink:href="http://megasun.bch.umontreal.ca/ogmp/ogmpid.html">http://megasun.bch.umontreal.ca/ogmp/ogmpid.html</ext-link>) and searched for similarity using FASTA. In addition, we searched for interspersed identical repeats in individual genomes using FINDREP (<ext-link ext-link-type="uri" xlink:href="http://megasun.bch.umontreal.ca/ogmp/ogmpid.html">http://megasun.bch.umontreal.ca/ogmp/ogmpid.html</ext-link>) with minimum repeat subsequence lengths of 20 bp and 100 bp respectively.</p>
<title>Results</title> <title>Genome organization and nucleotide composition</title> <p>All twenty-two analyzed mtDNAs of demosponges were circular-mapping molecules, each containing a conserved set of thirteen protein-coding and two rRNA genes identical to that found in the mtDNA of most bilaterian animals ##REF##10101183##[18]##. In addition, <italic>atp9</italic>, a gene for subunit 9 of ATP synthase was identified in mtDNA of all demosponges except <italic>Amphimedon queenslandica</italic>
##REF##17053047##[15]##, and <italic>tatC,</italic> a gene for twin arginine translocase subunit C, was found in <italic>Oscarella carmela</italic>
##REF##17090697##[14]##. The number of tRNA genes showed more variation. Although 24 or 25 tRNA genes were present in most analyzed demosponge mitochondrial genomes, as few as 2 and as many as 27 tRNA genes were found in mitochondrial genomes of some demosponge species (##FIG##1##Figure 2##, see below). In addition, a sequence with a potential to form a tRNA-like structure, named <italic>trnX</italic>, was located downstream of <italic>cox1</italic> in <italic>Xestospongia muta</italic> and <italic>Ephydatia muelleri</italic> mtDNA. Inferred tRNA(X) had a well-conserved primary (65.3% nucleotide identity) and secondary structure, except for the putative anticodon arm, which differed both in length and in sequence between the two species.</p> <p>The sampled demosponge mitochondrial genomes displayed moderate size variation (16–26 kb; mean = 19.7 kb), most of which could be attributed to the expansions of non-coding regions usually caused by the presence of repetitive elements (##FIG##1##Figure 2##). We detected no obvious phylogenetic pattern associated with this variation, and no similarity in the sequence of repetitive elements among different species. Most demosponge mitochondrial genomes were larger than their counterparts in bilaterian animals. However, even the largest demosponge mitochondrial genomes were dwarfed in comparison to those in the choanoflagellate <italic>Monosiga brevicollis</italic> and the placozoan <italic>Trichoplax adhaerens</italic>, which have a much higher percentage of non-coding DNA and, in the case of <italic>M. brevicollis</italic>, an expanded gene set (##FIG##1##Figure 2##).</p> <p>All analyzed mitochondrial genomes were relatively uniform in the overall nucleotide composition (A+T content between 56–72%) and, on average, displayed negative AT- and positive GC-skews of the coding strand (##FIG##2##Figure 3##). The sense strand of protein and tRNA genes had a negative AT-skew in all species, that of rRNA genes had a positive AT-skew, while non-coding regions and 3<sup>rd</sup> codon positions showed a large variation in AT-skew both among and within major demosponge groups (##FIG##2##Figure 3B##). All types of sequences in demosponge mtDNAs showed positive GC-skews except for the tRNA genes in <italic>Igernella notabilis</italic> and the non-coding regions in <italic>Ephydatia muelleri</italic> and <italic>Aplysina fistularis</italic>. The genomic values for AT- and GC-skews correlated more strongly with those for protein genes (R<sup>2</sup> = 0.89 and 0.95, respectively) and rRNA genes (R<sup>2</sup> = 0.61 and 0.93) than those for tRNA genes (R<sup>2</sup> = 0.06 and 0.57) and non-coding regions (R<sup>2</sup> = 0.13 and 0.34), while genomic A+T content correlated most strongly with that of rRNA genes (R<sup>2</sup> = 0.89) comparing to non-coding regions (R<sup>2</sup> = 0.78), tRNA genes (R<sup>2</sup> = 0.65), and protein genes (R<sup>2</sup> = 0.44). Interestingly, non-coding regions and 3<sup>rd</sup> codons (that are usually assumed to experience similar mutational pressure) showed little correlation in all three types of measurements (R<sup>2</sup> values are 0.05, 0.3 and 0.58 for A+T content, AT- and GC-skews, respectively).</p> <p>All but four pairs of sampled mitochondrial genomes had unique gene orders, with the number of shared gene boundaries between individual genomes ranging from 1 to 41. The extent of gene order variation and the type of gene rearrangements differed among major groups of demosponges (##FIG##3##Figure 4##). Gene arrangements of protein and rRNA genes were generally well conserved within G2, G3 and G4 and the predominant type of change within these groups was tRNA transposition. By contrast, more rearrangements were found within G0 (13% of shared boundaries between two sampled genomes) as well as within G1 (59% of shared boundaries among three genomes). Still, most of the rearrangements were transpositions and only two inversions were found in the whole dataset (in <italic>Oscarella carmela</italic> and <italic>Aplysina fistularis</italic>).</p> <title>Protein coding genes</title> <p>The protein coding genes identified in all 22 demosponge mtDNAs showed 0.33–11.81% variation in size and 31.9–87.3% average pairwise identity calculated based on inferred amino acid sequences (##SUPPL##3##Table S1##). <italic>Atp8</italic> was the least conserved gene both in terms of size (11.81% variation), pairwise sequence identity among demosponges (31.9% on average, range 8.5–85.7%), and genetic distance to cnidarian homologues (##FIG##4##Figure 5##), followed by <italic>nad6</italic>. By contrast, <italic>atp9</italic>, a gene encoding another subunit of the ATP-synthase complex, was the most conserved, with an average pairwise identity of 87.3% (range 76.9%–100%). Other genes were relatively uniform both in their average pairwise identities across the demosponges and the calculated rates of sequence evolution (##FIG##4##Figure 5##).</p> <p>Codon usage in all analyzed demosponge mitochondrial genomes was consistent with the minimally modified genetic code inferred in our previous study ##REF##15703239##[12]##. All 22 mtDNAs share similar codon usage bias with an effective number of codons equivalent to 41.8±3.5. Synonymous codons ending with A or T were clearly preferred (56–85% for individual species; 73.6% on average), while the codon CGC was not used at all in mitochondrial coding sequences of 12 species. <italic>Tethya actinia</italic> displayed the most biased mitochondrial codon usage with no AAC, CGC, CTC, CTG, and TGC codons present.</p> <p>ATG was the most common initiation codon, followed by GTG, which occurred frequently in <italic>nad6</italic> (15 out of 22 species) and occasionally in other genes (##SUPPL##4##Table S2##). The unusual start codon ATT was inferred for <italic>cox2</italic> in <italic>Hippospongia lachne</italic>, <italic>nad3</italic> in <italic>Cinachyrella kuekenthali</italic> and <italic>nad6</italic> in <italic>Vaceletia</italic> sp. and a TGG start codon was inferred for <italic>nad2</italic> in <italic>Ephydatia muelleri</italic>, <italic>nad6</italic> in <italic>Tethya actinia</italic>, <italic>Axinella corrugata</italic>, <italic>Amphimedon queenslandica</italic> and <italic>tatC</italic> in <italic>Oscarella carmela</italic> (##SUPPL##4##Table S2##). Such initiation codons are common in mitochondrial coding sequences of bilaterian animals ##REF##1282405##[4]##, but are rare, although not unprecedented, in non-bilaterian animals and non-animal outgroups ##REF##7592355##[19]##, ##REF##12591875##[20]##. The stop codons TAA and TAG were inferred for all coding sequences except <italic>nad5</italic> in <italic>Amphimedon compressa</italic>, <italic>Ectyoplasia ferox</italic>, <italic>Ephydatia muelleri</italic>, and <italic>Callyspongia plicifera</italic> as well as <italic>nad4L</italic> in <italic>Cinachyrella kuekenthali.</italic> No standard or abbreviated stop codons were found for the latter genes and the mechanism of their translational termination remains unclear.</p> <p>Among the five major clades within the Demospongiae (G0–G4), a significant acceleration in the rates of evolution was found in G1, especially in the lineage leading to <italic>Vaceletia</italic> sp. and <italic>Hippospongia lachne</italic> (##FIG##4##Figure 5##; RRTree P = 1.00E−07). We tested whether the G1 accelerated rates could have been the result of positive selection as suggested by Bazin <italic>et al.</italic>
##REF##16645093##[21]## but did not find significant support for this hypothesis by either the M1–M2 test in PAML or by the synonymous vs. non-synonymous substitution rate test with the DNASP program ##REF##14668244##[22]##.</p> <title>Introns in cox1</title> <p>Although introns are common in mtDNA of two groups of non-bilaterian animals, Cnidaria and Placozoa, only one mitochondrial intron (in <italic>cox1</italic> of <italic>Tetilla</italic> sp.) has been reported so far in demosponges ##REF##16972986##[23]##. Among the 22 demosponge mitochondrial genomes analyzed for this study, we found two additional group I introns, both of them in <italic>cox1</italic> of <italic>Plakortis angulospiculatus</italic>. These introns were 388 bp and 1118 bp in size (henceforth intron 1 and 2, respectively), and separated by only 9 nucleotides (3 codons) in the gene. Intron 2 in <italic>P. angulospiculatus</italic> was found after position 726 in <italic>cox1</italic>, at the same location as the intron reported for <italic>Tetilla</italic> sp. ##REF##16972986##[23]##. Intron 2 in <italic>P. angulospiculatus</italic> and its counterpart in <italic>Tetilla</italic> sp. share 81.2% nucleotide sequence identity, have a similar secondary structure, and both contain an ORF homologous to LAGLIDADG-type homing endonuclease with identical LAGLIDADG motifs (LAGLIEGDG and LAGFLDADG). By contrast, introns 1 and 2 in <italic>P. angulospiculatus</italic> share only 43.5% sequence identity in the aligned overlap regions and intron 1 does not contain any ORF.</p> <p>Recently, group I introns highly similar to, and in the same position as intron 2 in <italic>P. angulospiculatus</italic> and its homolog in <italic>Tetilla</italic> sp. were reported in <italic>cox1</italic> of 20 scleractinian corals ##REF##17437148##[24]##. Phylogenetic analysis of amino-acid sequences derived from intronic LAGLIDADG ORFs in <italic>P. angulospiculatus</italic>, <italic>Tetilla</italic> sp., scleractinian corals, and several outgroup taxa grouped introns found in <italic>Tetilla</italic> sp. and <italic>P. angulospiculatus</italic> with 72% bootstrap support and placed them as a sister group to Scleractinian corals with 100% bootstrap support (##SUPPL##1##Figure S2##). The results of this analysis are consistent with the vertical evolution of this intron in cnidarians and sponges and suggest that its sporadic presence among sampled taxa is due to independent losses rather than the horizontal intron transfer proposed earlier ##REF##16972986##[23]##. This inference is reinforced by the observations that the genetic distance between LAGLIDADG ORFs in <italic>P. angulospiculatus</italic> and <italic>Tetilla</italic> sp. is similar to that between their host genes and that both ORFs contain a TGA codon at the same position (data not shown). The latter finding makes it highly unlikely that the two introns have been transferred in parallel from the nucleus, because TGA signifies a stop codon in cytoplasmic translation.</p> <title>rRNA genes</title> <p>Genes for the small and large subunit ribosomal RNAs (<italic>rns</italic> and <italic>rnl</italic>) were located in close proximity of each other (separated by 1–3 tRNA genes) in most analyzed genomes, with the most common gene order being <italic>+rns+trnG+trnV+rnl</italic> (##FIG##3##Figure 4##). The two exceptions to this pattern were found in <italic>Igernella notabilis,</italic> where the two genes were separated by <italic>atp9</italic>, and <italic>Oscarella carmela,</italic> where <italic>rnl</italic> and <italic>rns</italic> were separated by multiple genes and had opposite transcriptional orientations. The size of <italic>rns</italic> ranged between 828 (<italic>Hippospongia lachne</italic>) and 1516 bp (<italic>Ephydatia muelleri</italic>), with the average size being 1224 bp. The size of <italic>rnl</italic> varied between 2166 (<italic>Hippospongia lachne</italic>) and 3487 bp (<italic>Axinella corrugata</italic>), with the average size being 2589 bp. The size differences in rRNA genes were due to two factors. First, some helices outside the core region of each rRNA were shortened or lost in several lineages, especially G1 (##FIG##5##Figure 6##). Second, unusual repetitive elements (see below) were inserted in rRNA genes in several distantly related species, in particular <italic>Axinella corrugata</italic>, <italic>Ephydatia muelleri</italic>, <italic>Igernella notabilis</italic>, and <italic>Vaceletia</italic> sp. (##SUPPL##2##Figure S3##).</p> <title>tRNA genes</title> <p>Sampled demosponge mitochondrial genomes contained as few as 2 and as many as 27 tRNA genes. The variation in the number of tRNA genes was due to the loss of all but two mitochondrial tRNA genes (<italic>trnM(cau)</italic> and <italic>trnW(uca)</italic>) in G1, partial losses of tRNA genes in <italic>Agelas schmidti</italic> (at least one gene), <italic>Amphimedon queenslandica</italic> (at least 7 genes), and <italic>Plakortis angulospiculatus</italic> (at least 18 genes), the sporadic presence of <italic>trnM(cau)e</italic> among sampled species, and duplication of <italic>trnT(ugu)</italic> and <italic>trnV(uac)</italic> in <italic>Oscarella carmela</italic> mtDNA. Given that at least 24 species of tRNAs are needed for mitochondrial translation in demosponges ##REF##15703239##[12]##, we expect that the loss of tRNA genes from mtDNA is compensated by the import of required tRNAs from the cytoplasm.</p> <p>In accord with our previous study ##REF##15703239##[12]##, tRNA genes in all studied demosponge mtDNA were well conserved in terms of size, primary sequence and inferred secondary structure. All inferred mt-tRNA structures had well conserved D- and T-loops (7–11 and 7 nucleotides in length, respectively) with a potential to form the standard tertiary interactions G18-U55 and G19-C56. Variable or semi-invariable nucleotide positions, and secondary and tertiary interactions known for prokaryotic and nuclear tRNAs were also well conserved (##FIG##6##Figure 7##). At the same time, an unusual A11-T24 pair in  and an unusual G11-C24 pair in  were present among all sampled demosponges. The first of them is characteristic for demosponges, glass sponges, and placozoans ##REF##15703239##[12]##, ##REF##17090697##[14]##, while the second–for all bilaterian animals ##REF##1452431##[25]##. The R11-Y24 pair is otherwise a distinctive feature of bacterial, archaeal, and organellar initiator  that is strongly counter-selected in elongator tRNAs ##REF##12403461##[26]##.</p> <p>Among individual tRNA genes, <italic>trnW(uca)</italic> had the most conserved primary structure (84.9% pairwise sequence identity on average) while <italic>trnS(uga)</italic> was the least conserved (66.7% identity on average). The inferred gene for elongator tRNA(M) (<italic>trnM(cau)e</italic>) that is present in 11 out of 22 analyzed genomes also displayed high sequence conservation (average pairwise identity 79.6%), an observation that suggests its intermittent occurrence among sampled genomes is due to multiple losses rather than <italic>de novo</italic> evolution through gene duplication and/or recruitment ##REF##15734570##[e.g., 27]##. Interestingly, the other gene for methionine tRNA (<italic>trnM(cau)f</italic>) is more conserved among the species where <italic>trnM(cau)e</italic> is present, than among species were it is absent (78.1% vs. 67.8% pairwise identity on average).</p> <p>Our previous analysis discovered several cases of tRNA gene recruitment in <italic>Axinella corrugata</italic>
##REF##15734570##[27]##. The more expanded dataset of demosponge mitochondrial tRNA genes assembled for this study revealed several additional instances of tRNA gene recruitment in demosponge mtDNA (to be described elsewhere).</p> <title>Intergenic regions and repeats</title> <p>The combined size of non-coding regions in the 22 demosponge mtDNAs analyzed in this study varied from 371 bp in <italic>Geodia neptuni</italic> to 6077 bp in <italic>Axinella corrugata</italic> or from 2 to 24% of the total genome size. In contrast to bilaterian animals, the distribution of non-coding nucleotides was more even in demosponge mtDNA, with the largest intergenic region usually containing &lt;15%, and at most 39% (in <italic>Iotrochota birotulata</italic>), of all non-coding nucleotides. We found little conservation in the position of the largest intergenic regions among the sampled genomes, even for the species that share identical gene arrangements, such as <italic>Chondrilla nucula</italic> and <italic>Halisarca dujardini</italic>, <italic>Geodia neptuni</italic> and <italic>Cinachyrella kuekenthali</italic>, and <italic>Hippospongia lachne</italic> and <italic>Vaceletia</italic> sp. Furthermore, we detected little sequence conservation either among individual regions within each mtDNA or between identically located non-coding regions in different species, except for the presence of repetitive elements in some genomes, as described below.</p> <p>Multiple repetitive elements were found in several analyzed genomes. Repeats larger than 100 bp were found only in <italic>Vaceletia</italic> sp., with the two biggest repetitive elements (229 bp) located in the intergenic regions that flank <italic>nad2</italic>, while 20–100 bp repeats were discovered in multiple species. The most abundant repeats were found in <italic>Vaceletia</italic> sp., <italic>Igernella notabilis</italic>, <italic>Ephydatia muelleri</italic>, and <italic>Axinella corrugata,</italic> where they have been located in most intergenic regions, as well as in ribosomal RNA genes and even some protein coding genes. The presence of repeated elements was very sporadic in respect to phylogeny, with repeats often present/absent in closely related species. Overall, repeats were very rare in sampled species from G0, G2 and G3, but more common in G1 and G4.</p>
<title>Discussion</title> <p>Our analysis of 22 complete mtDNA sequences representing all 14 orders of demosponges revealed both remarkable conservation and also an extensive diversity in mitochondrial genome organization within this group. Among the features shared among all sampled demosponge mitochondrial genomes are compact organization of the genetic material, similar gene content, well conserved structures of encoded tRNAs, a minimally modified genetic code for mitochondrial translation, and the absence of a single large “control” region characteristic of mtDNA in bilaterian animals. Genomic features that showed substantial variation include the number of tRNA genes, rRNA structures, the presence/absence of introns, and gene arrangements. In particular, two groups clearly stand out in our analysis with respect to their genome organization: G0 (order Homosclerophorida) and G1 (orders Dictyoceratida, Dendroceratida, and Verticillitida).</p> <p>As reported previously, the mitochondrial genome of the homosclerophorid <italic>Oscarella carmela</italic> contains 44 genes–the largest complement of genes in animal mtDNA–including <italic>tatC,</italic> a gene for subunit C of the twin arginine translocase that has not been found in any other animal mtDNA, and genes for 27 tRNAs ##REF##17090697##[14]##. By contrast, the mtDNA sequence of the homosclerophorid <italic>Plakortis angulospiculatus</italic> determined for this study contains only 20 genes and lacks <italic>tatC</italic> as well as 19 of the 25 tRNA genes typical for demosponges. Other differences between these two genomes include distinct gene arrangements (only 4 shared gene boundaries) and the presence of two group I introns in <italic>P. angulospiculatus cox1</italic>. Furthermore, the estimated genetic distances between these two species are greater than those between many orders of demosponges, indicating an ancient radiation and the presence of extensive genetic diversity within this group.</p> <p>Mitochondrial genomes of the three species within the G1 group are also unusual. These genomes lack all but two tRNA genes (for methionine and tryptophan tRNAs)–a feature previously associated with cnidarian mtDNA ##UREF##6##[28]##. Furthermore, this is the only group of demosponges where a significant acceleration in the rates of mitochondrial sequence evolution has been detected. There appears to be no causal connection between these two observations, as the loss of all but two tRNA genes is shared by all three species in the group, while the accelerated sequence evolution is much more pronounced in Dictyoceratida and Verticillitida. The retention of <italic>trnW(uca)</italic> and <italic>trnM(cau)</italic> as the only tRNA genes in the genome supports our previous inference ##REF##16945488##[29]## that these genes are difficult to replace because of the unique role of their products in mitochondrial translation:  is used for the initiation of mitochondrial translation with formylmethionine ##REF##5760639##[30]## while  must translate the TGA in addition to the TGG codons as tryptophan. The presence of such constraints can cause a parallel genomic evolution in independent lineages.</p> <p>An unusual mitochondrial genome has been previously reported for the haplosclerid demosponge <italic>Amphimedon queenslandica</italic>
##REF##17053047##[15]##. This genome lacks <italic>atp9</italic> and at least seven tRNA genes, contains deletions in several protein coding genes, and displays accelerated rates of sequence evolution in both protein and RNA genes. Our analysis of three additional species from the same order, <italic>Amphimedon compressa</italic>, <italic>Callyspongia plicifera,</italic> and <italic>Xestospongia muta</italic>, found no similar features in the latter taxa. These results most likely indicate that <italic>A. queenslandica</italic> mitochondrial genome has undergone an unusual evolution and is a poor representative of the G3 group, although incorporation of nuclear sequences, such as nuclear Numts ##REF##11369110##[31]##, in the mtDNA assembly cannot be ruled out. Given that <italic>A. queenslandica</italic> has become a model system for the study of demosponge biology, the evolution of its unusual mtDNA should be investigated in more details.</p> <p>Another interesting result that came out of this study is the discovery of two group I introns in <italic>cox1</italic> of <italic>P. angulospiculatus</italic>. Several lines of evidence, including phylogenetic analysis, the identical location in <italic>cox1</italic>, a similar extent of genetic divergence to their host genes, and the presence of TGA codons at the same position, support the vertical evolution of one of these introns from the common ancestor shared not only with <italic>Tetilla</italic> sp. (order Spirophorida), but also with scleractinian corals. This in turn suggests that the absence of this intron in most demosponge lineages is due to massive parallel loss. While examples of such losses are well known in nuclear genomes ##REF##17470438##[32]##–##REF##15231741##[37]##, an interesting question posed by this result is why mitochondrial introns are retained so scarcely in demosponges but so commonly in cnidarians?</p> <p>Finally, this study is interesting in what we did not find–any structures and/or sequences potentially involved in the maintenance and expression of mtDNA. Obviously, replication and transcription initiation/termination signals do exist in these genomes, but they were not detected by our comparative genomic analysis. Further data collection and experimental work will be essential to elucidate the mechanisms of these processes in demosponge mitochondria.</p>
<p>Conceived and designed the experiments: DVL. Performed the experiments: XW. Analyzed the data: XW DVL. Contributed reagents/materials/analysis tools: DVL. Wrote the paper: XW DVL.</p> <p>Two major transitions in animal evolution–the origins of multicellularity and bilaterality–correlate with major changes in mitochondrial DNA (mtDNA) organization. Demosponges, the largest class in the phylum Porifera, underwent only the first of these transitions and their mitochondrial genomes display a peculiar combination of ancestral and animal-specific features. To get an insight into the evolution of mitochondrial genomes within the Demospongiae, we determined 17 new mtDNA sequences from this group and analyzing them with five previously published sequences. Our analysis revealed that all demosponge mtDNAs are 16- to 25-kbp circular molecules, containing 13–15 protein genes, 2 rRNA genes, and 2–27 tRNA genes. All but four pairs of sampled genomes had unique gene orders, with the number of shared gene boundaries ranging from 1 to 41. Although most demosponge species displayed low rates of mitochondrial sequence evolution, a significant acceleration in evolutionary rates occurred in the G1 group (orders Dendroceratida, Dictyoceratida, and Verticillitida). Large variation in mtDNA organization was also observed within the G0 group (order Homosclerophorida) including gene rearrangements, loss of tRNA genes, and the presence of two introns in <italic>Plakortis angulospiculatus</italic>. While introns are rare in modern-day demosponge mtDNA, we inferred that at least one intron was present in <italic>cox1</italic> of the common ancestor of all demosponges. Our study uncovered an extensive mitochondrial genomic diversity within the Demospongiae. Although all sampled mitochondrial genomes retained some ancestral features, including a minimally modified genetic code, conserved structures of tRNA genes, and presence of multiple non-coding regions, they vary considerably in their size, gene content, gene order, and the rates of sequence evolution. Some of the changes in demosponge mtDNA, such as the loss of tRNA genes and the appearance of hairpin-containing repetitive elements, occurred in parallel in several lineages and suggest general trends in demosponge mtDNA evolution.</p>
<title>Supporting Information</title>
<p>We thank Alexander Ereskovsky, Michelle Kelly, Franz Lang, Sally Leys, Shirley Pomponi, Antonio Solé-Cava, Robert Thacker, and Gert Wörheide for samples contributed to this project; Natalia Frishman, Megan Rasmuson, Brandi Sigmon and Zhiyong Shao, for assistance in laboratory work; Karri Haen and Ehsan Kayal for their help with earlier versions of this manuscript.</p>
<fig id="pone-0002723-g001" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002723.g001</object-id><label>Figure 1</label><caption><title>Phylogenetic analysis of demosponge relationships using mitochondrial genomic data.</title><p>Posterior majority-rule consensus tree obtained from the analysis of 2,558 aligned amino acid positions under the CAT+F+Γ model is shown. Other methods of phylogenetic reconstruction produced similar topologies ##UREF##4##[11]##. The numbers at each node are Bayesian posterior probabilities. Nodes with ≥95% support are marked with an asterisk. For simplicity, non-demosponge clades were collapsed to triangles. The full tree is presented in ##SUPPL##0##Figure S1##.</p></caption></fig> <fig id="pone-0002723-g002" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002723.g002</object-id><label>Figure 2</label><caption><title>The size (A) and gene content (B) of demosponge mtDNA. Demosponge species are subdivided into five major groups (G0–G5).</title><p>Selected species from other animal groups and the outgroup <italic>Monosiga brevicollis</italic> are included for comparison. Species are abbreviated as following: <italic>mb</italic>, <italic>M. brevicollis</italic>; <italic>ta</italic>, <italic>Trichoplax adhaerens</italic>; <italic>ms</italic>, <italic>Metridium senile</italic>; <italic>hs</italic>, <italic>Homo sapiens</italic>; <italic>oc</italic>, <italic>Oscarella carmela</italic>; <italic>pa</italic>, <italic>Plakortis angulospiculatus</italic>; <italic>hl</italic>, <italic>Hippospongia lachne</italic>; <italic>in</italic>, <italic>Igernella notabilis</italic>; <italic>vs</italic>, <italic>Vaceletia</italic> sp.; <italic>af</italic>, <italic>Aplysina fistularis</italic>; <italic>cn</italic>, <italic>Chondrilla nucula</italic>; <italic>hd</italic>, <italic>Halisarca dujardini</italic>; <italic>ac</italic>, <italic>Amphimedon compressa</italic>; <italic>aq</italic>, <italic>Amphimedon queenslandica</italic>; <italic>cp</italic>, <italic>Callyspongia plicifera</italic>; <italic>xm</italic>, <italic>Xestospongia muta</italic>; <italic>as</italic>, <italic>Agelas schmidti</italic>; <italic>ck</italic>, <italic>Cinachyrella kuekenthali</italic>; <italic>ef</italic>, <italic>Ectyoplasia ferox</italic>; <italic>em</italic>, <italic>Ephydatia muelleri</italic>; <italic>gn</italic>, <italic>Geodia neptuni</italic>; <italic>to</italic>, <italic>Topsentia ophiraphidites</italic>; <italic>ib</italic>, <italic>Iotrochota birotulata</italic>; <italic>pw</italic>, <italic>Ptilocaulis walpersi</italic>; <italic>ax</italic>, <italic>Axinella corrugata</italic>; <italic>te</italic>, <italic>Tethya actinia</italic>.</p></caption></fig> <fig id="pone-0002723-g003" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002723.g003</object-id><label>Figure 3</label><caption><title>Nucleotide composition of mtDNA in five major groups of demosponges.</title><p>(A) A+T content; (B) AT-skew; (C) GC-skew. The values are shown for the sense (non-template) strand of the whole genome (genome), its concatenated genetic components (protein genes, rRNA genes, and tRNA genes), 3<sup>rd</sup> codon positions in protein genes, and for the corresponding strand in intergenic regions. Colored bars indicate the mean value for each group of demosponges; error bars show standard deviation.</p></caption></fig> <fig id="pone-0002723-g004" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002723.g004</object-id><label>Figure 4</label><caption><title>Mitochondrial gene arrangements in demosponges.</title><p>Protein and rRNA genes (larger boxes) are: <italic>atp6</italic>, <italic>8-9</italic>–subunits 6, 8 and 9 of the F<sub>0</sub> ATPase, <italic>cox1-3</italic>–cytochrome <italic>c</italic> oxidase subunits 1-3, <italic>cob</italic>–apocytochrome <italic>b</italic> (<italic>cob</italic>), <italic>nad1-6</italic> and <italic>nad4L</italic>–NADH dehydrogenase subunits 1-6 and 4L, <italic>rns</italic> and <italic>rnl</italic>–small and large subunit rRNAs, <italic>tatC</italic>–twin-arginine translocase component C. tRNA genes (smaller boxes) are abbreviated using the one-letter amino acid code. The two arginine, isoleucine, leucine, and serine tRNA genes are differentiated by subscripts with <italic>trnR(ucg)</italic> marked as <italic>R<sub>1</sub></italic>, <italic>trnR(ucu)</italic>–as <italic>R<sub>2</sub></italic>, <italic>trnI(gau)</italic>–as <italic>I<sub>1</sub></italic>, <italic>trnI(cau)</italic>–as <italic>I<sub>2</sub></italic>, <italic>trnL(uag)</italic>–as <italic>L<sub>1</sub></italic>, <italic>trnL(uaa)</italic> as <italic>L<sub>2</sub></italic>, <italic>trnS(ucu)</italic>–as <italic>S<sub>1</sub></italic>, and <italic>trnS(uga)</italic>–as <italic>S<sub>2</sub></italic>. All genes are transcribed from left to right except those underlined to indicate an opposite transcriptional orientation. Genes are not drawn to scale and intergenic regions are not shown.</p></caption></fig> <fig id="pone-0002723-g005" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002723.g005</object-id><label>Figure 5</label><caption><title>Relative rates of evolution of individual species (A) and individual genes (B).</title><p>Rates are estimated by average genetic distances to orthologous genes from four cnidarians. Each boxplot represents data for 13 individual genes in (A) and 22 demosponge species in (B). Lower horizontal bar, non-outlier smallest observation; lower edge of rectangle, 25 percentile; central bar within rectangle, median; upper edge of rectangle, 75 percentile; upper horizontal bar, non-outlier largest observation; open circle, outlier.</p></caption></fig> <fig id="pone-0002723-g006" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002723.g006</object-id><label>Figure 6</label><caption><title>Inferred secondary structure of <italic>Hippospongia lachne</italic> mitochondrial small subunit RNA in comparison to that of <italic>Oscarella carmela.</italic>
</title><p>The helices are numbered in boldface as in Brimacombe <italic>et al.</italic>
##REF##7607205##[52]##. Structural regions present in <italic>O. carmela</italic> srRNA but absent in <italic>H. lachne</italic> srRNA are shown in red.</p></caption></fig> <fig id="pone-0002723-g007" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002723.g007</object-id><label>Figure 7</label><caption><title>Secondary structures and consensus sequences of demosponge mitochondrial tRNAs.</title><p>The secondary structure of each type of tRNAs was folded based on sequence and structure alignment. Nucleotides in uppercase letters indicate &gt;90% sequence conservation, lowercase letters indicate &gt;75% sequence conservation, and the dots represent &lt;75% conservation.</p></caption></fig>
<inline-formula></inline-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <inline-formula></inline-formula>
<supplementary-material content-type="local-data" id="pone.0002723.s001"><label>Figure S1</label><caption><p>(0.29 MB EPS)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pone.0002723.s002"><label>Figure S2</label><caption><p>(0.24 MB EPS)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pone.0002723.s003"><label>Figure S3</label><caption><p>(0.69 MB EPS)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pone.0002723.s004"><label>Table S1</label><caption><p>(0.04 MB XLS)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pone.0002723.s005"><label>Table S2</label><caption><p>(0.04 MB XLS)</p></caption></supplementary-material>
<fn-group><fn fn-type="COI-statement"><p><bold>Competing Interests: </bold>The authors have declared that no competing interests exist.</p></fn><fn fn-type="financial-disclosure"><p><bold>Funding: </bold>The funding for this work came from the College of Liberal Arts and Sciences at Iowa State University.</p></fn></fn-group>
<graphic xlink:href="pone.0002723.g001"/> <graphic xlink:href="pone.0002723.g002"/> <graphic xlink:href="pone.0002723.g003"/> <graphic xlink:href="pone.0002723.g004"/> <graphic xlink:href="pone.0002723.g005"/> <graphic xlink:href="pone.0002723.g006"/> <inline-graphic xlink:href="pone.0002723.e001.jpg" mimetype="image"/> <inline-graphic xlink:href="pone.0002723.e002.jpg" mimetype="image"/> <inline-graphic xlink:href="pone.0002723.e003.jpg" mimetype="image"/> <graphic xlink:href="pone.0002723.g007"/> <inline-graphic xlink:href="pone.0002723.e004.jpg" mimetype="image"/> <inline-graphic xlink:href="pone.0002723.e005.jpg" mimetype="image"/>
<media xlink:href="pone.0002723.s001.eps"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pone.0002723.s002.eps"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pone.0002723.s003.eps"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pone.0002723.s004.xls"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pone.0002723.s005.xls"><caption><p>Click here for additional data file.</p></caption></media>
[{"label": ["6"], "element-citation": ["\n"], "surname": ["Sollas"], "given-names": ["WJ"], "year": ["1885"], "article-title": ["A Classification of the Sponges."], "source": ["Sci Proc R Dublin Soc"], "volume": ["5"], "fpage": ["12"]}, {"label": ["7"], "element-citation": ["\n"], "surname": ["Vacelet", "Boury-Esnault"], "given-names": ["J", "N"], "year": ["1995"], "article-title": ["Carnivorous Sponges."], "source": ["Nature"], "volume": ["373"], "fpage": ["333"], "lpage": ["335"]}, {"label": ["8"], "element-citation": ["\n"], "surname": ["Hooper", "Van Soest", "Hooper", "Van Soest"], "given-names": ["JNA", "RWM", "JNA", "RWM"], "year": ["2002"], "article-title": ["Class Demospongiae Sollas, 1885."], "source": ["Systema Porifera: a guide to the classification of sponges"], "publisher-loc": ["New York"], "publisher-name": ["Kluwer Academic/Plenum Publishers"], "fpage": ["15"], "lpage": ["18"]}, {"label": ["9"], "element-citation": ["\n"], "surname": ["van Soest", "Boury-Esnault", "Janussen", "Hooper"], "given-names": ["R", "N", "D", "J"], "year": ["2005"], "article-title": ["World Porifera database."], "comment": ["Available online at "], "ext-link": ["http://www.vliz.be/vmdcdata/porifera/"]}, {"label": ["11"], "element-citation": ["\n"], "surname": ["Lavrov", "Kelly", "Wang"], "given-names": ["DV", "M", "X"], "year": ["2008"], "article-title": ["Reconstructing Ordinal Relationships in the Demospongiae Using Mitochondrial Genomic Data."], "source": ["Mol Phylogenet Evol"], "comment": ["doi:10.1016/j.ympev.2008.05.014"]}, {"label": ["17"], "element-citation": ["\n"], "surname": ["Pisera"], "given-names": ["A"], "year": ["2006"], "article-title": ["Palaeontology of sponges\u2014a review."], "source": ["Can J Zool"], "volume": ["84"], "fpage": ["242"], "lpage": ["261"]}, {"label": ["28"], "element-citation": ["\n"], "surname": ["Beagley", "Macfarlane", "Pont-Kingdon", "Okimoto", "Okada"], "given-names": ["CT", "JL", "GA", "R", "N"], "year": ["1995"], "article-title": ["Mitochondrial genomes of Anthozoa (Cnidaria)."], "source": ["Progress in Cell Research"], "volume": ["5"], "fpage": ["149"], "lpage": ["153"]}]
PLoS One. 2008 Jul 16; 3(7):e2723