A broadly resolved molecular phylogeny of New Zealand cheilostome bryozoans as a framework for hypotheses of morphological evolution

Larger molecular phylogenies based on ever more genes are becoming commonplace with the advent of cheaper and more streamlined sequencing and bioinformatics pipelines. However, many groups of inconspicuous but no less evolutionarily or ecologically important marine invertebrates are still neglected in the quest for understanding species- and higher-level phylogenetic relationships. Here, we alleviate this issue by presenting the molecular sequences of 165 cheilostome bryozoan species from New Zealand waters. New Zealand is our geographic region of choice as its cheilostome fauna is taxonomically, functionally and ecologically diverse, and better characterized than many other such faunas in the world. Using this most taxonomically broadly-sampled and statistically-supported cheilostome phylogeny comprising 214 species, when including previously published sequences, we tested several existing systematic hypotheses based solely on morphological observations. We find that lower taxonomic level hypotheses (species and genera) are robust while our inferred trees did not reflect current higher-level systematics (family and above), illustrating a general need for the rethinking of current hypotheses. To illustrate the utility of our new phylogeny, we reconstruct the evolutionary history of frontal shields (i.e., a calcified bodywall layer in ascus-bearing cheilostomes) and asked if its presence has any bearing on the diversification rates of cheilostomes.


Introduction
Large and broadly-sampled phylogenies are vital to robustly answering many different classes of evolutionary questions, including those involving trait evolution, origins and evolution of biogeographic distributions and rates of taxonomic diversification. While megaphylogenies with hundreds to thousands of species (Smith et al., 2009) are available for many groups of vertebrates (Meredith et al., 2011;Prum et al., 2015) and plants (Zanne et al., 2014), and also for some non-vertebrate terrestrial groups (Varga et al., 2019), the molecular phylogenetics of many marine invertebrate groups remains relatively neglected (Arrigoni et al., 2017;Kocot et al., 2018;O'Hara et al., 2017).
In this contribution, we begin to rectify the paucity of large and/or taxonomically broadly sampled molecular phylogenies for marine invertebrates, targeting a phylum whose rich fossil record can be subsequently integrated for evolutionary analyses. Our focal group is including 867 cheilostomes (of which 285 species remain to be formally described). About 61% of New Zealand's marine Bryozoa are endemic (Gordon et al., 2019), making New Zealand a doubtless diversity hotspot for cheilostome bryozoans. Complementing Recent diversity, the published Cenozoic record of cheilostome bryozoans is also rich, though relatively less studied (Brown, 1952;Gordon and Taylor, 2015;Rust and Gordon, 2011), comprising 531 species (of which 240 are in open nomenclature). This complementarity of living and fossil species renders a molecular phylogeny of New Zealand taxa amenable to modern statistical methods that integrate molecular and fossil data for inferring evolutionary processes (Heath et al., 2014). Last, but not least, New Zealand is one of the better-studied marine regions taxonomically and ecologically for Bryozoa (e.g. Gordon, 1984Gordon, , 1986Gordon, , 1989Gordon et al., 2009;Schack et al., 2020), a phylum that is somewhat neglected in many other parts of the world. Bryozoan research has been continuously conducted in New Zealand since 1841 ) and a governmental agency, the National Institute of Water and Atmospheric Research (NIWA), is both the data manager and custodian for fisheries and invertebrate research data, hence assuring knowledge curation. All of this means that a cheilostome phylogeny with New Zealand species broadly represented allows us to begin to ask evolutionary and ecological questions while controlling for phylogenetic nonindependence.
Here we apply a genome-skimming approach to New Zealand cheilostome bryozoans and present a robustly supported molecular phylogeny based on 15 mitochondrial and 2 rRNA genes. The molecular sequences of 199 cheilostome colonies sampled in New Zealand are presented here for the first time. Using 180 species and 96 genera from New Zealand and previously sequenced, non-New Zealand species, we construct the largest and most taxonomically broadly sampled cheilostome phylogeny to date, with 263 in-group colonies, representing 214 species and 120 genera. The inclusion of non-New Zealand taxa allows us to explore the robustness of the inferred relationships among New Zealand species but also reduces phylogenetic inference errors by nature of a broader taxonomic sampling (Pollock et al., 2002). To illustrate the utility of our inferred tree for understanding cheilostome evolution, we reconstruct the evolutionary history of a morphological trait (the calcified frontal shield) and ask if its evolution might have changed the diversification rates of cheilostomes that have such a shield (ascophoran-grade) versus those that do not (anascangrade). We also discuss several other key taxonomic traits widely thought to be evolutionarily stable and the consequences our highly resolved cheilostome phylogeny has for these. Our contribution is a first step towards a global cheilostome megaphylogeny, needed for answering biological questions that go beyond those probing genealogical relationships.

Sampling & SEM
Sequences are provided here for 207 New Zealand cheilostome colonies that were collected during several field expeditions by NIWA and University of Otago, New Zealand. While we have newly sequenced 199 colonies, we also supply unpublished sequences for 8 extra colonies we previously presented (see Supplementary Table S2). Samples were sorted, preserved in 70-96% ethanol, then shipped to the University of Oslo, Norway, for processing.
Each bryozoan colony, preliminarily identified to the lowest possible taxonomic level (usually genus but sometimes species) using a stereoscope, was subsampled for DNA isolation, and also for scanning electron microscopy. The scanning electron micrographs (SEMs), taken with a Hitachi TM4040PLus after bleaching to remove tissue (where appropriate), are required for species-level confirmation. All SEM digital vouchers are supplied as a supplementary data file. Taxonomic identifications are made independently of the phylogenetic inference and metadata to avoid identification bias.

DNA isolation, sequencing and assembly
The 199 subsamples of colonies (henceforth "samples") were dried before genomic DNA isolation using the DNeasy Blood and Tissue kit (QIAGEN, Germantown, MD, USA).
Samples were homogenized in lysis buffer, using a pestle, in the presence of proteinase-K.
Genomic DNA were sequenced at the Norwegian Sequencing Centre (Oslo, Norway) using Illumina HiSeq4000 150 bp paired-end (PE) sequencing with a 350 bp insert size.
Congruence between the topological signal of the ML and Bayesian trees for both the New Zealand (Supplementary Figs. S1 and S2) and global (Supplementary Figs. S3 and S4) inferences was tested using the I cong index (de Vienne et al., 2007).

Ancestral state reconstruction and BiSSE analyses
The tips states of whether the sampled species is anascan (0, having a non-calcified frontal membrane) or ascophoran (1, having a calcified frontal shield), both states decipherable from SEMs, is given in in Fig. 3. We use a standard Markov model of binary character evolution (Pagel, 1994) implemented in ape (Paradis and Schliep, 2018) to estimate the ancestral states of the nodes on our inferred phylogeny. We use a standard binary state speciation and extinction model (Maddison et al., 2007) implemented in diversitree (FitzJohn, 2012) to investigate any differences in diversification rates due to the anascan or ascophoran frontal shield state of the species involved. As input for this latter analysis, we estimate that of the 1876 anascans and 3358 ascophoran species in Bock 2020, we have sampled 4.4% and 3.9% respectively to account for biases due to the sampling of species given the trait. We perform ancestral state reconstruction and BiSSE analyses for both ML and Bayesian "global" trees (Supplementary Figs. S3 and S4 respectively) to account for minor differences in the topological signal (see Results). In cases where there are multiple representatives within a species, we choose the colony with the highest number of nucleotides/amino-acids/genes to represent the species for these analyses.

Sequencing and concatenation
We successfully sequenced and assembled 199 New Zealand cheilostome colonies, representing 165 species (SEM vouchers in Supplementary file) that have never been presented previously (Supplementary Table S1). We supply additional sequence data for a further eight species previously presented (Supplementary Table S2 and Orr et al., 2019b).
The final 17 gene and 267 taxa "global" supermatrix constitutes 77% total character completeness for the dataset used to infer Fig. 2. For the convenience of future workers interested in only the New Zealand taxa, we supply also trees based on these data (Supplementary Figs. S1 (ML) and S2 (Bayesian), where character completeness is 78%).
The assembled rRNA and mitogenomes are deposited at NCBI with accession numbers (Supplementary Table S2).

Broad taxon-sampling
Our inferred "global" cheilostome phylogeny, encompassing 214 species and 120 genera, from 56 families ( Fig. 2) of which 229 colonies, 186 species and 96 genera, currently distributed in 48 families, are from New Zealand (Supplementary Figs. S1 and S2). The New Zealand and global trees represent c. 21% described species of cheilostomes from New Zealand and c. 15% of the described cheilostome genera globally, respectively. Both phylogenies ( Fig. 2 and Supplementary Fig. S1) are robustly resolved with most branches and relationships receiving either high (>90 bootstrap (BS) / >0.99 Posterior Probability (PP)) or full support (100 BS / 1.00 PP). Our ingroup cheilostome taxa form a fully supported monophyletic clade, when we infer the global tree including a ctenostome outgroup (Fig. 2).
We summarize only general ingroup observations while referring the reader to topological details in Fig. 2 and Supplementary Fig. S3 that are not discussed here or in the Discussion.
We also refrain from summarizing results above the family-level for reasons stated in the Discussion.

Family relationships
Several families for which we have four or more genera represented form supported monophylies (Fig. 2)

Genus relationships
In contrast to family-level systematics, the 50 currently morphologically defined genera for which we have two or more representatives in general form monophyletic groupings (e.g. Because there are indications that some species are phenotypically highly variable and others have morphologies that are not yet well-understood, we also sequenced multiple colonies of the same species in several cases even though our goal was to sequence one colony of each species. Morphologically identified species match genetic species inferred by phylogenetic inferences in these cases, including Parasmittina aotea, Parkermavella punctigera, Chiastosella longaevitas, C. enigma, Microporella agonistes and M. intermedia. For more details, please see individual SEM cards in the Supplementary data file.

Congruent trees and a single incongruent branch
We show the ML and Bayesian trees for both the New Zealand (Supplementary Figs. S1 and S2) and global (Supplementary Figs. S3 and S4) inferences to have more congruent topologies than expected by chance with an I cong index of 38.72 and 10.05, respectively.
The probability that they are topologically unrelated are 0 and 2.43e-124 for the New Zealand and global trees, respectively.

We highlight the incongruent placement of the Euthyroides, Figularia and Valdemunitella
clade. Note that this clade is highly supported as a monophyly in both sets of trees, but its placement within the trees is contested; the ML trees, whether based only on the New Zealand taxa or all taxa (Fig. 2, Supplementary Figs. S1 and S3), place this clade in a basal position with an affinity to the Macropora/Monoporella grouping. The Bayes trees, however, infer a more derived position. In all instances (ML and Bayes), support for the inferred placement is lacking.

Ancestral state reconstruction and BiSSE
A different rates model for the transition of the anascan to ascophoran state has a less negative log-likelihood (-29.24) than that for an equal rates model (-35.16), suggesting that it describes our ML tree better. Parameter estimates indicate that the ascophoran state never goes to anascan, and anascan state goes to ascophoran at rate of 0.207 (std err 0.0273), in our ML tree. The estimated node states are shown in Fig. 3. Plots of posterior distributions of speciation and extinction rates show there is no detectable difference for either, given the frontal shield trait (Fig. 4). Note that BiSSE is prone to type II errors (Rabosky and Goldberg, 2015) but that we actually cannot reject the null hypothesis here and are hence on safe grounds. Ancestral state reconstruction for the frontal shield states and BiSSE analyses for It has long been known that molecular and morphological approaches (the latter including fossil taxa) must be simultaneously embraced for robust phylogenetic inferences (Pyron, 2015). In this contribution, we have taken a substantial step in contributing new molecular data and a greatly expanded and robustly supported phylogeny for an understudied but ecologically and evolutionarily important phylum. Although we are interested primarily in New Zealand cheilostome bryozoans for reasons stated in our introduction, we have also now

Higher-level cheilostome systematics needs revision
Cheilostome systematics is in a state of flux as molecular studies, coupled with the introduction of genome-skimming, are starting to take off for this diverse clade (Orr et al., 2019a(Orr et al., , b, 2000. In providing a broadly sampled and robustly supported framework to evaluate evolutionary hypotheses we find that less than half of the 29 currently recognized families for which we have multiple genera represented are phylogenetically coherent. Our result emphasizes that much of the current family and higher-level bryozoan systematics, based largely on morphology, is unreliable, and further corroborates previous studies with statistically well-supported, but less broadly sampled, phylogenies (Orr et al., 2019a(Orr et al., , b, 2000. One implication of this observation is that higher-level systematics (involving families, superfamilies and suborders) likely require substantial revision. We have hence refrained from detailing the mismatches of higher-level systematics (Bock, 2020) prematurely, but highlight new evolutionary hypotheses that have emerged, that are potentially supportable by morphological traits, given our molecular inferences (Fig. 2 Notwithstanding some discrepancies between morphology-based hypotheses (Bock, 2020) and molecular data (this study), there is frequently mutual support. Take, for example, the basal grouping of Scrupariidae (Scruparia) as sister taxa to Electridae (Electra), Membraniporidae (Biflustra and Membranipora) and Aeteidae (Aetea) plus Steginoporellidae (Steginoporella) and Calpensiidae (Calpensia) (Fig. 2): all of these families are understood to have diverged prior to the evolution of ovicells (brooding structures) produced from a distal zooid (Ostrovsky, 2013). Our tree now corroborates this hypothesis with full statistical support.
A closely positioned clade formed by Monoporella (Monoporellidae) and Macropora (Macroporidae) shares the presence of large ooecia that evolved from basally articulated spines or costae (Ostrovsky, 2013) (but see next paragraph). The fully supported Arachnopusiidae (Arachnopusia) + Foveolariidae (Foveolaria and Odontionella) relationship is not indicated in present classification schemes (Bock, 2020), as species of Arachnopusia have an ascophoran state, while the Foveolariidae has an anascan state. However, we note that not only is the arachnopusiid frontal shield a straightforward structure to form (unlike other ascophoran structures), but some species in Arachnopusiidae (e.g. A. gigantea) are anascan-like, where the frontal shield is practically non-existent (Hayward, 1995).

A need for even broader taxon sampling to fill gaps
Our ML (Fig. 2, Supplementary Fig. S3) and Bayesian ( Supplementary Fig. S4) trees are largely in agreement with only one clade demonstrating incongruence. This is the fully supported clade comprising Valdemunitella (currently Calloporidae), Figularia (currently Cribrilinidae) and Euthyroides (currently Euthyroididae). Based only on morphology, we might have hypothesized that the Valdemunitella clade (based on 4 species represented by 6 colonies) is closely associated with other representatives of the family Calloporidae (e.g. Crassimarginatella or Callopora), but neither of our trees inferred this position. Rather, our ML tree places this clade (including Figularia and Euthyroides, both currently belonging to other families) in a position close to Monoporellidae and Macroporidae (see paragraph above) and our Bayesian tree places it in a more derived position. However, note that nodes subtending this clade in both trees are poorly supported. Rather than speculating on evolutionary and/or morphological arguments for either or both of these placements, we argue, rather that this indicates that there are many crucial unsampled taxa that would potentially allow a more robust placement of this clade, such as other cribrilinids in addition to Figularia and other calloporids such as Cauloramphus which, similarly to Valdemunitella, has spines encircling the frontal uncalcified membrane, forming a costate shield in some species (Dick et al., 2011). In the event, Valdemunitella, Figularia and Euthyroides are morphologically united, not by a costate shield, but by identical bilobate ooecia with a median suture, and the presence of vicarious avicularia in most of their species.

The evolution of the cheilostome frontal-shield
Historical studies of cheilostome body-wall development and morphology led to the conclusion that ascophoran frontal shields were phylogenetically informative (Banta, 1970;Gordon and Voigt, 1996;Sandberg, 1977). Our results substantiate the observation that characters considered to have deep phylogenetic information such as frontal shields are more evolutionarily labile than previously thought, and sometimes may even be convergent rather than homologous traits (Knight et al., 2011;Orr et al., 2019a). It has already been suggested, for instance, that anascan and ascophoran states, respectively regarded as stemward and crownward, have evolved more than once (Dick et al., 2009;Gordon, 2000;Waeschenbach et al., 2012). We show here that the anascan state is very likely basal in the cheilostome tree and that the change from an anascan to ascophoran state has happened multiple times independently (seven times in Fig. 3), hence likely more times in the history of cheilostome evolution. It is also striking that an ascophoran-state never reverts back to the anascan-state, suggesting that it is evolutionarily unproblematic to evolve a more complex calcified skeleton, but that once this structure is in place, it has not been lost again (Fig. 3). This could be due to genetic or developmental constraints, and/or because the advantages conferred by a calcified frontal shield vastly outweighs its disadvantages. Testing a classic idea that morphological complexity may predict diversification rates (Schopf et al., 1975), we found that the (potentially) more morphologically complex ascophoran-grade cheilostomes have indistinguishable speciation and extinction rates compared with anascan-grade ones.
The frontal shield clearly contains phylogenetic information, but more research is needed to understand when it is informative, and why. As a further example, frontal shields produced by different developmental processes (e.g., lepralioid or umbonuloid (Hayward and Ryland, 1999;Taylor, 2020)) leave such distinct morphological tell-tale signs that it was commonly assumed that members within families constituted only a single type of frontal shield development. Our tree, however, places ascophoran taxa with lepralioid frontal shields (e.g., Powellitheca/Cyclicopora; Celleporina, Galeopsis, Osthimosia) and umbonuloid ones (Exochella; Celleporaria) in the same clades (Fig. 2), as already shown to a lesser extent in earlier extensive studies (Dick et al., 2009;Orr et al., 2019a;Waeschenbach et al., 2012).
Yet, at the more derived part of our inferred tree, the structure of the frontal shield seems to be more phylogenetically informative than seemingly distinct features such as the lyrula (Berning et al., 2014). This is an anvil-shaped tooth-like structure projecting from the orifice that functions in water compensation. Specifically, the clade containing Parasmittina to Hemismittoidea (containing four and three genera of the families Smittinidae and Bitectiporidae respectively) has a non-pseudoporous umbonuloid frontal shield (Gordon, 2000), while the next one containing Schizosmittina to Bitectipora (containing two smittinid and two bitectiporid genera) has a pseudoporous lepralioid shield. The presence of a lyrula seems haphazard among these genera, where those in the Smittinidae have lyrula and those in the Bitectiporidae have a sinus. Our tree suggests new ways of partitioning some of the families and genera of Smittinoidea, which unexpectedly also includes the genera Porella (Bryocryptellidae) and Oshurkovia (Umbonulidae). To summarize, it is clear that a much more thorough and systematic investigation of the development and evolution of frontal shields is necessary for a deeper understanding of ascophoran cheilostomes.

Molecules suggest morphological hypotheses and pinpoint research needs
Another example of traits thought to be phylogenetically related and hence informative is the sinus versus the ascopore, pertaining to the ascophoran plumbing system. Because Microporella, Fenestrulina and Calloporina all have ascopores, they were historically united in the Microporellidae. A previous molecular study has clearly shown that Fenestrulina does not belong in the same clade as Microporella (Orr et al., 2019b). Here, we give molecular support to the hypothesis that Calloporina is not a microporellid and further suggest that Chiastosella (having a sinus, currently belonging to the Escharinidae) and Calloporina (having a slit-like ascopore) belong in the same clade, a relationship supported also by their shared distinctive ovicell (Brown, 1954;Cook et al., 2018). Supporting the long-held hypothesis that an ascopore should evolve by the cutting-off of a sinus, Chiastosella should be basalwards of Calloporina (Cook et al., 2018, p. 218). This is supported by our tree, which also suggests that Chiastosella may be paraphyletic with respect to Calloporina.
In multiple cases, taxa that are considered unique or unusual have placed in phylogenetic positions that suggest hypotheses of their evolutionary relationships based on morphology.
For instance, Rhabdozoum, currently placed in its own family because of its highly distinctive morphology, is basal to Candidae, suggesting that they are closely related and that Candidae sensu stricto may have been derived from a Rhabdozoum-like ancestor. In fact, the initial zooid of the colony (ancestrula) of Rhabdozoum resembles those in some Scrupocellaria species and several of its mature zooidal features such as its ooecia, frontal avicularia and spines are reminiscent of species of Amastigia and Menipea (all Candidae s.s.). Margaretta, another rather distinct genus, is in a family with only one other monospecific genus (Tubucella). Here, Margaretta is inferred to be basal to Catenicellidae, suggesting that Catenicellidae s.s. may have been derived from a Margaretta-like ancestor, although it has always been thought that catenicellids are derived from cribrimorphs (Gordon, 2000;Gordon and Braga, 1994). Much research is required to unravel the mystery of this grouping, given that they are both so distinctive, sharing apparently only rhizoids, rootlets that attach the colony to the substrate. Note that we infer two distinct clades of Catenicellidae, one represented by Catenicella and Cornuticella, which are vittate (frontal pore chambers are long and narrow) and the second including Orthoscuticella and Pterocella, which are foraminate (frontal shield has numerous windows in the gymnocyst). Yet another example is the erect and branching calwelliid Malakosaria whose zooidal features resemble Fenestrulina (Fenestrulinidae), the genus in which Malakosaria nests in our tree.
One taxonomically challenging family deserves special mention. The speciose Celleporidae, with at least 252 described living taxa globally, is mostly characterized by nodular/massive colonies as a result of rapid frontal budding (the building of zooids on top of existing ones).
As a consequence, autozooids are somewhat irregularly disposed and difficult to characterize morphologically. These genera are currently distinguished by the morphology of their ooecia (development of endooecium/tabula) and orifices (always sinuate but the sinus varies from a narrow slit to a broad and shallow concavity). Genus-level hypotheses based on these characters are problematic as indicated by our tree, in which Celleporina, Galeopsis and Osthimosia are non-monophyletic. Buffonellaria is excluded from the family and allied with Buffonellodidae, whereas Celleporaria, historically included in Celleporidae but subsequently split off because of its umbonuloid frontal shield (Harmer, 1957;Cook et al. 2018 p. 182), is reinstated.

Lower-level cheilostome systematics are very robust
Although higher-level systematics are in need of revision, we report that lower-level morphological hypotheses (i.e., species and genera) are very robust, supporting inferences based on common-garden experiments, to put forward the idea that "morphological species" are as good as "genetic species" in cheilostome bryozoans (Jackson and Cheetham, 1990). While Jackson and Cheetham experimented only with a handful of species, we now confirm their hard-earned insight implies that many more species and genera can be treated as distinct evolutionary lineages. This is an important result as many evolutionary and paleontological studies use morphospecies or even morpho-genera as the unit of analyses (Alroy, 2010;Heim et al., 2015). We also note that there are many New Zealand species in our tree that are yet undescribed (c. 20% of those newly sequenced here), indicating that continued exploration in the EEZ of New Zealand is crucial even for such a geographically well-characterized marine clade.

Conclusions
Our work shows that lower-level taxonomic sampling in phylogenetics is vital for understanding higher-level systematics, especially in an understudied group like cheilostome bryozoans. While we have contributed a substantial number of sequences from diverse species, many more must be included for the phylogenetic inferences and reliable systematic groupings for cheilostomes. By contributing molecular data and robustly supported phylogenetic inferences, we have supplied the basis for evolutionary (including phylogenetic) hypotheses that can be further examined. Once we are confident in the topology of at least parts of the cheilostome tree, we can start asking further questions on evolutionary processes.

CRediT authorship contribution statement
LHL conceived the project, RJSO performed the lab, bioinformatic and phylogenetic inference. LHL performed the ancestral state and diversification analyses. MHR and EDM carried out the scanning electron microscopy, EDM and DPG identified the bryozoans, and AMS led the major ship-based collecting expeditions in which HLM and DPG took part in.
LHL, EDM and RJSO drafted the paper. All coauthors contributed to editing and revisions.

Acknowledgements
We thank Seabourne Rust for helping to sort and identify samples from the University of