Tackling the phylogenetic conundrum of Hydroidolina (Cnidaria: Medusozoa: Hydrozoa) by assessing competing tree topologies with targeted high-throughput sequencing

Bait design for targeted sequencing

Biotinylated RNA baits of 120 bp length were designed using 355 orthogroup partitions (loci) from the phylogenomic OF-PTP_75tx matrix (Kayal et al., 2018). This dataset contains the coding sequences from a broad taxonomic sample of cnidarians, including anthozoans, medusozoans, and endocnidozoans. Due to their large sequence divergences, endocnidozoans were excluded from the dataset used for RNA bait design. Nucleotide sequences were aligned per locus guided by their predicted amino acid translations using the Geneious aligner (version 9; Biomatters, Auckland, New Zealand). Several sequences were flagged after manual audit of alignments due to their apparent divergence and identified as non-homologous sequences after BLAST searches against NCBI’s GenBank. Following curation and removal of sequences shorter than 80 bp, putative repetitive elements in the remaining 5,523 sequences were soft-masked using the cross_match search option against the Cnidaria repeats database in RepeatMasker (version 4.06; Smit, Hubley & Green, 2013-2015).

Next, sequence stretches of 10 or fewer ambiguous nucleotides (N) were replaced with thymine (T) repeats to allow bait design across short regions of undetermined sequences. Similarly, short sequences less than 120 bp were padded with T repeats to allow for 120 bp long baits to be mapped against reference sequences. Baits were designed at approximately 2.5× tiling density with some 50 bp spacing between the start of neighboring baits, yielding 38,102 candidate baits. Each candidate bait was verified against six cnidarian genome assemblies (Table 1) using BLAST searches to evaluate the bait’s fit to the reference. Hits with a length greater than 45 bp in length and identity greater than 75% were retained for further analysis. Melting temperatures (Tm; defined as the temperature at which 50% of molecules hybridize) were estimated for each BLAST hit assuming standard Mybaits (Mycroarray, Arbor Biosciences, Ann Arbor, MI, USA) hybridization conditions. Based on the distribution of inferred Tms, baits of moderate or higher specificity were retained. That is, candidate baits with at most 10 BLAST hits between 62.5 °C to 65 °C, two BLAST hits above 65 °C, and fewer than 2 baits on each flanking region were retained, yielding a total of 37,546 baits (Supplemental File 1).

Target capture and sequencing

Genomic DNA was extracted from ethanol-preserved tissue samples (Table 2) using a standard organic phenol-chloroform extraction protocol (Green & Sambrook, 2012). Extracted DNA was quantified with a Qubit 4 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). 50 µl of DNA in TE buffer were chilled to 4 °C in the water bath of a Q800 Sonicator (Qsonica, Newton, CT, USA) and acoustically sheared for nine minutes using an amplitude of 25 with sonication pulses of 15s on/15s off. Illumina sequencing libraries were constructed from sheared DNA samples using an NEBNext Ultra II DNA library preparation kit (New England Biolabs, Ipswich, MA, USA) with dual indexes for multiplexing following the manufacturer’s protocol. After library amplification and magnetic bead purification using Ampure beads (Promega, Madison, WI, USA), amplicons longer or shorter than approximately 200 bp were removed using a BluePippin size select gel eletrophoresis system (Sage Science, Beverly, MA, USA). Concentrations of size-selected libraries were equilibrated, followed by pooling of libraries three to four samples deep for target enrichment.

Hybridization of RNA baits to pooled libraries followed the Mybaits version 3.02 protocol (Mycroarray, Arbor Biosciences, Ann Arbor, MI, USA) with the following modifications. After initial denaturation and blocking with Illumina adapter-specific oligonucleotides, RNA baits were allowed to hybridize for 19 h at 65 °C, 19 h at 60 °C, and 10 h at 55 °C. We used this touchdown procedure on the newly developed and untested bait set in an effort to increase on-target specificity while allowing for sensitivity of reactions. Following hybridization, biotinylated baits were bound to streptavidin-coated magnetic beads (Dynabeads MyOne Streptavidin; Thermo Fisher Scientific, Waltham, MA, USA), followed by stringent washing to remove unbound DNA library molecules. Captured libraries were amplified while bound to beads using KAPA HiFi DNA polymerase and HotStart ReadyMix (Roche, Basel, Switzerland) following the manufacturers protocol. The annealing temperature during the 15 amplification cycles was 60 °C. PCR reactions were cleaned using Ampure magentic beads (Thermo Fisher Scientific, Waltham, MA, USA), washing twice with 80% ethanol. DNA concentrations were quantified fluorometrically, followed by equilibration of target-enrichment pools to equimolar concentrations. All reactions were pooled and the size range of the pool selected for an average length of 450 bp using the BluePippin size select system (Sage Science, Beverly, MA, USA). Following quantification using qPCR, 300bp paired-end reads were generated on the Illumina MiSeq platform with the v3 reagent kit (Illumina, San Diego, CA, USA).

Sequence assembly, alignment, and gene tree discordance

Sequencing adapters were removed from paired-end sequencing reads using Trimmomatic (version 0.22; Bolger, Lohse & Usadel, 2014). Reads were trimmed further using a sliding window of size four, with an average quality of 15 or greater required for bases within the window to be retained. Quality trimmed reads shorter than 75 bp were discarded. The HybPiper pipeline (version 1.2; Johnson et al., 2016) was used to identify target sequences from enriched sequencing libraries by comparing all quality trimmed reads against the cnidarian reference protein collection from Kayal et al. (2018) using translated BLAST queries. After extraction of target reads and binning by locus, reads were assembled for each species using the SPAdes assembler (version 3.10.1; Bankevich et al., 2012), followed by alignment and scaffolding of contigs against the reference proteins using Exonerate (version 2.2.0; Slater & Birney, 2005). In-frame coding sequences of scaffolds (super-contigs) were translated into amino acids. All loci were individually aligned using MAFFT (version 7.271; Katoh & Standley, 2013) and ambiguous alignment positions were removed using Gblocks (version 0.91b; Talavera & Castresana, 2007) via Gblockwrapper (version 0.03; https://goo.gl/fDjan6). Gene trees and their support from 1,000 non-parametric bootstrap replicates were inferred using RAxML (version 8.2.12; Stamatakis, 2006) under the best fit model of sequence evolution (LG, WAG, or MtZoa; see phylogenetic inference below) using the computational resources provided by the Open Science Grid (Pordes et al., 2007; Sfiligoi et al., 2009). Gene trees were summarized by maximizing the number of shared quartet trees using ASTRAL-III (version 5.7.7; Zhang et al., 2018). Using the resulting phylogenetic tree and quartet frequenceies enabled us to identify areas of gene tree discordance. Following alignment and removal of ambiguous alignment positions, all loci were concatenated for further phylogenetic analyses.

Phylogenetic inference

Phylobayes (version 4.1c; Lartillot et al., 2013) was used to run eight independent MCMC chains and the posterior probability distribution was sampled until chains converged and a large sample of trees was generated. To account for site-specific differences in the evolutionary rates within and among alignment partitions, site-specific rates were inferred from the data using the CAT-GTR model during MCMC runs (Lartillot & Philippe, 2004). PartitionFinder (version 2.1.1; Lanfear et al., 2017) was used to determine the best partitioning scheme for the concatenated alignment, with maximum likelihood trees for the partitioning analysis inferred using RAxML (version 8.2.12; Stamatakis, 2006). In the absence of prior information on possible partitioning schemes, the relaxed clustering algorithm (rclusterf; Lanfear et al., 2014) was used to identify partitioning schemes that fit the data well. Initial clustering analyses included all substitution models implemented in PartitionFinder. These exploratory analyses failed to finish after more than a month of run-time but indicated that models with rate heterogeneity modeled by drawing from the gamma distribution (+G) and amino acid residue equilibrium frequencies estimated from the data (+F) fit data partitions best. Due to their prevalence in preliminary results and tractability of partitioning analysis, the final partitioning scheme was inferred using the LG (Le & Gascuel, 2008), WAG (Whelan & Goldman, 2001), and MtZoa (Rota-Stabelli, Yang & Telford, 2009) substitution matrices. The resulting partitioning scheme was used to infer the maximum likelihood phylogeny using RAxML (Stamatakis, 2006). The best tree was chosen from a set of 10 trees inferred from independent searches, starting from different random starting trees. Robustness of the resulting maximum likelihood phylogeny was assessed using 681 non-parametric bootstrap replicates.

Tree topology hypothesis testing

Minimally constrained maximum likelihood phylogenetic tree searches were conducted to evaluate which lineage of Hydroidolina is the sister group to the remainder of Hydroidolina: Aplanulata, Capitata (Cartwright & Nawrocki, 2010), or Siphonophorae (Kayal et al., 2015). In addition, tree topologies were fully constrained following the phylogenetic hypotheses proposed in this contribution and previous phylogenetic hypotheses of Hydroidolina (Cartwright & Nawrocki, 2010; Kayal et al., 2015). To discriminate between competing hypotheses of hydroidolinan relationships inferred under these comprehensive constraints, phylogenetic analyses were conducted using maximum likelihood followed by topological hypothesis tests. In short, the backbone of each tree, including the best trees found in unconstrained searches presented in this contribution, was fixed at nodes that define higher-level relationships of in- and outgroups. Under these constraints, tree inferences were able to rearrange the topology of unconstrained nodes and make adjustments to branch lengths to maximize the likelihood of the tree given the data in the concatenated amino acid residue alignment. These constrained maximum likelihood phylogenies were reconstructed using 10 independent partitioned RAxML searches, retaining the tree with the highest log likelihood (lnL). Further, the two competing hypotheses presented in this contribution were compared using fully constrained Bayesian inferences. By using fully constrained starting phylogenies, we set informative priors on tree searches sensu Bergsten, Nilsson & Ronquist (2013). Samples of the Bayesian posterior probability distribution were obtained from four independent MCMC chains under both the CAT-GTR and default CAT-Poisson (F81) models implemented in PhyloBayes, as described under phylogenetic inference above. Sampling the posterior of tree searches under two different models allowed us to identify the sensitivity of the phylogenetic inference to model misspecifications.

The Bayesian inference presented herein was used as the null hypothesis (T₀) against which all alternate trees were evaluated quantitatively using the likelihood ratio. The likelihood ratio statistic was calculated as ${\displaystyle \delta lnL=2\left(lnL{T}_A-lnL{T}_0\right)}$ where lnL T₀ represents the likelihood of the tree under the null hypothesis and lnL T_A the likelihood of the tree under the alternate hypothesis. Here, a positive δlnL indicates a better fit of T₀ to the data while a negative δlnL indicates that T_A explains the data better. Resampling of site log likelihoods (RELL; Kishino, Miyata & Hasegawa, 1990; Hasegawa & Kishino, 1994) was used to generate 10,000 bootstrap samples for each tree hypothesis, estimating the variance of tree likelihoods. Using these RELL bootstrap distributions, the null hypothesis that trees of the candidate set have the same lnL was tested via the approximately unbiased test (AU; Shimodaira, 2002), as implemented in Consel (version 0.2; Shimodaira & Hasegawa, 2001). Phylogenetic tree fit to the data was further evaluated using the posterior probability of all candidate trees using the Bayesian Information Criterion (BIC; Schwarz, 1978) approximation implemented in Consel.

Bayes factors (reviewed in Morey, Romeijn & Rouder, 2016) were calculated to compare the Bayesian majority rule consesus topology (T₀) with the alternate maximum likelihood topology (T₁) inferred in this constribution under the CAT-GTR and CAT-Poisson (F81) models in pairwise comparisons. These comparisons address the question of how well the best tree hypothesis predicts the observed amino acid alignment compared to the alternate tree hypotheses under a given substitution model. Bayesian model evaluation requires quantifying model evidence by a marginal likelihood function through integration of the product of the likelihood and the prior (Fourment et al., 2020). Model comparisons with Bayes factors use ratios of marginalized likelihoods, which are difficult to compute exactly. To address this issue, the harmonic mean of the posterior is widely used as an estimator of the marginal likelihood (Kass & Raftery, 1995). Combining the posterior distributions of four MCMC chains for a given constrained tree search, we calculated the moving harmonic mean of the posterior using a sliding window with a step length of one and a size equal to the current MCMC cycle times 0.33. The natural logarithm of the ratio of harmonic means was used for tree comparisons and interpreted following a modified version of Jeffreys’ (Jeffreys, 1961) categories of evidence (cf. Kass & Raftery, 1995). In particular, Bayes factors were calculated as ${\displaystyle B{F}_{10}=2\ast ln\left(lnL{T}_1-lnL{T}_0\right)}$ where lnL T₀ is the marginal likelihood of the MCMC search under the null hypothesis and T₁ the marginal likelihood of the alternate MCMC tree search. Note that the likelihood ratio may be negative leading to an undefined result of the logarithmic function, an issue we took into account by taking the absolute of the ratio if negative, followed by multiplying the resulting Bayes factor by negative one.

Recovery of target loci and alignment

The final concatenated alignment (Fig. 1A; Supplemental File 2) contained 44 medusozoan cnidarians, sequence data for 18 of which were generated in this study (Table 2). This alignment contained 134 of the 355 targeted loci, with a total alignment length of 21,816 character columns. Loci targeted but not included in the final alignment failed to generate sequence data from target enrichment reactions and were excluded from further analysis. Alignment lengths per locus ranged from 56 to 363 positions (average 164), with a combined total of 333 to 5,818 (average 2,891) amino acids contained in each alignment (Fig. 1B). Taxon occupancy per locus ranged from 5 to 36 species with an average of 21 species included in each alignment partition for phylogenetic analysis (Fig. 1C). Among ingroup hydrozoans, coverage varied between and within clades, with matrix completeness being highest for taxa of the Kayal et al. (2018) reference dataset (Fig. 1A). In particular, sequence data for Aplanulata, Siphonophorae, and Filifera III largely represent previously published sequence data that were used for bait design (Table 2). Consequently, alignment coverage in these groups ranges from some 40% to almost 80%. Despite being closely related to at least some of the taxa used in bait design, success in recovering loci from Podocoryne martinicana (Filifera III) and Athorybia rosacea (Siphonophorae) was limited (Table 2). In species of Filifera I, Filifera II, and Capitata target capture yields ranged from as few as some 500 amino acid residues spread over five loci (Myrionema hargitti) to as high as some 4,400 residues contained in 46 separate loci (Pennaria disticha). By contrast, the RNA bait set and hybridzation protocol employed here was highly successful in recovering sequence data for leptothecate hydrozoans despite being somewhat distantly related to any of the reference taxa. On average, 49 target loci comprising close to 5,000 amino acid residues were recovered for leptothecate species (Table 2). Here, the maximum yield was some 8,000 residues spread across 80 separate loci (Kirchenpaueria sp.). Consistent with the high degree of missing sequence data, gene trees displayed a large amount of discordance for both ingroup and outgroup taxa (Fig. 2).

Phylogeny of Hydroidolina inferred from target-capture sequencing

Bayesian MCMC chains using the CAT-GTR model were terminated after some 14,000 iterations, of which the first 3,000 were discarded as burn-in. Chains were thinned by sampling every tenth generation, yielding a mean difference between chains of 9e−4 and a maximum difference of 0.0469. The majority rule consensus of the posterior shows a well-resolved phylogeny with high posterior probabilities for the monophyly of the major ingroup and outgroup clades (Fig. 3). Aplanulata represents the sister group to the remainder of Hydroidolina. As in previous analyses, Filifera represents a polyphyletic taxon. Interestingly, Filifera I plus II are each other’s closest relative as are Filifera III plus IV. Filifera I plus II are the closest relatives to Capitata, albeit with a posterior probability slightly less than 0.95; Leptothecata is inferred as sister to the clade of Filifera I plus II and Capiata. Filifera III plus IV are the closest relatives of Siphonophorae.

For ML phylogenetic inference, the best partitioning scheme of the 134 locus dataset contained 25 partitions. In particular, the ML phylogeny was inferred using 22 partitions comprising 128 loci under the LG+G+F model, two partitions with a total of 5 loci under the MTZOA+G+F model, and 1 partition with a single locus under the WAG+G+F model. The ML tree largely agrees with the Bayesian inference (Fig. 3), but differs in the placement of Filifera I plus II relative to Capitata and Leptothecata. In contrast to the Bayesian inference, Filifera I plus II are the sister to Leptothecata rather than Capitata. However, this result lacks support with a bootstrap of less than 50. While the majority of the other nodes in both the in- and outgroup received high bootstrap support, the deep ingroup node uniting Filifera III plus IV and Siphonophorae is weakly supported with a bootstrap of 65. Despite these differences between Bayesian and ML phylogenies, both analyses agree in confidently inferring Aplanulata to be the sister lineage of the rest of Hydroidolina.

Phylogenetic hypothesis testing

Aplanulata represents the most likely sister taxon to the remainder of Hydroidolina (Fig. 4 left) with the alternate hypotheses of Capitata (Fig. 4 center) or Siphonophorae (Fig. 4 right) both receiving low BIC posterior probabilities and being rejected by the AU test. Fully constrained phylogenetic inferences provided further insight into the relationships within Hydroidolina (Fig. 5A). Two phylogenies, T₂ (cf. Cartwright & Nawrocki, 2010) and T₃ (cf. Kayal et al., 2015), were rejected by the AU test (Table 3). Interestingly, T₂ had a higher likelihood than T₃, indicating a worse fit of T₃ to the data compared to T₂. T₀ and T₁ were constrained following the Bayesian and Maximum likelihood phylogenies inferred herein (Fig. 3). Both tree hypotheses are roughly equally likely, with the maximum likelihood phylogeny T₁ being a slightly better fit to the data using maximum likelihood tree searches (Table 3). Constrained Bayesian searches using T₀ and T₁ as priors generated a sample of 27,000 log likelihoods from the posterior, with chains mixing well after burn-in (Fig. 5B). The marginal likelihood of T₀ was higher than the marginal likelihood of T₁ under the GTR model while the marginal likelihood of T₁ was higher under the F81 model (Fig. 5C). Under the GTR substitution model, T₀ is better at predicting the data than T₁ (BF₁₀ = 10.48; Fig. 5D). By contrast, T₁ is favored under the F81 model (BF₁₀ =  − 8.84; Fig. 5D).

A new phylogenetic hypothesis for Hydroidolina

Despite highly variable recovery rates of target loci that leave large gaps in the final alignment (Fig. 1A), the dataset analyzed here provides high resolution of the deep phylogeny of Hydroidolina (Fig. 3) even though discordance among gene trees is apparent (Fig. 2). Considering the ancient radiation of medusozoan cnidarians, the gene tree discordance observed here is unlikely a result of incomplete lineage sorting but rather the result of gene tree estimation error driven by the limited information contained in individual amino acid alignments. Coalescent-based summary methods lack accuracy in the presence of substantial gene tree estimation error (Warnow, 2015) and, under these conditions, phylogenetic analyses of fully partitioned concatenated alignments are preferable over summary methods that implement multi-species coalescent models (cf. Molloy & Warnow, 2018; Bryant & Hahn, 2020). While missing data could affect tree topology inference from concatenated alignments, in practice accurately placing taxa despite missing information is often not a major concern (Wiens & Morrill, 2011). Indeed, we were able to infer a well-resolved phylogeny by combining a publicly available data-rich amino acid sequence dataset with new data from key hydrozoan taxa that have so far been absent from multi-gene phylogenomic datasets. Overall topologies recovered in our phylogenetic inferences are consistent with previous phylogenomic analyses (Zapata et al., 2015; Kayal et al., 2018). While neither Zapata et al. (2015) nor Kayal et al. (2018) included sufficient taxa representing Hydroidolina to allow for much insight into the evolutionary history of the group, both recovered Aplanulata as the sister clade of the remaining hydroidolinan groups, a placement we confirmed after comprehensively sampling higher-level taxa of Hydroidolina (Fig. 3). This placement is at odds with previous phylogenetic treatments of Hydroidolina that inferred a more recent origin of Aplanulata within Hydroidolina as sister of Filifera I plus II (Cartwright & Nawrocki, 2010; Kayal et al., 2015). Placing Aplanulata as the sister-group of Hydroidolina has potentially far-reaching consequences for our understanding of hydrozoan evolution. Aplanulata is dominated by solitary polyp-forming taxa, a trait common in other medusozoans but rare across Hydroidolina in which colonial hydroids, many of which are polymorphic, displaying a reproductive division of labor, dominate (Cartwright & Nawrocki, 2010; Cartwright, Travert & Sanders, 2020). As such, the placement of Aplanulata as sister to the remainder of Hydroidolina suggests that the last common ancestor of Hydroidolina possessed a solitary polyp that produces medusae as part of a metagenetic life cycle in which an asexually reproducing polyp generation alternates with a sexually reproducing medusa generation (cf. Bentlage et al., 2018). Pending formal analyses using ancestral state reconstructions with greater taxon sampling, the evolution of coloniality and division of labor, including losses of metagenesis (lack of free-swimming medusae), are likely derived traits in Hydroidolina.

Ingroup relationships within Hydroidolina appear largely congruent across the Bayesian posterior tree set and the maximum likelihood phylogeny with the exception of the placement of Filifera I plus II with respect to Capitata and Leptothecata (Fig. 3). One possible explanation for this disagreement may be that the substitution matrix-based models used for maximum likelihood inference here are potentially less accurate in reflecting the substitution process that led to the observed data. Exploratory partitioning analyses with PartitionFinder suggested that substitutions in a significant number of partitions should be inferred using the GTR model. However, benchmarking of partitioning analyses under the GTR model indicated that computational time would be prohibitive to search across this large parameter space. To address this issue, partitioning schemes were inferred using simpler substitution matrices (i.e., LG, WAG, and MtZoa). Regardless of this poorly supported topological difference in the placement of Filifera I plus II (Fig. 3), the overall relationships recovered are consistent across analyses but at odds with previous phylogenetic hypotheses of Hydroidolina (cf. Cartwright & Nawrocki, 2010; Kayal et al., 2015). Despite conflicts between the phylogenetic hypotheses presented here and previous treatments of Hydroidolina, a consistent result that has emerged over the last decade is the polyphyly of Filifera. In agreement with prior phylogenetic hypotheses (Cartwright & Nawrocki, 2010; Kayal et al., 2015), we find that Filifera I plus II forms a monophyletic clade, as does Filifera III plus IV. As of now, these clades remain without clear definition based on morphological or life history characteristics but Cartwright et al. (2008) provide the most detailed discussion of possible synapomorphies of hydroidolinan clades to date.

Phylogenetic tree selection

Our results appear unequivocal on the position of Aplanulata as sister to the remaining lineages of Hydroidolina (Figs. 2–4). To facilitate discriminating between competing hypotheses of hydroidolinan relationships, we quantified the evidence for the set of alternate phylogenetic hypotheses (Fig. 5A) in light of the multi-locus dataset analyzed here. Likelihood ratio statistics (Table 3) indicate that the phylogenies proposed in this contribution (Fig. 3; Fig. 5A T₀ and T₁) are a better fit than either the previously published rDNA (Fig. 5A T₂; Cartwright & Nawrocki, 2010) or mitochondrial genome-based (Fig. 5A T₃; Kayal et al., 2015) phylogenetic hypotheses. Indeed, T₂ and T₃ were rejected by the AU test while neither T₀ nor T₁ could be rejected (Table 3). In addition, the majority rule consensus of the Bayesian posterior (T₀) and the maximum likelihood topology (T₁) received high BIC-based posterior probabilities compared to the very small posterior probabilities of the alternate hypotheses (Table 3).

The hypothesis tests above reduced the set of credible trees by eliminating unlikely tree topologies. The procedures employed rely on the likelihood of a single phylogeny, implicitly assuming that the tree is known without error. Resampling of site likelihoods from this best tree was used to estimate variances of tree likelihoods for comparisons. Bayesian approaches that estimate the posterior probability distribution of phylogenetic hypotheses allow averaging across tree topologies and branch lengths for more comprehensive incorporation of uncertainty into comparisons of phylogenetic trees but may be affected adversely by misspecifications of the prior (e.g., Bergsten, Nilsson & Ronquist, 2013). Bayes factors incorporate phylogenetic uncertainty in model comparisons by employing ratios of marginal likelihoods and allowed us to quantify the weight of evidence for a topological hypothesis under a given amino acid substitution model. We found that T₀ predicts the alignment data best under the GTR model while T₁ predicts the alignment best under the simpler F81 model (Figs. 5C & 5D). The difference between these two topologies is limited to the placement of Filifera I plus II with respect to Capitata and Leptothecata (Fig. 3) and we suspect that the sensitivity to model choice may be driven by the lack of data for Filifera I plus II (Fig. 1) and the resulting gene tree discordance (Fig. 2).

Despite being unable to confidently choose a single phylogenetic hypothesis as fitting the data best using maximum likelihood and Bayesian hypothesis testing frameworks, the phylogenetic hypotheses presented in this contribution (Fig. 3) consistently display high predictive power for the alignment data compared to the alternate hypotheses proposed previously. Consequently, we suggest that our phylogenetic framework for Hydroidolina represents the most viable hypothesis of hydroidolinan relationships to date. That said, further studies with additional taxa and characters are still needed to assess the validity of this working hypothesis and clarify the affinities of Filifera I plus II in particular. Such analyses will also allow for more thorough analyses of the evolution of life history characteristics across this most speciose of medusozoan clades.