Dataset: Siphonophore Phylogeny

Statistical analyses
We used a series of phylogenetic comparative methods to test the evolutionary hypotheses presented in this study. We reconstructed ancestral states using ML (R phytools::anc.ML [@revell2012phytools]), and stochastic character mapping (R phytools::make.simmap) for categorical characters. For more details on the data wrangling prior to these analyses, please see the Supplementary Methods. R scripts are available in the Dryad repository.

In order to study the evolution of predatory specialization, we reconstructed components of the diet and prey selectivity on the phylogeny using ML (R phytools::anc.ML). To identify evolutionary associations of diet with tentillum and nematocyst characters, we compared the performance of a neutral evolution model to that of a diet-driven directional selection model. First, we collapsed the diet data into the five feeding guilds mentioned above (fish specialist, small crustacean specialist, large crustacean specialist, gelatinous specialist, generalist), based on which prey types they were observed consuming most frequently. Then, we reconstructed the feeding guild ancestral states using the ML function ace (package ape [@paradis2019package]), removing tips with no feeding data. The ML reconstruction was congruent with the consensus stochastic character mapping (S15). Then, using the package *OUwie* [@beaulieu2012ouwie], we fitted an OU model with multiple optima and rates of evolution (OUm) matched to the reconstructed ancestral diet regimes, a single optimum OU model, and a BM null model, inspired by the analyses in [@cressler2015detecting]. We then ranked the models in order of increasing parametric complexity (BM, OU, OUm), and compared the corrected Akaike Information Criterion (AICc) support scores [@sugiura1978further] to the lowest (best) score, using a cutoff of 2 units to determine significantly better support. When the best fitting model was not significantly better than a less complex alternative, we selected the least complex model (S9). In addition, we calculated and reported the model adequacy scores using the R package *arbutus* [@pennell2015model].

In order to study correlations between the rates of evolution between different characters, we fitted a set of evolutionary variance-covariance matrices [@revell2009phylogenetic] (R phytools::evol.vcv). For more details on the data wrangling preceding these analyses, please see Supplementary Methods. To test whether phenotypic integration changed across selective regimes determined by the reconstructed feeding guilds, we carried out character-pairwise variance-covariance analysis comparing alternative models (R phytools::evolvcv.lite), including those where correlations are the same across the whole tree and models where correlations differ between selective regimes (S19). Number of taxa used in each pairwise comparison is reported in S20. Finally, we compared regime-specific variance-covariance matrices to the general matrix and to their preceding regime matrix to identify the changes in character dependences unique to each regime (S21-22). 

We carried out a linear discriminant analysis of principal components (DAPC) using the dapc function (R adegenet::dapc) [@jombart2010discriminant]. This function allowed us to incorporate more predictors than individuals. We generated discriminant functions for feeding guild, and for the presence of copepods, fish, and shrimp (large crustaceans) in the diet (S10-13). From these DAPCs we obtained the highest contributing morphological characters to the discrimination (characters in the top quartile of the weighted sum of the linear discriminant loadings controlling for the eigenvalue of each discriminant). In order to identify the sign of the relationship between the predictor characters and prey type presence in the diet, we then generated generalized logistic regression models (as a type of generalized linear model, or GLM using R stats::glm) and phylogenetic generalized linear models (R phylolm::phyloglm) with the top contributing characters (from the corresponding DAPC) as predictors (S14). We also carried out these GLMs on the Ivlev’s selectivity indices for each prey type calculated from [@purcell1981dietary]. In addition, we ran a series of comparative analyses to address hypotheses of diet-tentillum relationships posed in the literature. Additional details on the DAPC optimization are available in the Supplementary Methods.

Supplementary Methods
Phylogenetic inference:
We aligned the sequences using MAFFT [@katoh2002mafft] (alignments available in Dryad). We inferred a Maximum Likelihood (ML) phylogeny (S2) from 16S and 18S ribosomal rRNA genes using IQTree [@nguyen2014iq] with 1000 bootstrap replicates (iqtree -s alignment.fa -nt AUTO -bb 1000). We used ModelFinder [@kalyaanamoorthy2017modelfinder] implemented in IQTree v1.5.5. to assess the relative model fit. ModelFinder selected GTR+R4 for having the lowest Bayesian Information Criterion score. Additionally, we inferred a Bayesian tree with each gene as an independent partition in RevBayes [@hohna2016revbayes] (S5), which was topologically congruent with the unconstrained ML tree. The *alpha* priors were selected to minimize prior load in site variation.

Data wrangling for comparative analyses:
For comparative analyses, we removed species present in the tree but not represented in the morphology data, and *vice versa*. Although we measured specimens labeled as *Nanomia bijuga* and *Nanomia cara*, we are not confident in some of the species-level identifications, and some specimens were missing diagnostic zooids. Thus, we decided to collapse these into a single taxonomic concept (*Nanomia* sp.). All *Nanomia* sp. observations were matched to the phylogenetic position of *Nanomia bijuga* in the tree. We carried out all phylogenetic comparative statistical analyses in the programming environment R [@team2017r], using the Bayesian ultrametric species tree (S8), and incorporating intraspecific variation estimated from the specimen data as standard error whenever the analysis tool allowed it. R scripts and summarized species-collapsed data available in the Dryad repository. For each character (or character pair) analyzed, we removed species with missing data and reported the number of taxa included. We tested each character for normality using the Shapiro-Wilk test [@shapiro1965analysis] and log-transformed those that were non-normal.

Data wrangling for the variance-covariance analyses:
When fitting all variance-covariance terms simultaneously (S16-18), we selected the largest set of characters that would allow the analysis to run without computational singularities. This excluded many of the morphometric characters which are linearly dependent on other characters. Since the functions do not tolerate missing data, we ran the analyses in two ways: One including all taxa but transforming absent states to zeroes, and another removing the taxa with absent states. These analyses could only be carried out on the subset of taxa for which diet data is available, and only among character pairs that are not computationally singular for that taxonomic subset. Gelatinous specialist correlations could only be estimated for a small subset of characters present in *Apolemia* (S21F, S22E, S23D) and should be interpreted with care.

Comparative tools used to test character associations:
To test for correlated evolution among binary characters, we used Pagel’s test [@pagel1994detecting]. To characterize and evaluate the relationship between continuous characters, we used phylogenetic generalized least squares regressions (PGLS) [@grafen1989phylogenetic]. To compare the evolution of continuous characters with categorical aspects of the diet, we carried out a phylogenetic logistic regression (R nlme::gls using the 'corBrownian' function for the argument ‘correlation’).

DAPC optimization:
Some taxa have inapplicable states for certain absent characters (such as the length of a nematocyst subtype that is not present in a species), which are problematic for DAPC analyses. We tackled this by transforming the absent states to zeroes. This approach allows us to incorporate all the data, but creates an attraction bias between small character states (*e.g.* small tentilla) and absent states (*e.g.* no tentilla). Absent characters are likely to be very biologically relevant to prey capture and we believe they should be accounted for in a predictive approach. We limited the number of linear discriminant functions retained to the number of groupings in each case. We selected the number of principal components retained using the a-score optimization function (R adegenet::optim.a.score) [@jombart2010discriminant] with 100 iterations, which yielded more stable results than the cross validation function (R adegenet::xval). This optimization aims to find the compromise value with highest discrimination power with the least overfitting. The discriminant analysis for feeding guild (7 principal components, 4 discriminants) produced 100% discrimination, and the highest loading contributions were found for the characters (ordered from highest to lowest): Involucrum length, heteroneme volume, heteroneme number, total heteroneme volume, tentacle width, heteroneme length, total nematocyst volume, and heteroneme width (S10).