Evolution of Gene Expression across Species and Specialized Zooids in Siphonophora

Abstract Siphonophores are complex colonial animals, consisting of asexually produced bodies (zooids) that are functionally specialized for specific tasks, including feeding, swimming, and sexual reproduction. Though this extreme functional specialization has captivated biologists for generations, its genomic underpinnings remain unknown. We use RNA-seq to investigate gene expression patterns in five zooids and one specialized tissue across seven siphonophore species. Analyses of gene expression across species present several challenges, including identification of comparable expression changes on gene trees with complex histories of speciation, duplication, and loss. We examine gene expression within species, conduct classical analyses examining expression patterns between species, and introduce species branch filtering, which allows us to examine the evolution of expression across species in a phylogenetic framework. Within and across species, we identified hundreds of zooid-specific and species-specific genes, as well as a number of putative transcription factors showing differential expression in particular zooids and developmental stages. We found that gene expression patterns tended to be largely consistent in zooids with the same function across species, but also some large lineage-specific shifts in gene expression. Our findings show that patterns of gene expression have the potential to define zooids in colonial organisms. Traditional analyses of the evolution of gene expression focus on the tips of gene phylogenies, identifying large-scale expression patterns that are zooid or species variable. The new explicit phylogenetic approach we propose here focuses on branches (not tips) offering a deeper evolutionary perspective into specific changes in gene expression within zooids along all branches of the gene (and species) trees.


Introduction
Colonial animals, consisting of genetically identical bodies that are physically and physiologically connected, can be found across the metazoan tree of life (Hiebert et al. 2021). Functional specialization of such bodies has evolved multiple times within colonial animals, with siphonophores in particular showing the greatest diversity of functionally specialized bodies (Beklemishev 1969). Siphonophores are highly complex, colonial "superorganisms" consisting of asexually produced bodies (termed zooids) that are homologous to solitary free-living polyps and medusae (the typical body forms in Cnidaria), but that share a common gastrovascular cavity (Mackie 1963(Mackie , 1986Totton 1965;Mackie et al. 1987;Dunn and Wagner 2006). The extreme specialization of siphonophore zooids has been of central interest to zoologists since the 19th century, in part because these zooids are highly interdependent (Mackie 1963;Beklemishev 1969).
Siphonophore zooids have been likened to animal organs, with each functionally specialized zooid performing distinct roles within the colony ( fig. 1) (Mackie 1963). Although this analogy works well for understanding the function of zooids within the colony as a whole, the analogy falls short in terms of explaining the evolutionary origin and development of these biological units: functionally specialized zooids are evolutionarily homologous to free living organisms, they are multicellular, possess distinct zooid-specific cell types, and show regional subfunctionalization (Mackie 1960;Totton 1965;Church et al. 2015). Although the developmental mechanisms generating zooids are very different in different clades of colonial animals (Carr e 1967(Carr e , 1969; Carr e and Carr e 1991; Dunn and Wagner 2006;Siebert et al. 2015), the evolutionary processes acting on zooids may be similar to those acting on other modular biological units such as cell type, tissue, and organ (Hiebert et al. 2021). The cellular and molecular processes underlying the patterning and molecular function of functionally specialized cnidarian zooids remain an open biological question. In recent years, there has been a focus in particular on differential gene expression (DGE) patterns found in different functionally specialized zooids (Siebert et al. 2011;Plachetzki et al. 2014;Sanders et al. 2014. Efforts to investigate the functional specialization of siphonophores have been limited, in part because there have been few detailed investigations of zooid structure. In the last half century, the microanatomy of siphonophore zooids and tissues has been investigated in only a handful of siphonophore species (Mackie 1960; Carr e 1969; Bardi and Marques 2007;Church et al. 2015). This leaves many unknowns about how zooid structure and function differ across zooid types and species. Recent in situ gene expression analyses in siphonophores have described where a small number of preselected genes are expressed at high spatial resolution (Siebert et al. 2011Church et al. 2015), but these methods are limited since they require a large number of specimens per gene and siphonophores are relatively difficult to collect. RNA-seq analyses of hand-dissected specimens (Siebert et al. 2011;Sanders et al. 2014Macrander et al. 2015), in contrast, can describe the expression of a very large number of genes at low spatial resolution. The fact that so many data are obtained from each specimen is particularly advantageous in difficult-to-collect organisms like siphonophores. An earlier RNA-seq study of gene expression in two zooid types in a single siphonophore species showed the potential of this method to better understand differences between zooids (Siebert et al. 2011).
In this study, we use RNA-seq data to investigate the functional specialization and evolution of zooids across siphonophore species. We address this question at three comparative levels, namely within species, between species, and across species incorporating an explicit phylogenetic approach. To carry out these analyses, we also address three key challenges to study the evolution of gene expression in a comparative framework (Dunn, Luo, et al. 2013). First, we introduce a new metric to normalize gene expression which is valid for comparison across species. A common metric used to quantify gene expression is transcripts per million (TPM) (Li and Dewey 2011;Wagner et al. 2012). TPM is a relative measure of expression which depends on the number of genes present in a reference. Hence, TPM values are a valid metric when comparing libraries within a single species, but are not directly comparable across species when the references are incomplete, or genes have been gained and lost in the course of evolution. To address this issue, we introduce transcripts per million 10K (TPM10K), a metric which normalizes TPM to account for different sequencing depths among species (see Materials and Methods for details). We use TPM10K in our between species and phylogenetic-based analyses.
Second, we further account for gene and species effects to compare gene expression across genes and species. Normalized read counts produced by RNA-seq experiments are proportional to gene expression, but they are also impacted by factors that can vary across genes and species such as gene sequence and length. Therefore, to estimate and compare expected gene counts it is necessary to model such factors with unknown species-and gene-specific counting-efficiency coefficients (see Dunn, Luo, et al. [2013]). A direct comparison of expected gene counts among species, however, can be misleading if differences in counts are simply due to differences in counting efficiency and not due to differences in expression. To address this issue, we use ratios of expected counts. Using ratios of counts, we are able to eliminate species-and gene-specific counting efficiency within species, as the unknown counting efficiency factor is in both the numerator and the denominator, and is thus removed prior to comparisons among species (Dunn, Luo, et al. 2013). We use ratios of expected counts in our between species and phylogenetic-based analyses.
Third, we use a novel phylogenetic approach to investigate the evolution of gene expression on gene trees with complex evolutionary histories. Most evolutionary studies of gene expression focus exclusively on comparing expression values of strict orthologs (i.e., gene lineages related only by speciation events) which are shared across all species. Analyses of strict orthologs are limited to a subset of genes that have no evidence of duplication, and we call this approach species tree filtering (STF) (Brawand et al. 2011;Yang and Wang 2013;Levin et al. 2016;Cardoso-Moreira et al. 2019). Because the history of genes is characterized by more complex scenarios involving gene duplication and speciation, we introduce a new method that takes advantage of this rich history to examine the evolution of gene expression across more genes and species. We call this method species branch filtering (SBF) ( fig. 2), as we map expression data to gene phylogenies and identify equivalent branches in the species tree which are descended from speciation events in order to make valid comparisons across species.
Using STF and SBF, combined with classical DGE analysis, we conducted DGE analyses among siphonophore zooids, and one specialized tissue (the pneumatophore, a gas-filled float) (fig. 1) in seven siphonophore species, to address three questions. First, we analyzed differences in gene expression within species to identify patterns that define particular zooid types, and gain insight into the genes that may be playing a role in the structure, maintenance, and functioning of zooids. We also examined whether gene expression patterns could distinguish novel, distinct species-specific zooid types within particular species. Second, we compared gene expression between species, using STF followed by linear models to identify zooid-and species-variable genes. And third, we used phylogenetic comparative methods and SBF to identify putative clade/lineage or zooid-specific expression patterns.

Results
Within Species Analyses: Which Genes Are Specific to Particular Zooid Types?
We sequenced mRNA from microdissected zooids from seven different siphonophore species (at least two replicate colonies each) and mapped these short-read libraries to previously published transcriptomes (Munro et al. 2018). We collected RNA-seq data from, where possible, five different zooids and one specialized tissue, the pneumatophore, as well Munro et al. . doi:10.1093/molbev/msac027 MBE as unique zooids specific to Agalma elegans (B palpons), Physalia physalis (tentacular palpon), and Bargmannia elongata (yellow and white gastrozooids), and where possible developing and mature zooids (supplementary table S1, Supplementary Material online). Due to species availability, we were not able to sample more than one replicate for some zooids-single replicates were excluded from downstream DGE analyses (see supplementary fig. S1, Supplementary Material online).
The first component of variation that we assessed was among technical replicates. The technical replicates consist of resequenced developing nectophores and developing gastrozooids from the same Frillagalma vityazi individual that were spiked in across multiple lanes and runs. Lane and run effects have been proposed as major sources of technical variability in RNA-seq data that may confound observations of biological variation (Auer and Doerge 2010;McIntyre et al. 2011). The differences between technical replicates (supplementary fig. S2, Supplementary Material online) were much smaller (0.39% variance of expression distance) than the differences between zooids (98.32% variance of expression distance). Differences among technical replicates of the same zooid were correlated with library size and run, not by lane.
The second component of variation we considered was biological variation among sampled colonies (supplementary figs. S3-S9, Supplementary Material online). Specimens were collected in the wild at different depths and over different time periods, but despite these environmental factors there was remarkably little variation among sampled colonies. Some samples did show greater variation across replicates, such as a developing palpon replicate in Nanomia bijuga The third component of variation we considered was among zooids/specialized tissues within species. This was the greatest component of variation (supplementary figs. S3-S9, Supplementary Material online). We identified genes that were significantly differentially expressed in particular zooids, and in the pneumatophore, for each of the sampled SGZ NGZ

Gastrozooid
Zooid specialized for feeding. Homologous to polyps. Ingests food, and extracellularly digests prey. Typically has a tentacle for prey capture. Has a nematogenic region that supplies nematocysts to the tentacle.

Nectophore
Zooid specialized for swimming. Homologous to medusae. Highly muscular and propels the colony via contractions that generate water jets out of the ostium

Male gonodendron
Compound structures with multiple gonophores (zooids homologous to medusae) that contain male gametes and are the site of spermatogenesis.

Female gonodendron
Compound structures with multiple gonophores (zooids homologous to medusae) that contain female gametes and are the site of oogenesis.

Pneumatophore
Tissue specialized in buoyancy, flotation, and geotaxis. Thought to generate carbon monoxide. Not a zooid.

Palpon
Zooid specialized for feeding. Derived gastrozooid. Considered to be an accessory digestive zooid. Hypothesized to play defensive and sensory functions. Typically has a reduced tentacle (palpacle).
Gene Expression across Species and Specialized Zooids in Siphonophora . doi:10.1093/molbev/msac027 MBE species, based on pairwise comparisons (supplementary file 1, Supplementary Material online). In each case, significantly differentially expressed refers to genes with higher transcript abundance relative to the other zooid. As it is possible for the same gene to be significantly differentially expressed in pairwise comparisons of several zooids, we also identified a subset of genes that were significantly differentially expressed in only one zooid but not in any other zooid within a given species (supplementary fig. S10   Step 1, we label each of the nodes in the species tree, and identify equivalent speciation nodes across every gene tree (an exemplar is shown here).
Step 2, we map expression values (TPM10K) to the tips (expression values are mapped and reconstructed for each homologous zooid separately).
Step 3, we reconstruct ancestral trait expression values at all internal nodes where expression data are available.
Step 4, we calculate scaled change in gene expression (child node expres-sionÀparent node expression/branch length). Branch length is calibrated to the species tree branch lengths.
Step 5, we identify branches in gene trees that correspond to equivalent branches in the species tree. There may be more than one branch in a gene tree that corresponds to the same branch in the species tree. Munro et al. . doi:10.1093/molbev/msac027 MBE feeding and digestion, for example, and we found these zooids to be enriched for genes involved in chitin, glutathione, and peptide catabolism, proteolysis, as well as metabolism of carbohydrates. Likewise, for male gonodendra, we found GO term enrichment for biological processes such as sperm flagellum, mitotic cell cycle process, DNA replication. For female gonodendra, there was GO term enrichment for a number of biological processes including mitotic cell cycle process, DNA replication (in A. elegans), as well as a number of signaling pathways and developmental processes (in F. vityazi).
Within the four best sampled species, we also identified 92 higher abundance putative transcription factors (

Novel Zooids within Species: Can Expression Patterns Distinguish Distinct Zooid Types?
In siphonophores, there are several instances of lineagespecific zooid diversification events. We investigated gene expression patterns between the novel zooid type and the hypothesized most closely related zooid type in three species. In B. elongata, there are two morphologically distinct gastrozooids, that we termed "white" and "yellow" gastrozooids (supplementary fig. S12A and B, Supplementary Material online). The "yellow" gastrozooid is larger and darker and occurs as the seventh to tenth gastrozooid on the stem (Dunn 2005). In the Portuguese man of war, P. physalis, the gastrozooid is unique compared with other gastrozooids in other speciesit has a mouth, but no tentacle, and the basigaster region is greatly reduced (Mackie 1960;Totton 1960). Meanwhile, the tentacle is associated with another zooid, the tentacular palpon (supplementary fig. S12C, Supplementary Material online) (Haeckel 1888;Totton 1960;Bardi and Marques 2007;Munro et al. 2019). In P. physalis, both the gastrozooid and the tentacular palpon are considered to be subfunctionalized from an ancestral gastrozooid type (Munro et al. 2019). Finally, in A. elegans, there are thought to be at least two different palpon types: gastric palpons that arise at the base of the peduncle of the gastrozooid, and a palpon called the Bpalpon (supplementary fig. S12D, Supplementary Material online) (Dunn and Wagner 2006). The distinction between these two types of palpon is based on the location of these zooids-the gastrozooid is typically the last element of each repeating pattern of zooids along the stem (cormidium), but based on the budding sequence, Dunn et al. proposed that the enlarged B-palpon is the last element in A. elegans (Dunn and Wagner 2006). Each of these cases represents a different type of novelty: in B. elongata, the distinction between zooids was made based on size and color but not on obvious differences in function, in P. physalis, the gastrozooids and tentacular palpons differ structurally and functionally, and finally in A. elegans, gastric palpons and B palpons differ only in colony location, development, and possibly size.
In P. physalis, 976 genes showed significant differential expression (in all cases, higher transcript abundance relative to the other zooid) in the mature tentacular palpon, compared with 606 genes in the mature gastrozooid (supplementary file 4, Supplementary Material online). A number of genes were significantly differentially expressed in the mature tentacular palpon relative to all other tissues, of which, 670 genes were significantly differentially expressed relative to other tissues that were not shared with the gastrozooid. In the gastrozooid, 849 genes were significantly differentially expressed relative to other tissues that were not shared with the tentacular palpon. A number of genes significantly differentially expressed in the tentacular palpon are uncharacterized, however, we identified 46 putative toxins in the tentacular palpon, including hemostasis interfering and platelet aggregation activating toxins, phospholipases, serineproteases, hydrolases, metalloendoproteases, calglandulin-like genes, and a neurotoxin (supplementary file 5, Supplementary Material online). By contrast, in the gastrozooid, we found 59 significantly differentially expressed putative toxins (supplementary file 6, Supplementary Material online), including pore-forming Conoporin-Cn1like and Tereporin-Ca1-like toxins, multiple neurotoxins, hydrolases, serine proteases, toxins likely involved in the promotion of blood coagulation and inhibition of platelet aggregation. Reflecting the role of the gastrozooid in digestion, we also find significant differential expression of digestive enzymes, including chymotrypsin-like genes.
Between the white mature gastrozooid and the yellow mature gastrozooid in B. elongata, few significantly differentially expressed genes were identified (eight genes were up in "white" mature gastrozooids relative to 36 genes up in "yellow" gastrozooids) (supplementary file 4, Supplementary Material online). Among genes that were significantly differentially expressed in either "white" or "yellow" gastrozooids relative to all other tissues, 276 genes were unique to "yellow" gastrozooids and not found in "white" gastrozooids, and 886 genes were found in "white" gastrozooids and not found in "yellow" gastrozooids.
Finally, in A. elegans, very few significantly differentially expressed genes were identified between the B palpon and gastric palpons (one was significantly differentially expressed in B palpons and two were significantly differentially expressed in gastric palpons) (supplementary file 4, Supplementary Material online). All three genes have no significant BLAST hit and did not map to any gene trees). Genes were identified that are significantly differentially expressed in B palpons relative to all other zooids (gastric palpons were excluded). Of this, 928 genes were differentially expressed in the B palpons relative to all other zooids. Most of these genes overlapped with those differentially expressed in the gastric palpons relative to all other zooids (746 genes).

Classical Analysis of Gene Expression between Species: How Different Is Zooid-Specific Expression between Species?
Most comparative studies of gene expression focus exclusively on strict 1:1 orthologs, requiring an assessment of orthology across all species. We call this type of analysis Gene Expression across Species and Specialized Zooids in Siphonophora . doi:10.1093/molbev/msac027 MBE STF, as it limits analyses to a subset of gene trees with very specific evolutionary histories. Following this approach, we used Orthofinder 2 (v2.4.0) (Emms and Kelly 2019) to identify strict orthologs. For all seven species, we identified 1,173 strict orthologs, of which 952 orthologs had expression data across all seven species. In order to increase the number of genes for analysis, we focused on the four best sampled species (A. elegans, B. elongata, F. vityazi, and N. bijuga). Using data from four taxa, we identified 4,009 strict orthologs, of which 3,174 orthologs had expression values in four zooids/tissues: gastrozooids (developing and mature), nectophores (developing), palpons (mature), and the pneumatophore. Using TPM10K values directly, we found that orthologs clustered largely by species rather than zooid/tissue ( fig. 3A). Following Breschi et al. (2016) (see Materials and Methods), we used linear models to identify the proportion of variance that could be explained by species or zooid/tissue ( fig. 3B), and found that of 3,174 total orthologs, 2,125 are species-variable genes (SVG), and 168 are zooid/tissue-variable genes (TVG).
The strong species-dominated clustering observed when comparing expression values directly may be due to differences in counting efficiency between species, especially as we used reference transcriptomes (Dunn, Luo, et al. 2013). To account for this, we used ratios of expression values, using TMP10K expression values of the most commonly sampled zooid-mature gastrozooid-as the denominator (Dunn, Luo, et al. 2013). Using this approach, we found that orthologs clustered largely by zooid/tissue rather than species ( fig. 3C). Using linear models, we identified, out of 3,174 orthologs, 494 are SVGs, and 349 are TVGs ( fig. 3D). GO terms enriched among TVG included embryonic morphogenesis, embryo development, cartilage development, and a number of metabolic and biosynthetic processes. Meanwhile, SVG were enriched in GO terms such as cellular response to stress, DNA repair, and cell cycle processes. A list of TVG and SVG can be found in supplementary files 7 and 8, Supplementary Material online, respectively.
Phylogenetic Analysis of Gene Expression with SBF: What Are the Evolutionary Changes in Zooid Expression along Branches in the Species Tree? For our SBF analyses (see fig. 2; Materials and Methods), we used transcriptome and genome data from 41 cnidarian species to generate a total of 7,070 gene trees, of which 3,831 gene trees passed filtering criteria (see Materials and Methods for filtering approach). We used Orthofinder 2 (Emms and Kelly 2019) to reconcile the species tree and the gene trees to annotate each gene tree node as either a speciation or duplication event. The number of genes represented in these gene trees is shown in supplementary table S1, Supplementary Material online. The internal nodes on these gene trees represent 20,088 speciation events and 9,082 duplication events. Expression values (TPM10K) were mapped to the tips of the gene trees for each zooid/tissue separately, and ancestral expression values were inferred at internal nodes in a maximum likelihood framework. We focused on a subset of zooids and tissues that are common across the sampled species: gastrozooids (developing and mature), nectophores (developing), palpons (mature), and the pneumatophore. The distribution of expression changes along species-equivalent branches are shown in supplementary figure S13, Supplementary Material online.
As with the STF and linear model analyses, we investigated the impact of counting efficiency on our comparative analyses. We used the TPM10K values at the tips and nodes to calculate expression ratios and calculated changes in expression ratios across a branch for pairs of tissues, using the mature gastrozooid values as the denominator. The distribution of expression ratio changes along species-equivalent branches (branches in the gene tree that correspond to equivalent branches in the species tree, and are descended from speciation events) are shown in supplementary figure S14, Supplementary Material online. We found that the variance of change across a given species-equivalent branch is considerably higher for raw TPM10K values as compared with ratios of expression (supplementary figs. S15 and S16, Supplementary Material online). We also identified the number of gene tree branches that had a negative (change À2), positive (change !2), or neutral (change >À2 and <2) change in either TPM10K or TPM10K ratios across the branch ( fig. 4). The mean number of species-equivalent branches considered in these analyses are shown in table 1-this value reflects the mean number of branches across each zooid or zooid ratio (due to incomplete sampling, some zooids may have fewer or more branches; as seen in fig. 4). For each species-equivalent branch, we found more negative and positive branches than neutral branches when we used TPM10K values, whereas for expression ratios, we found that the majority of branches show no change across the branch. This suggests that as for the STF analyses, counting efficiency has a large impact on expression change and in turn leads to large gene/branch-specific signal. For expression ratios, these results indicate that the vast majority of branches show close to 0 (neutral) change across the branch, suggesting that for closely related genes, expression ratios tend not to differ.
Out of a total of 3,357 final gene trees considered in these analyses, 3,329 gene trees contained branches with neutral changes in expression ratios, 1,294 gene trees contained branches with positive changes in expression ratios, and 1,041 gene trees contained branches with negative changes in expression ratios. Expression ratio changes across speciesequivalent branches from all gene trees is available in supplementary file 9, Supplementary Material online (note: in this file, "BLAST hit" is the most frequent BLAST hit for the gene tree, gene identity should still be confirmed for each gene, particularly in large gene trees).
The vast majority of gene trees contained speciesequivalent branches with neutral changes in expression ratios, including a number of transcription factors and morphogenic signaling pathway genes. For example, in ratios of developing and mature gastrozooids, among the relevant speciesequivalent branches in the Wnt gene phylogeny (Wnt gene identity based on Condamine et al. [2019]), changes were neutral across Wnt3 and Wnt2 branches, with the exception of a slight positive change across the branch leading to a Diphyes dispar Wnt2 paralog, suggesting a higher relative expression of this gene in the developing gastrozooid of Munro et al. . doi:10.1093/molbev/msac027 MBE D. dispar ( fig. 5 and supplementary fig. S17, Supplementary Material online). For all zooid ratios, we found very consistent Wnt3 expression patterns with neutral or very small positive or negative changes in expression ratio across the branches. Indeed, across all cnidarians examined, Wnt3 has consistent localized expression at the oral pole, likely playing a role in axis formation and maintenance (Hobmayer et al. 2000;Guder et al. 2006;Momose et al. 2008;Nawrocki and Cartwright 2013;Hensel et al. 2014;Bagaeva et al. 2019). However, for Pneumatophore/ Gastrozooid ratios, some large changes were seen across branches K and L, leading to an A. elegans and N. bijuga Wnt2 paralog, respectively, likewise, for Palpon/Gastrozooid, a large change was seen across branch J, leading to a F. vityazi Wnt2 paralog. In Nematostella vectensis and Hydractinia echinata, Wnt2 is expressed in the middle of the polyp (Kusserow et al. 2005;Hensel et al. 2014). Without more detailed spatial expression patterns, it is difficult to know whether these differences in expression between species reflect speciesspecific differences in axial patterning and morphogenesis, such as an expansion or reduction of the expression domain.
We also looked at the 277 gene trees with species-relative branches that had very large changes (5 < change < À5), representing putative lineage-specific expression patterns   . Bottom right panel shows number of branches (ordinate) with negative, neutral, or positive (colors) changes in gene expression ratios for different tissues/zooids (abscissa) across each branch in the species tree (species-equivalent branches). Neutral change (green) is defined as 2 and !À2, negative change (red) is defined as <À2, and positive change (blue) is defined as >2. These units correspond to TPM10K or TPM10K ratios/substitutions per site. GasdevGasmat, developing gastrozooid/mature gastrozooid; NecdevGasmat, developing nectophore/mature gastrozooid; PalmatGasmat, mature palpon/mature gastrozooid; PneGasmat, Pneumatophore/mature gastrozooid. MBE across branches L (leading to N. bijuga) and J (leading to F. vityazi) and negative change across branch K (leading to A. elegans). In ratios of developing nectophores/gastrozooids, we also found very large positive changes in branch J for Wnt-7blike, Wnt-4, (putatively) Hox-B8 and NKX1-2. Likewise, in pneumatophore/gastrozooid ratios, we also found large positive changes in branches J and L for Frizzled5/8, and a negative change across branch K. Across J, we also saw positive changes in pneumatophore/gastrozooid in Hox-B8 and NKX1-2. Finally, for branch K, we saw large positive changes in Wnt4 MBE and Wnt7b-like. Many of these genes have been shown to have specific expression domains patterning cnidarian bodies (Kusserow et al. 2005;Ryan et al. 2007;Sinigaglia et al. 2013;Hensel et al. 2014), and further detailed expression analyses of these genes and others within the zooid are required to determine whether the patterns identified here reflect differences in expression domain between species and zooid. The nature of these analyses makes it difficult to investigate GO term enrichment, as each gene tree contains multiple gene lineages often with diverse molecular function and that play diverse roles in biological processes-additionally, our approach focuses on species-relative branches rather than tips. Nevertheless, we assigned GO terms at a gene tree level, and investigated GO term enrichment among gene trees that contain branches with negative and positive changes. Among gene trees that contained species-relative branches with positive changes, we found an enrichment for a number of biological processes, including detection of stimulus, neuroblast proliferation, stem cell proliferation, and axon guidance (supplementary files 10 and 11, Supplementary Material online). Among gene trees that contained species-relative branches with negative changes, we found in particular an enrichment for a number of metabolic processes including amine, NAD, tryptophan, and indolalkylamine metabolic process (supplementary files 12 and 13, Supplementary Material online).
To assess the extent to which missing data may impact these analyses, we identified BUSCO scores for each of the reference transcriptomes (supplementary table S2, Supplementary Material online). BUSCO completeness score ranged from 88.7% to 70.2% against the Metazoa BUSCO data set (Manni et al. 2021). Of this, 92.65% of identified BUSCO genes are present in the initial expression data set. This suggests that our de novo transcriptome assemblies captured the vast majority of known metazoan genes, and therefore are likely not sensitive to missing data. After filtering for the ortholog STF analyses, 38.65% of identified BUSCO genes were retained. By contrast, in the SBF analyses, 44.02% of identified BUSCO genes were retained. These findings show that our novel phylogenetic approach SBF may improve evolution of gene expression studies as it retains more gene trees for downstream analyses.
The vast majority of zooid libraries mapped well to the reference transcriptomes, with the majority of libraries having between 85% and 95% of reads aligned to the reference, with an average of 85.84% of reads aligned (supplementary file 14, Supplementary Material online). Only eight libraries had poor alignment scores (<70%), notably all three developing nectophore libraries from N. bijuga (56.53-63.08%), male and female gonodendron libraries from a single F. vityazi replicate (68.18% and 68.89%, respectively), an Apolemia lanosa developing nectophore replicate (60.21%), and a developing bract and pneumatophore replicate in B. elongata (64.46% and 49.61%, respectively). The fact that all three N. bijuga developing nectophore replicates had lower mapping scores, points to a reduced presence of likely nectophore-specific genes in the N. bijuga reference transcriptome, which also had the lowest BUSCO score of all of the reference transcriptomes.

Evolution of Gene Expression in Siphonophores
In this study, we used RNA-seq to investigate the functional specialization and evolution of zooids across siphonophore species. In our analysis of differential expression within species, we found a large number of differentially expressed genes across siphonophore zooids, reflecting the distinct anatomy and function of these zooids. In addition, we identified potential transcription factors that are significantly differentially expressed within particular zooids, with potential homologs found in several species, which are interesting candidates for future study in siphonophores or related cnidarian colonial groups. Using these within-species DGE analyses, we identified distinct expression patterns that may be used to define particular zooid types.
We also explored gene expression patterns in three zooid types that are unique to three species, P. physalis, A. elegans, and B. elongata. In siphonophores, different zooid types are typically defined based on morphological and functional differences/similarities, and on the location of the zooid within the colony (which is determined in most species by the asexual budding process in the growth zone, which gives rise to a repeating pattern of zooids along the stem) (Totton 1965;Mackie et al. 1987;Dunn 2005;Dunn and Wagner 2006). Based on both morphology and development, P. physalis tentacular palpons are considered to be a distinct and unique zooid type (Munro et al. 2019), and the DGE data indicate that there are clear morphological and functional differences between gastrozooids and tentacular palpons. We also identified several putative toxin genes with distinct expression profiles between these two zooids, which matches tissuespecific venom observed in other cnidarian species (Ames et al. 2016;Macrander et al. 2016;Klompen et al. 2021).
For "yellow" and "white" gastrozooids in B. elongata, the differential expression data point to few morphological or functional differences between these two zooids. Through pairwise differential expression patterns between "yellow" or "white" gastrozooids and other zooids within the colony, we were able to identify hundreds of genes that are significantly expressed in one zooid and not the other. This suggests that these two gastrozooids may be distinct zooid types, although they are functionally very similar to one another. Sequencing at greater depth, and also functional work within this species may help clarify the nature of these differences. By contrast, there is no strong evidence in our data that the B palpon and gastric palpons in A. elegans are sufficiently different to constitute a novel zooid type. These findings suggest that location within the colony is not necessarily sufficient to designate a novel zooid type.
With our STF and linear model analyses, we found that the vast majority of identified orthologs were neither exclusively species-or tissue-variable, with only a subset of genes being either tissue-or species-variable. As has been observed for vertebrate organs, we find based on GO terms that identified Munro et al. . doi:10.1093/molbev/msac027 MBE SVGs tend to play a role in basic cellular functions (so-called housekeeping genes), as compared with TVGs (Breschi et al. 2016).
What can we learn from the SBF analyses, using ratios of expression data? It is important to stress that changes in expression ratios cannot tell us about changes in expression magnitude. A TVG (identified via linear models in classical analyses) that is either highly or lowly expressed in, for example, nectophores relative to gastrozooids in all species, may show largely neutral changes across all species-equivalent branches. The sign (positive or negative) of change provides an indication of whether expression in the child node is relatively higher in the denominator or the numerator, relative to the parent node. We find overall that most expression ratios show very little change across species-equivalent gene tree branches, suggesting that gene expression patterns tend to be largely consistent among species. Positive or negative changes across branches, especially very large changes, represent putative linage-specific shifts in expression. Whether these reflect distinct morphological or functional changes in the tissues of a particular species or clade requires further validation-a difficulty in this system, where the species are difficult to collect and maintain in the lab.

Challenges and Solutions in the Analysis of the Evolution of Gene Expression
In this study, we used three different analyses to investigate gene expression patterns within homologous zooids across species, each are complementary and enable different possibilities for biological discovery. Within-species analyses focus on expression among zooids and tissues within a single species-a very large number of genes identified in the transcriptome or genome can be considered, and the broader across species gene-tree need not necessarily be considered. This approach is especially useful for investigating expression patterns within novel zooids that do not have clear homologs in other species. The classical between species analysis as implemented here, using STF followed by linear models, relies on phylogenetic assumptions yet is a nonphylogenetic tipfocused approach where a single gene is identified per species, and where the identified differences in values across tips are due to changes along all branches in the tree. With linear models, we can identify genes with expression values that are common to a particular zooid across all species, or with expression that is specific to a particular species. However, the vast majority of genes are neither zooid nor species variable, and show more complex patterns of expression. We proposed a third approach, called SBF, that gives access to changes in expression across gene phylogenies, which may differ significantly from the species phylogeny.
Where traditional approaches identify strict ortholog genes before conducting analyses, SBF uses information from all gene copies within a gene tree regardless of their evolutionary history of duplication or speciation. STF is a gene tree filtering approach, where entire trees or subtrees are discarded due to even a single duplication event, SBF is a branch filtering approach. This means that SBF retains many informative branches from trees that would be removed by STF, and preserves many more evolutionary comparisons for analyses. By filtering our data by branches, we consider expression patterns at both the tips and internal nodes of gene phylogenies. We are thus able to identify shifts in expression leading to tips as well as shifts in expression leading to particular clades. Although we focus here on expression following speciation events, it is also possible to compare expression following duplication and speciation events.
Using ancestral trait reconstruction, we also overcame sampling issues at the tips of the gene and species phylogenies, as expression values of different treatments can be reconstructed at deep internal nodes, even where there may be inconsistent sampling at the tips. That is, even if expression values are missing for a given tissue and gene in a gene phylogeny, we are still able to examine expression patterns within this gene tree. By contrast, in the STF analyses, incomplete sampling in the expression matrix leads to the elimination of the ortholog from the analysis.
As with all methods that rely on mapping to reference transcriptomes rather than genomes (including STF analyses), this approach is limited by the quality of the reference transcriptomes. Ratios of expression helped significantly to improve issues of differences in count efficiency among species (Dunn, Luo, et al. 2013). However, reference transcriptome quality also has an impact on the quality of the gene trees used for expression mapping. Not all reference transcriptomes were sequenced to equal depth among species (supplementary table S2, Supplementary Material online), and this has important effects on the presence or absence of genes from particular species within the gene tree, as well as within the broader expression data set. This not only has an effect on the representation of expression values, but also impacts the power to investigate patterns of expression among branches within a gene tree. With genome sequencing becoming cheaper and more readily available, the widespread availability of reference genomes will help alleviate many of these issues. Reference genomes will also improve gene models, enabling the distinction of different alleles of the same gene from duplicated gene copies, this in turn will improve the quality of the gene trees. However, gene loss will nevertheless present a challenge to these analyses.

Conclusions
With the expansion of functional genomic tools, including RNA-seq and single cell-sequencing methods, there is considerable interest in looking not only at how genomic variation gives rise to phenotypic diversity in a single species or organism, but also at how functional genomic variation shapes phenotypic diversity across multiple closely and distantly related species to understand broader evolutionary patterns and processes (Brawand et al. 2011;Barbosa-Morais et al. 2012;Merkin et al. 2012;Perry et al. 2012;Yang and Wang 2013;Necsulea et al. 2014;Zhang et al. 2014;Sudmant et al. 2015;Breschi et al. 2016;Levin et al. 2016;Macrander et al. 2016;Clarke et al. 2017;Ma et al. 2018;Cardoso-Moreira et al. 2019; Darbellay and Necsulea 2020; Fukushima and Pollock 2020). Many of these analyses have the goal of identifying Gene Expression across Species and Specialized Zooids in Siphonophora . doi:10.1093/molbev/msac027 MBE shared expression patterns among modular biological units (cell type, tissue, organ, zooid) across species, in order to identify commonalities in expression patterns among species. This is of particular interest for medically orientated fields interested in understanding the extent to which expression results can be extrapolated from one model organism to another. Another goal is to identify expression patterns in a particular biological unit that are unique to a particular species or even clade. For these questions, within-species analyses provide the greatest depth, in terms of number of genes investigated, however comparisons between species on the basis of within-species DGE are limited and largely qualitative. Here, we use two approaches to investigate expression patterns in a quantitative manner across homologous tissues. Classic between species analyses, using STF in conjunction with linear models, enabled the investigation of a smaller number of strict orthologs that vary in a tissue-or speciesspecific manner, however as this analysis focuses on the tips of gene trees, our ability to investigate lineage or especially cladespecific patterns of expression in a given zooid are more limited. Meanwhile, phylogenetic analysis using SBF focuses on branches rather than tips, also with specific evolutionary histories (descended from speciation events), but it enables the identification of expression patterns that vary little among genes/species, as well as expression patterns that show strong lineage-or clade-specific patterns for a given zooid.

Materials and Methods
Collecting Specimens were collected in the north-eastern Pacific Ocean in Monterey Bay and, in the case of P. physalis, the Gulf of Mexico. Specimens were collected by remotely operated vehicle or during blue-water SCUBA dives. Physalia specimens were collected by hand from the beach after being freshly washed on-shore by prevailing winds. Available physical vouchers have been deposited at the Peabody Museum of Natural History (Yale University), New Haven, CT (supplementary file 15, Supplementary Material online). Specimens were relaxed using 7.5% MgCl 2 hexahydrate in Milli-Q water at a ratio of 1/3 MgCl 2 and 2/3 seawater. Zooids were subsequently dissected from the colony and flash frozen in liquid nitrogen. Colonies were cooled to collection temperatures (e.g., 4 C for deep sea species) while the dissections took place. Dissections took no longer than 15-20 min. In the case of large colonies, the stem was cut and only partial sections of the colony were placed under the microscope at a given time. Each replicate individual represents a genetically distinct colony from the same species. Replicate specimens were of an equivalent colony size, and zooid replicates were also equivalent sizes. Larger zooid types, such as gastrozooids, were sampled as a single zooid, but smaller zooids were pooled. Pooled zooids were of a comparable maturity and sampled from the same location in a single colony. Sampling data, including time, date, depth, and voucher ID, can be found in supplementary file 15, Supplementary Material online. Sequencing mRNA was extracted directly from tissue using Zymo Quick RNA MicroPrep (Zymo No. R1050), including a DNase step, and subsequently prepared for sequencing using the Illumina TruSeq Stranded Library Prep Kit (Illumina, No. RS-122-2101). 50 base-pair single-end libraries were all sequenced on the HiSeq 2500 sequencing platform. Three sequencing runs were conducted, representing three full flow cells. To avoid potential run/lane confounding effects, where possible, libraries of multiple zooids/tissue of a single individual in a species were barcoded and pooled in a single sequencing lane, and replicate lanes of zooids/tissue from different individuals of the same species were sequenced in separate runs. Additionally, two libraries were run as technical replicates across all runs and many lanes, for a total of 20 technical replicates.

Analysis
Differential Gene Expression Short-read libraries were mapped to previously published transcriptomes (Munro et al. 2018) using Agalma v 2.0.0 (Dunn, Howison, et al. 2013;Guang et al. 2021), which uses a number of existing tools for transcript quantification, including RSEM (which uses Bowtie) (Langmead et al. 2009;Li and Dewey 2011). Using the agalmar package (https://github. com/caseywdunn/agalmar, last accessed February 7, 2022), we filtered out genes that were flagged as being rRNA, and selected only protein coding genes. We also only considered genes that had greater than 0 counts in at least two libraries. DGE analyses, including normalization, were conducted in R, using the DESeq2 package (Love et al. 2014). Libraries that were found to be outliers based on mean Cook's distance were removed from the DESeq object and from downstream analyses and normalization. Testing for differential expression was conducted using the results() function in DESeq2. Genes were considered to be significantly differentially expressed if adjusted P values (Bonferroni correction) were less than 0.05. Differential expression analyses were only conducted on zooids/tissue with two or more replicates.
GO annotations were retrieved for each of the referencetranslated transcriptomes (Munro et al. 2018) using the PANNZER2 web server (Törönen et al. 2018). The PANNZER2 format was modified to match the gene2GO format required for the package topGO (Alexa and Rahnenfuhrer 2020). Gene set enrichment analyses were carried out within species using the R package GOseq (Young et al. 2010), which takes gene length into account. Over and underrepresented categories were calculated using the Wallenius approximation, and P values were adjusted using the Benjamini and Hochberg method. Categories with an adjusted P value below 0.05 were considered enriched. Gene set enrichment analyses were also conducted at the gene tree level, considering representative GO terms for particular gene trees. Representative GO terms were selected based on the frequency of occurrence among genes in the gene tree. As gene lengths vary among species and genes in the gene tree, the GOseq approach could not be used, and topGO was used to detect GO terms that are enriched based on Fisher's exact test. This approach assumes that each gene Munro et al. . doi:10.1093/molbev/msac027 MBE tree has an equal probability of having genes shared among species that are detected as differentially expressed, however results may be biased by a number of factors, including mean gene length among genes in the gene tree (Young et al. 2010

TPM10K
For all comparative expression analyses, expression values were normalized using a new method we call TPM10K. For gene i of a given species, TPM is typically calculated as (Li et al. 2010): where h i is the number of the mapped reads to gene i, ' i is the effective length of the gene, and n is the number of genes in the reference. The intent of this measure is to make libraries comparable within a single species. The sum of TPM values within a library is 10 6 , and the mean is 10 6 n . One implication of this is that TPM values are not directly comparable across species, since in practice, n differs across species. If this were not accounted for, then it could appear, for example, that genes all have lower expression in a species with a more complete reference transcriptome and higher n. To account for differences in means among species, we use a new measure, TPM10K, that accounts for differences in n: where the sum of TPM10K values within a library is 10 2 Â n and the mean is 10 2 . By multiplying by n; we are able to account for different sequencing depths among species, and ensure a common mean. As n is large, we divide by an arbitrary number (in this case 10 4 ) in order to reduce the magnitude of the expression value.

Species Tree Filtering
For STF analyses, single copy orthologs were obtained using Orthofinder 2 (v2.4.0) (Emms and Kelly 2019). Linear models used in STF analyses were constructed using lm(), following methods and code developed by Breschi et al. (2016), https:// github.com/abreschi/Rscripts/blob/master/anova.R (last accessed February 7, 2022). All SVGs and TVGs are genes for which both species and tissue explain greater than 75% of variance. Additionally, in SVGs, the proportion of variance explained by species is two times greater than that explained by tissues; whereas in TVGs, the proportion of variance explained by tissues is two times greater than that explained by species.

Species Branch Filtering
The principles of SBF are as follows: first, we identify speciation and duplication events in a given gene tree, specifically labeling nodes that correspond to speciation events in the species tree ( fig. 2, step 1). Next, we map gene expression values to the tips of the gene tree ( fig. 2, step 2), and use phylogenetic methods to reconstruct expression values at the internal nodes ( fig. 2, step 3). Expression values are mapped and reconstructed for each zooid/tissue separately, with the assumption that the structure is homologous across species. Then, we calculate scaled expression values across branches, by subtracting the expression value at the child node from the parent node, and dividing by branch length (branch lengths are calibrated against branches in the species tree) ( fig. 2, step  4). Finally, we identify branches in the gene trees that correspond to equivalent branches in the species tree (hereafter species-equivalent branches) ( fig. 2, step 5). Speciesequivalent branches are gene tree branches with parent and child nodes that are both speciation events and that correspond to branches in the species tree. Each branch in the species tree is given a unique identifier (i.e., a letter) and the species-equivalent branches in gene trees are given the same identifier ( fig. 2, step 5). This method enables the selection of branches within gene trees that are equivalent to branches within the species tree, and are thus comparable with one another across all gene trees. Unlike the STF approach, this approach considers equivalent branches that are descended from speciation events, but that have more complex evolutionary histories. For example, due to deeper gene duplication events, gene trees often contain multiple branches that correspond to the same branch in the species tree ( fig. 2, step 5). Our method allows us to consider all of these branches.
Gene trees were generated using the transcriptomes and genomes from 41 species (Munro et al. 2018) using Agalma v 2.0.0 (Dunn, Howison, et al. 2013;Guang et al. 2021). Although we have expression data for seven species, we used a broader species sampling in order to infer more complete gene histories and more easily assign speciation and duplication events. Following the treeinform step of Agalma v 2.0.0, amino acid data were exported and supplied to Orthofinder 2 (v2.3.8) (Emms and Kelly 2019) for simultaneous co-estimation of gene trees with the published maximum likelihood species tree (Munro et al. 2018). Within Orthofinder 2, the selected multiple sequence alignment method was MAFFT (Katoh and Standley 2013) and maximum likelihood tree inference method was IQ-tree with the LGþFþR4 substitution model (Nguyen et al. 2015). BUSCO scores were generated from the reference transcriptomes using the metazoa_odb10 BUSCO data set with BUSCO v 5.2.2 (Manni et al. 2021).
Gene Expression across Species and Specialized Zooids in Siphonophora . doi:10.1093/molbev/msac027 MBE For each gene tree, we generated summary statistics of the branch lengths in the tree (maximum, minimum, SD) as well as the fraction of branches that have a default length value of 10 À6 (that is indicative of branch length ¼ 0). The majority of gene trees have branches with a maximum length around 1. Gene trees were filtered to exclude trees with branches whose maximum length is >2 and that had more than 0.25 branches with the default length value, representing 15.37% of gene trees. The goal of this is to exclude trees that include very long branches, or a large proportion of branches with very short branch lengths. Gene tree nodes were annotated as speciation events or duplication events, based on assignment by Orthofinder 2 (Emms and Kelly 2019). Speciation nodes were subsequently assigned a node ID equivalent to the species tree node, using species tip names from the gene tree to determine the most common recent ancestor in the species tree (Munro et al. 2018), using the phytools package. Due to the use of the species-overlap method by Orthofinder 2, some clades of single copy genes were assigned as speciation events, although the topology is inconsistent with the species tree. To avoid time calibration issues, due to descendant nodes being assigned the same species node ID, descendant speciation nodes with the same species node ID were marked as null, and gene trees with nodes greater than 0.3 null nodes per internal node were excluded-indicating widespread topological differences with the species tree. Additionally, only trees with one or more speciation events were retained, as speciation events are used for time calibrations. The gene trees were then time calibrated to the species tree using chronos() in the ape package, so that the branch lengths were scaled to the same equivalent length across all gene trees (Paradis et al. 2004). Some gene trees could not be calibrated against the node constraints from the species tree and were discarded. Additionally, we calculated the maximum node age for each gene tree, and excluded gene trees with roots that exceeded a maximum root depth of 5, which can be indicative of calibration problems (the root of the species tree is 1). Only a single tree was excluded as a result of a very deep root. Tips without expression values were then pruned out of the tree. Gene trees with fewer than three expression values at the tips were discarded, retaining only trees with three or more values. Pruned and unpruned calibrated gene tree files are available as Supplementary Material online.
We then took the mean TPM10K value for each gene across replicates of the same zooid/tissue within a species and applied a log transformation. Using gene trees with expression values for each gene within a species at the tips, maximum likelihood ancestral trait values were generated at the nodes using the anc.recon() function in Rphylopars assuming a Brownian motion (BM) model of evolution ( fig. 2, steps 2 and 3) (Goolsby et al. 2017). As not all zooids are present in all of the species, the trees were pruned down to the subset of tips with expression values for ancestral trait reconstructions. Node values were then added back to the unpruned tree with all of the reconstructed expression values. Change in expression was measured across a branch by taking the difference between a parent node and a child node, and then this difference is scaled by branch length (fig. 2, step 4).
To examine if expression values in our data set evolve under BM, we simulated data set of random expression values using fastBM() from the phytools package with empirically derived mean and SD values for each gene tree. Under this model, expression variance accumulates among linages on the gene tree as a function of time, and is used as a model of evolution under drift, as well as some forms of natural selection (Felsenstein 1973;).

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.
available as a supplementary zipped folder, Supplementary Material online.