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

Higher-level relationships of the Hydrozoan subclass Hydroidolina, which encompasses the vast majority of medusozoan cnidarian species diversity, have been elusive to confidently infer. The most widely adopted phylogenetic framework for Hydroidolina based on ribosomal RNA data received low support for several higher level relationships. To address this issue, we developed a set of RNA baits to target more than a hundred loci from the genomes of a broad taxonomic sample of Hydroidolina for high-throughput sequencing. Using these data, we inferred the relationships of Hydroidolina using maximum likelihood and Bayesian approaches. Both inference methods yielded well-supported phylogenetic hypotheses that largely agree with each other. Using maximum likelihood and Baysian hypothesis testing frameworks, we found that several alternate topological hypotheses proposed previously may be rejected in light of the genomic data generated for this study. Both the maximum likelihood and Bayesian topologies inferred herein consistently score well across testing frameworks, suggesting that their consensus represents the most likely phylogenetic hypothesis of Hydroidolina. This phylogenetic framework places Aplanulata as sister lineage to the remainder of Hydroidolina. This is a strong deviation from previous phylogenetic analyses that placed Capitata or Siphonophorae as sister group to the remainder of Hydroidolina. Considering that Aplanulata represents a lineage comprised of species that for the most part possess a life cycle involving a solitary polyp and free-swimming medusa stage, the phylogenetic hypotheses presented herein have potentially large implications for clarifying the evolution of life cycles, coloniality, and the division of labor in Hydrozoa as taxon sampling for phylogenetic analyses becomes more complete.


INTRODUCTION
While the fossil record of medusozoan cnidarians is scant, the origin of the group has been inferred to lie near the end of the Ediacaran, approximately 550-580 million years ago (Han, Zang & Komiya, 2016). Plausible crown-group hydrozoans have been described from some 500 million year old Cambrian deposits (Cartwright et al., 2007), suggesting an ancient origin of extant hydrozoans likely dating back to the period of rapid diversification of metazoan life during which all major modern animal phyla emerged (Valentine, Jablonski & Erwin, 1999;Erwin, 2020). Hydrozoans are of particular interest in the study of the evolution of development, as their radiation gave rise to diverse life cycle strategies, diverse forms of coloniality and the division of labor (Cartwright & Nawrocki, 2010;Bentlage et al., 2018;Cartwright, Travert & Sanders, 2020). This diversity is concentrated in the hydrozoan subclass Hydroidolina, the medusozoan clade that contains the vast majority of the 3,800 nominal hydrozoan species (Daly et al., 2007;Schuchert, 2020). Elucidating the evolutionary history and patterns of complex character evolution in Hydroidolina requires a solid understanding of the phylogenetic history of the group (e.g., Cartwright, Travert & Sanders, 2020).
However, the goal of inferring the deep phylogeny of Hydroidolina has been elusive, possibly as a result of the early origin and likely rapid diversification of the group. The most comprehensive phylogenetic hypothesis (Cartwright & Nawrocki, 2010) of higher-level relationships within Hydroidolina was inferred using ribosomal DNA (rDNA) from a broad taxonomic sample. While shallow nodes were well resolved with high confidence, higher-level relationships generally received weak support. In particular, Leptothecata, Siphonophorae, Capitata, and Aplanulata were inferred to be monophyletic groups while Filifera was polyphyletic (Cartwright et al., 2008;Cartwright & Nawrocki, 2010). Previously, the latter three taxa were united under the name Anthoathecata but their non-monophyly had been demonstrated earlier (Collins et al., 2006). Whole mitochondrial genomes have been employed previously to address the issue of reconstructing deep nodes within the phylogeny of Hydroidolina (Kayal et al., 2015). While this approach led to a well-resolved and highly supported phylogenetic hypothesis, several nodes of the resulting tree topology are at odds with rDNA-based phylogenies (Cartwright et al., 2008;Cartwright & Nawrocki, 2010) and recent phylogenomics-based hypotheses (Kayal et al., 2018).
Advances in understanding of medusozoan, and more broadly cnidarian relationships, were made by employing phylogenomic datasets derived from whole-transcriptome and genome-sequencing efforts (Zapata et al., 2015;Kayal et al., 2018). While these efforts provided answers to several long-standing questions of cnidarian evolutionary history, taxon sampling was insufficient for rigorous evaluation of hydrozoan relationships. We used the coding sequences generated by these phylogenomic studies as a backbone for targeted high-throughput sequencing, producing a multi-locus dataset to infer the phylogeny of Hydroidolina. For this purpose, we developed a set of custom baits to enrich target loci from a representative sample of hydroidolinan hydrozoans. This work provides a new framework for the phylogeny of Hydroidolina that will enable further phylogenetic comparative studies of character evolution in Hydrozoa. To evaluate multiple competing topological hypotheses of hydroidolinan relationships, both likelihood and Bayesian statistical frameworks (cf. Sober, 2008) were employed to discriminate between alternate tree topologies and evaluate the strength of evidence supporting prior phylogenetic hypotheses of Hydroidolina as well as those phylogenies inferred from the multi-locus dataset presented here.

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. 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  Putnam et al. (2007) 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   , 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 0 ) against which all alternate trees were evaluated quantitatively using the likelihood ratio. The likelihood ratio statistic was calculated as where lnL T 0 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 0 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 0 ) with the alternate maximum likelihood topology (T 1 ) 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 where lnL T 0 is the marginal likelihood of the MCMC search under the null hypothesis and T 1 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 2 (cf. Cartwright & Nawrocki, 2010) and T 3 (cf. Kayal et al., 2015), were rejected by the AU test (Table 3). Interestingly, T 2 had a higher likelihood than T 3 , indicating a worse fit of T 3 to the data compared to T 2 . T 0 and T 1 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 1 being a slightly better fit to the data using maximum likelihood tree searches (Table 3). Constrained Bayesian searches using T 0 and T 1 as priors generated a sample of 27,000 log likelihoods from the posterior, with chains mixing well after burn-in  5B). The marginal likelihood of T 0 was higher than the marginal likelihood of T 1 under the GTR model while the marginal likelihood of T 1 was higher under the F81 model (Fig. 5C). Under the GTR substitution model, T 0 is better at predicting the data than T 1 (BF 10 = 10.48; Fig. 5D). By contrast, T 1 is favored under the F81 model (BF 10 = −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 Figure 5 Phylogenetic hypothesis tests using fully constrained tree searches. (A) Phylograms inferred using fully constrained maximum likelihood searches to compare explicit hypotheses of tree topologies; constraints indicated by solid black circles; scale bar indicates substitutions per site. Tree topologies were constrained to reflect the Bayesian inference (T 0 ) and maximum likelihood phylogenies (T 1 ) presented in this contribution (Fig. 3). Two alternate tree searches were constrained to reflect previously published phylogenies (T 2 : Cartwright & Nawrocki, 2010; T 3 : Kayal et al., 2015); *** indicates a p-value of less than 0.001 for the approximately unbiased (AU) test (Table 3). Shallow nodes were collapsed and species labels removed for legibility. 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 multigene 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) Kayal et al., 2015) phylogenetic hypotheses. Indeed, T 2 and T 3 were rejected by the AU test while neither T 0 nor T 1 could be rejected (Table 3). In addition, the majority rule consensus of the Bayesian posterior (T 0 ) and the maximum likelihood topology (T 1 ) 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 0 predicts the alignment data best under the GTR model while T 1 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.

CONCLUSIONS
This study represents one of only a few attempts at using target-enriched high throughput sequencing to generate a multi-locus alignment for phylogenetic analyses of Cnidaria. Previous attempts had been limited to Anthozoa (Quattrini et al., 2018;Erickson et al., 2021), relying on the greater availability of genomic data for Anthozoa compared to Medusozoa. Despite mixed success in capturing target loci, we were able to generate an informative multi-gene alignment that produced a well-supported phylogenetic hypothesis for Hydroidolina. One drawback of bait-development based on coding sequences alone, as done here, is that baits may inadvertently be designed across intron-exon boundaries which is likely to reduce bait effectiveness. As availability of genomic resources increases, it will be possible to address these issues and apply target capture sequencing approaches routinely to phylogenetic studies of medusozoan cnidarians, building on the growing knowledge-base for designing targeted high-throughput sequencing experiments (reviewed in Andermann et al., 2020). Nonetheless, we were able to address long-standing questions in hydrozoan phylogenetics by applying relevant constraints for discriminating between alternate phylogenetic hypotheses for Hydroidolina. We suggest that the consensus of T 0 and T 1 (Fig. 5A) represents the most likely tree topology for hydroidolinan evolutionary relationships to date. While the likelihood-based AU test (Table 3) allowed rejection of several alternate topologies, Bayes factors ( Fig. 5D) suggested that the choice of substitution model affects the posterior distributions of the two alternate topologies inferred. As such, the most viable representation of phylogenetic relationships of Hydroidolina is the consensus of the two topologies inferred herein.