Sponge Classes

Analyses of the seven nuclear housekeeping gene dataset provide strong support for each of the four major clades of sponges assigned the rank of class (Calcarea, Demospongiae, Hexactinellida, and Homoscleromorpha). Because we did not include non-metazoan outgroups our results cannot be used to assess sponge monophyly. Concerning the relationships among the four sponge classes, support is generally poor. Our tree does not recover Silicea^p (Demospongiae^p+Hexactinellida^p), which has been supported in a great deal of other works based on disparate datasets [4], [18], [19], [28], [54], but instead places Calcispongiae^p with Hexactinellida^p (Figs. 1–2), most likely erroneously with low support (bs = 74%; pp = 0.68). Relationships within Calcispongiae^p and Hexactinellida^p are consistent with previous analyses [54]–[56]. As designed, our analyses do not provide any basis for inferring relationships among the sponge classes (as they do not include non-metazoan outgroups), but rather elucidate phylogenetic relationships within Demospongiae^p (Figs. 1–2).

Major Demosponge Clades

Hypotheses derived from our analyses of nuclear housekeeping gene data (Figs. 1–2) are fairly consistent with the so-called “G clades” originally derived from analysis of ribosomal DNA data [17], and largely recovered by mitochondrial genome [43] and nuclear housekeeping gene data [18]. G1 and G2 correspond to Keratosa^p and Myxospongiae^p, respectively, following the names of Borchiellini et al. [17]. One key difference between the results of these studies concerns the placement of the clade containing all freshwater sponges, Spongillida^p, phylogenetically defined in Cárdenas et al. [40]. Traditionally, these sponges were classified as the suborder Spongillina within the order Haplosclerida. However, ribosomal DNA and mitochondrial genome data suggested that Spongillida^p falls as the earliest diverging lineage of the “G4” clade. Sperling et al. [18] found a similar clade, for which they provided a phylogenetic definition and the name Democlavia^p ( = subclass Heteroscleromorpha of Cárdenas et al. [40]), with the exception that Spongillida^p was found as the sister group of the marine haplosclerids. The marine haplosclerid taxa have consistently been shown to be a well-supported clade that has recently been phylogenetically defined and named Haploscleromorpha^p [40].

This study finds strong support at nearly all deep nodes within Demospongiae^p (Figs. 1–2), even with our more diverse taxon sampling. The clear distinction of these clades indicates that the divergence among these groups is likely ancient [4]. Thus, future genomic exploration within Demospongiae^p will be guided by these emerging phylogenetic results so as to make best use of the comparative method. To be especially useful for rank-based taxonomy and nomenclature, type species within genera and type genera within families (e.g., our sampling of Spongia officinalis, Halisarca dujardini, and Desmacella pumilio) should be targeted whenever possible. Also, to the extent possible, type species should be collected from their respective type localities for maximum taxonomic and nomenclatural utility. For phylogenetic nomenclature, ‘specifiers’ (i.e., species, specimens or apomorphies used in PhyloCode definitions) should be targeted. Of course, when species are used as specifiers (which has so far usually been the case for poriferan names), their name-bearing type specimens are de facto specifiers (PhyloCode, Note 13.2.2.).

Nuclear housekeeping gene data strongly support an as yet unnamed clade containing the groups of demosponges with silica-mineralized skeletons: Democlavia^p, Haploscleromorpha^p, and Spongillida^p (Figs. 1–2), in accordance with other analyses of ribosomal genes [17], complete mitochondrial genomes [43], and a smaller dataset of nuclear housekeeping genes [18]. Our ML and Bayesian analyses provide equivocal results concerning the phylogenetic relationships of Keratosa^p and Myxospongiae^p. A sister group relationship between Keratosa^p and Myxospongiae^p has been suggested, with only modest support, based on analyses of 18S rRNA genes [17] and complete mitochondrial genomes [43] but has also been contradicted by earlier Bayesian analyses of nuclear housekeeping genes [4], [18], [44]. Our ML topology (Fig. 1) shows Keratosa^p and Myxospongiae^p [which both lack mineralized skeletons (with the exception of siliceous microscleres in Chondrilla within Myxospongiae^p: Chondrosida)] as a clade that is sister to the mineralized sponges. In contrast, the Bayesian analysis (Fig. 2) identifies Myxospongiae^p as the earliest diverging clade of Demospongiae^p, and shows Keratosa^p as the sister group to the mineralized groups. It is important to note, however, that all site-homogeneous models (LG and GTR) display the Keratosa^p+Myxospongiae^p clade, while the site-heterogeneous CAT and CAT-GTR models (which fit the data better) support Myxospongiae^p as the sister group of all the other demosponges. Thus, model selection is responsible for this disagreement. Because the best fitting models suggest Myxospongiae^p is sister to the remaining demosponges, the contradicting results obtained using LG, GTR and WAG (Keratosa^p+Myxospongiae^p) are likely artifactual.

Keratosa^p

This clade is composed of members of the demosponge orders Dictyoceratida and Dendroceratida. Our sampling includes members of five of the six families: Dysideidae, Irciniidae, Spongiidae and Thorectidae in the former, Dictyodendrillidae in the latter. Ribosomal data [17] indicate that Dendroceratida is monophyletic, but our results rely on a single genus (Igernella) so we cannot support or refute that result. The nuclear housekeeping gene data also fail to provide support for the monophyly of Dictyoceratida, a result that has also been obtained through the analysis of ribosomal data [35], [57]. We have conflicting results concerning Dictyoceratida, with our ML-topology (Fig. 1) suggesting that dendroceratids are derived from within a paraphyletic Dictyoceratida and the Bayesian tree having a poorly supported monophyletic Dictyoceratida. The key taxon, from the perspective of this analysis, is the representative of Dysideidae. All the other dictyoceratids in our study, representing Irciniidae, Spongiidae, and Thorectidae, always form a well-supported clade. It is interesting to note that when the worst performing markers (CAT, EF1A, and ATPB) are sequentially removed from analysis, Dictyoceratida, including our representative of Dysideidae, forms a monophyletic group with strong support (Figs. S15, S16, S17).

Myxospongiae^p

Members of the orders Chondrosida and Verongida make up Myxospongiae^p. Our sampling includes both families of Chondrosida (Chondrillidae and Halisarcidae), the latter of which was previously placed in its own order Halisarcida (e.g., [58]). Within Verongida, just one of the four families of Verongida (i.e., Aplysinidae) is sampled. With the present taxon sampling, our analyses support monophyly of Chondrosida, a result not obtained by some analyses of ribosomal data [17], [59], but found in others [35], [60]. However, our analysis lacks a representative of Chondrosia, which has proven to be a difficult taxon in relation to the question of Chondrosida monophyly [17], [59]. Similarly lacking a representative of the problematic Chondrosia, an analysis of complete mitochondrial genome data also supports a monophyletic Chondrosida^p [43], which has nevertheless recently been given a phylogenetic definition [40]. Within Verongida, nuclear housekeeping genes support monophyly of Aplysinidae, for which we were able to sample each of its component genera (Figs. 1–2). Relationships among the three aplysinid genera (Verongula, Aplysina, and Aiolochroia), however, are not well supported. Based on ribosomal data, Erwin and Thacker [61] found that Aplysinidae is not monophyletic because Verongula grouped with members of Pseudoceratinidae and members of Aiolochroia grouped with Ianthellidae and Aplysinellidae. The absence of pseudoceratinids, ianthellids and aplysinellids from our samples prevents our analyses from testing these hypotheses, but if Erwin and Thacker's [61] findings are true, they would suggest that our sampling represents a more disparate group of Verongida (Aplysina in Aplysinidae and Verongula in Pseudoceratinidae) than is suggested by current taxonomy (Aplysina and Verongula in Aplysinidae). Indeed, this phylogenetic result (i.e., that Aplysina and Verongula belong to distinct families) was recently verified with mitochondrial and nuclear markers by Erpenbeck et al. [59].

Haploscleromorpha^p & Spongillida^p

From a broad perspective, one of the most important outstanding questions in demosponge phylogenetics is the phylogenetic placement of the freshwater sponges, Spongillida^p, which is phylogenetically defined in Cárdenas et al. [40]. Traditional taxonomy based on morphology [62] and earlier analyses of nuclear housekeeping genes [18] suggest a close relationship between Spongillida^p and the marine haplosclerids, Haploscleromorpha^p. In contrast, both mitochondrial genome and ribosomal data suggest that Spongillida^p is sister to the rest of the Democlavia^p [17], [35], [43], [63]. The results here, for the most part, agree with the former hypothesis and specifically indicate that Spongillida^p is the earliest diverging lineage of the traditional order Haplosclerida (with high support, Figs. 1–2). An exception to this result is one of the single gene analyses (ALD, Fig. S8), which found Spongillida^p branching among democlaviid taxa, albeit with no support. Limited taxon sampling, and in particular, the fact that our analyses do not include any representatives of the democlaviid family Scopalinidae (which was recently suggested by Morrow et al. [38] to have a close relationship to the freshwater sponges), could explain these contradicting results. In any event, it is fairly clear that Spongillida^p is a distinct lineage from the marine haplosclerids.

Our sampling within Haploscleromorpha^p represents five of the six accepted families. Monophyletic haplosclerid suborders Petrosina and Haplosclerina were not recovered (although support values are somewhat low at some of the deeper branches of the clade), corroborating the results of McCormack et al. [64] and Redmond et al. [35], [37]. Not surprisingly, given that studies with denser taxon sampling have shown widespread polyphyly of subtaxa within this group [35], [37], [65], we find both Petrosiidae and Niphatidae to be polyphyletic. Even at the generic level, Amphimedon (Niphatidae) is revealed to be polyphyletic. Amphimedon queenslandica, whose genome has been sequenced [66], clusters with Callyspongia vaginalis (Callyspongiidae) with high support, suggesting that the taxonomy of this important model organism remains confused, corroborating evidence from ribosomal data [35], [37].

Democlavia^p

Democlavia^p is the most species-rich (roughly 75% of demosponge species; [38]) and diverse of the major demosponge clades, and includes the traditional orders Agelasida, Astrophorida, Hadromerida, Halichondrida, Poecilosclerida, and Spirophorida [48], several of which are already thought to not be monophyletic (as discussed below). As such, the systematics of Democlavia^p presents many challenges, but important breakthroughs are being made in understanding the phylogeny of this clade based on increasingly taxon-rich analyses of ribosomal RNA and mitochondrial CO1 data [38]. Our nuclear housekeeping gene dataset and analyses provide an opportunity to test hypotheses arising from these alternative sets of data and suggest new hypotheses where previous results have provided no resolution.

Our analyses reveal a well-supported clade containing members of Astrophorina and Spirophorina (suborder designations for these taxa, following [40]), including our only sampled lithistid (Microscleroderma sp. nov.). Other analyses of ribosomal and mitochondrial data have revealed the same clade [17], [35], [42], [67]–[69], the phylogenetically defined Tetractinellida^p [17], [40]. Although modest in support, our analyses always suggest that Tetractinellida^p is sister to the remaining members of Democlavia^p. Our sampling of sub-order Astrophorina includes two of the six families, Ancorinidae (Dercitus, recently transferred from Pachastrellidae by Cárdenas et al. [70]) and Geodiidae (Geodia tumulosa and Geostelletta^p fibrosa), as well as an incertae sedis taxon, Characella aff. connectens, which was also formerly assigned to family Pachastrellidae. The latter three species form a well-supported clade, but no specific position for our representative of Ancorinidae within Tetractinellida^p is supported (Figs. 1–2). The family Pachastrellidae sensu Maldonado [71] is based on a plesiomorphic character (streptasters; [70]) so it is no surprise that our results confirm that Characella and Dercitus do not have an especially close relationship.

Our analyses include two representatives of Spirophorina – Cinachyrella sp., representing the family Tetillidae, and the lithistid Microscleroderma sp. nov., representing the family Scleritodermidae – but there is no support for the group being monophyletic. The lithistids are a taxonomically rich group sharing a common growth form (skeleton of interlocked desmas), with 13 recognized families. Lithistids have always presented taxonomic challenges from morphological perspectives (see 72) and the redistribution of its members to different sponge clades has been proposed for quite some time [72], [73] and continues [40]. In this vein, the lithistid family Desmanthidae appears to be closely related to Dictyonellidae [38]. The presence of sigmaspires in Scleritodermidae [72] is consistent with this group being reallocated to Spirophorina within Tetractinellida^p [40].

Another well-supported alliance of taxa includes most members of order Poecilosclerida that we have sampled, specifically representatives of Coelosphaeridae, Crambeidae, Hymedesmiidae, Microcionidae, Mycalidae, and Tedaniidae (Figs. 1–2). Monophyly of Poecilosclerida has been found in several analyses of ribosomal data [17], [35], [42], [74], but more recent studies with greater taxon sampling have shown the group to be polyphyletic [38], [75], as found here. Morrow et al. [38] demonstrated that the families Desmacellidae and Raspailiidae should be removed from Poecilosclerida. Our results support this action, as our representatives of these families branch deeper within Democlavia^p (Figs. 1–2). Unfortunately, these data do not provide strong support for relationships within this poecilosclerid group, which remains the most species-rich order and therefore one of the more challenging clades within Demospongiae^p.

The sister group to Poecilosclerida (sensu 38) consists of most of our sampled hadromerids as well as the family Halichondriidae from the order Halichondrida. A similar relationship was derived in Morrow et al. [38]. Within this clade, three hadromerids, Cliona (Clionaidae), Placospongia (Placospongiidae), and Spirastrella (Spirastrellidae) form a well-supported clade. In turn, this clade is revealed to have a relatively well-supported relationship with the families Halichondriidae and Suberitidae. The latter two families, currently classified within Halichondrida and Hadromerida, respectively, have long been known to have a close relationship [27]. Interestingly, the hadromerid Tethya (Tethyidae) consistently branches with this alliance of Suberitidae, Halichondriidae, and the hadromerids (representing Clionaidae, Placospongiidae and Spirastrellidae) albeit with limited support. One other hadromerid in our analysis, Polymastia tenax, falls outside this clade, a peculiar result given that Polymastiidae is considered among the “core” components of Hadromerida [76]. In the 28S-based analysis of Morrow et al. [38], Polymastiidae emerged as a distinct clade, sister to Suberitidae plus Halichondriidae but with low support, whereas their analysis of CO1 data recovered a clade with Polymastiidae sister to the hadromerid families Tethyidae, Hemiasterellidae, and Clionaidae, but again with only low support.

The monophyly of Agelasida^p is well supported. This result is obtained only after taking into account recent findings made by Gazave et al. [36], who provided a phylogenetic definition of the clade, and corroborated by Morrow et al. [38]. In light of polyphyly of Axinella (order Axinellida), Gazave et al. [36] erected the taxon Cymbaxinella^p for those species, including Axinella corrugata sampled here, with a close relationship to Agelas (family Agelasidae). With broader taxon sampling, Morrow et al. [38] established the new family Hymerhabdiidae for this same clade within Agelasida^p. In contrast with this study [38], however, nuclear housekeeping gene data do not provide further support for a sister group relationship between Agelasida^p and the clade containing the core poecilosclerids, hadromerids and Halichondriidae. The only representative of order Axinellida in our analysis is Ectyoplasia; the species belongs to the family Raspailiidae, which was moved from Poecilosclerida to Axinellida by Morrow et al. [38]. That study [38] also found that representatives of Desmacellidae fell in two groups, a finding we also recovered given that Desmacella and Biemna did not exhibit a particularly close relationship. It is important to note that our analysis includes the type species of Desmacella. Nuclear housekeeping gene data provide modest support for a relationship between Desmacella and the family Dictyonellidae (Figs. 1–2).

Ethics Statement

In accordance with policy and legal requirements associated with specimens vouchered in the collections of the Smithsonian US National Museum of Natural History (NMNH), Harbor Branch Oceanographic Institute (HBOI), Harvard Museum of Comparative Zoology (MCZ), and Zoological Museum Bergen Norway (ZMBN), all collections involved in this study were obtained with all appropriate and relevant permits. Specifically, samples from Panama were collected under a Marine Collecting Permit provided by The Republic of Panama; samples from the State of Florida were collected under a Florida recreational resident saltwater fishing license issued from Florida Fish and Wildlife Conservation Commission; and one sample from Honduras was collected with the permission of Rosa del Carmen Garcia, Directora General de Pesca y Acuicultura. No permits were required to collect sponge specimens in US territorial waters outside state boundaries, the Catalan coast of Spain, Vancouver Island, Canada, or Norway.

Sample and sequence collection

Samples were collected from a variety of locations and stored as described below or obtained from frozen collections at the Harbor Branch Oceanographic Institute-Florida Atlantic University (Table 1; http://PorToL.org/NHK7data). To obtain RNA of sufficient quality and quantity, when possible, fresh material was collected and preserved via one of several methods. One involved placing fresh material in cold 75% ethanol with liquid changes occurring after 15 min, 1 hour and 4 hours. When available, material was also placed in RNAlater (Invitrogen), directly in TRIzol® (Invitrogen) reagent, following the manufacturer's instructions, or in liquid nitrogen. In most cases, the tissue placed directly in TRIzol® or frozen in liquid nitrogen yielded the highest quality and/or quantity of RNA. However, the most practical storage method in the field was 75% ethanol preservation and in most cases this was suitable for RNA extraction and subsequent polymerase chain reaction (PCR) amplifications from cDNA.

Following Sperling et al. [18], [44] total RNA was isolated using a one-step TRIzol® method (Invitrogen), and cDNA was synthesized from 1–2 µg RNA using RETROSCRIPT® (Ambion) reverse transcriptase using both random decamers and oligo dT primers, which were then pooled. PCR was used to amplify 7 nuclear-encoded genes: aldolase (ALD), ATP synthase beta chain (ATPB), catalase (CAT), elongation factor 1-alpha (EF1alpha), methionine adenosyltransferase (MAT), phosphofructokinase (PFK), and triose-phosphate isomerase (TPI). All primer sequences for initial PCR of housekeeping genes can be found in Sperling et al. [44]. In many cases, however, it was necessary to use nested PCR primers if amplification and re-amplification of housekeeping gene products was not possible. Table S3 provides primer sequences for nested amplifications of individual housekeeping genes. Primary or nested amplification products were cloned into PCR cloning vectors (pGEM®-T, Promega or TOPO TA®, Invitrogen) and individual clones were prepared for DNA sequencing using standard protocols.

After editing and trimming vector sequences with GENEIOUS [77], DNA sequences were assessed for gene and sponge identity via BLASTX or BLASTP queries [78], followed by preliminary single-gene phylogenetic analyses under the likelihood framework described below. The identification of likely paralogs followed standard procedures based on the generation of trees including all the members of each gene family that could be identified in GenBank (via reciprocal blasting). Within the context of these trees, paralogy groups were identified and only the sequences nesting within the selected orthology group were used. New sequences generated in this study have been submitted to GenBank (Table 1). Sequences are also available via the Porifera Tree of Life database (PorToL.org). In addition, voucher specimens for many of the sequences presented in Sperling et al. [18], [44] were examined, resulting in several instances of updated taxonomic identification and classification (Table 1).

Nucleotide sequences were translated and aligned using MUSCLE [79] and visualized in SEAVIEW (v. 4.3) [80]. In addition to the new sequences, the initial alignment included data for sponges that had already been published (Table 1). Also, five species for which transcriptome data exist were also added to the dataset. Both mRNA and cDNA from Corticium candelabrum, Petrosia ficiformis and Sycon coactum were obtained using protocols available in Riesgo et al. [81]. Sycon ciliatum and Leucosolenia complicata sequences are derived from current genome and transcriptome sequencing projects for these species [82] and Adamska, unpublished). De novo assemblies of the reads obtained with Illumina GA (Illumina, Solexa, USA) were built with CLC Genomics Workbench 4 (CLCbio, MA, USA). Local blasts against the contig lists generated were used to search for the housekeeping genes. Initially, 50 outgroup taxa representing Bilateria, Ctenophora, Cnidaria, Placozoa and non-metazoan Opisthokonta were included in the analyses. However, preliminary phylogenetic analyses, conducted as described below, indicated that inferred demosponge relationships were robust to outgroup choice and therefore outgroups in the final dataset were reduced to the cnidarian taxa (Acropora, Metridium and Nematostella) and the placozoan Trichoplax. Approximately 40 positions in the alignment were manually excluded from analyses because they represented insertions present in one or a small number (<5) of taxa.

Phylogenetic Analyses

For all gene trees we investigated the presence of significant clustering information using Maximum Likelihood Mapping [83] as implemented in Treepuzzle V. 5.2 [51]. The dataset was analyzed in both Bayesian and Maximum Likelihood (ML) frameworks. For the ML analyses, appropriate models of amino acid evolution were assessed using the Akaike Information Criterion (AIC), as implemented in ProtTest (v.2.4) [49]. The computing cluster of the Smithsonian's Laboratories of Analytical Biology was used to run the parallelized version of RaxML [84] to search for maximum likelihood (ML) topologies. We assumed the model that best fit our data according to the second-order AIC (AICc-1) with the exception that a proportion of invariant sites was not estimated (according to a recommendation in the RaxML manual). We also used RaxML to conduct bootstrap analyses (400 replicate searches) to assess nodal support. We searched for ML topologies using each gene separately as well as all genes combined. We analyzed the combined data a) assuming a single model for all the data and b) by assigning most appropriate models to each gene partition (mixed models).

Bayesian analyses were performed using the site-heterogeneous CAT-GTR+gamma in Phylobayes 3.3b [85]. This model was selected because Bayesian model selection, performed using 10-fold cross-validation [86], showed that CAT-GTR best fitted our dataset, outperforming CAT, GTR and LG. The considered models were: WAG, LG, GTR, CAT, and CAT-GTR (all models used a gamma correction to account for rate heterogeneity among sites). The CAT based models (in this case CAT and CAT-GTR [86]) are mixture models developed to better take into account site-specific features of protein evolution. These models are thus expected to fit the data better than homogeneous time reversible models like LG and GTR [86]. Indeed, CAT based models have previously been shown to fit amino acid datasets better than other models and they have been shown to be highly effective at reducing systematic biases, like long branch attraction, which are well known to be very pervasive in deep time phylogenetics. In Phylobayes two independent analyses were run for 30,000 cycles sampling every 100 points. The analyses were considered converged when the largest discrepancy observed across all bipartitions (i.e. the maxdiff statistics) dropped below 0.15, despite the Phylobayes manual's suggestion that a chain has reached convergence when maxdiff <0.3. Support values for the nodes recovered in the CAT-GTR analysis are expressed as posterior probabilities.

Comparisons were made between the different single-gene topologies and the Bayesian topology to the ML tree derived from the overall data. In addition, nodal differences were calculated, as measured by the root-mean-squared distance, in Topd (v.3.3) [87]. Taxa that were missing data for some genes were pruned from the combined tree prior to calculating nodal differences. Topd was also used to conduct randomization analyses to test whether similarities between the various topologies and the combined ML topology were not greater than expected by chance. Finally, further ML searches were conducted by sequentially excluding the three genes that subtend the trees that are most distant from the tree derived from the concatenated dataset, as measured by subtracting the random nodal difference from the actual nodal difference. To further investigate the extent to which the principal signal [53] in the single-genes corroborated the results of concatenated Bayesian and ML analyses, we performed a supertree analysis. The supertree was built using the Matrix Representation with Parsimony method [53]. Input trees used for this analysis were, for each gene, the 400 bootstrap trees derived (see above) under ML. This set of 2800 input trees was bootstrapped to generate 100 replicate datasets, each of which scored 2800 trees using the software CLANN [88]. For each bootstrapped dataset a bootstrap supertree was recovered and a majority rule consensus of the recovered bootstrap supertrees was built to estimate nodal support.

Finally, analyses were performed to test for tree reconstruction artifacts. More precisely we investigated the potential effect of long-branch attraction and compositional attraction on our results. We first investigated the effect of using alternative model of evolution on our results. We thus built trees (within a Bayesian framework) using models (WAG, LG, GTR, CAT, and CAT-GTR, each with a gamma correction) providing different levels of fit to the data and compared the trees we obtained. We tested whether the taxa in our dataset were compositionally heterogeneous performing a posterior predictive analysis (see for example [18]) of compositional heterogeneity using Phylobayes under the CAT-GTR model. The posterior predictive analysis indicated that several taxa displayed a biased composition of their sites. This, if not addressed, can cause compositional artifacts. To test whether our results were affected by compositional biases we performed two analyses. First we analysed (under CAT-GTR) a dataset from which all compositionally heterogeneous taxa were excluded. This experiment has the downside of excluding potentially important taxa. Accordingly, a second experiment was performed in which our dataset was recoded using the Dayhoff scheme. Dayhoff recoding can alleviate compositional artifact, and a posterior predictive analysis of our Dayhoff-recoded dataset was performed (under CAT-GTR) to evaluate whether further compositionally biased taxa remained after the application of Dayhoff recoding. Finally, our Dayhoff recoded dataset was analysed using both a site-homogeneous (GTR) and a site heterogeneous (CAT-GTR) model.

To test for the potential effect of long-branch attraction artifacts we identified fast evolving sites in our dataset using the program Tiger [89]. After that, sites that Tiger deemed as being fast evolving (bins 7 to 10) were excluded and the slowly evolving sites analysed in isolation. In addition to the site-stripping analysis, we also performed an analysis where all the outgroups to Demospongiae (including Hexactinellida) were removed.