The joint evolution of the Myxozoa and their alternate hosts: A cnidarian recipe for success and vast biodiversity

The relationships between parasites and their hosts are intimate, dynamic and complex; the evolution of one is inevitably linked to the other. Despite multiple origins of parasitism in the Cnidaria, only parasites belonging to the Myxozoa are characterized by a complex life cycle, alternating between fish and invertebrate hosts, as well as by high species diversity. This inspired us to examine the history of adaptive radiations in myxozoans and their hosts by determining the degree of congruence between their phylogenies and by timing the emergence of myxozoan lineages in relation to their hosts. Recent genomic analyses suggested a common origin of Polypodium hydriforme, a cnidarian parasite of acipenseriform fishes, and the Myxozoa, and proposed fish as original hosts for both sister lineages. We demonstrate that the Myxozoa emerged long before fish populated Earth and that phylogenetic congruence with their invertebrate hosts is evident down to the most basal branches of the tree, indicating bryozoans and annelids as original hosts and challenging previous evolutionary hypotheses. We provide evidence that, following invertebrate invasion, fish hosts were acquired multiple times, leading to parallel cospeciation patterns in all major phylogenetic lineages. We identify the acquisition of vertebrate hosts that facilitate alternative transmission and dispersion strategies as reason for the distinct success of the Myxozoa, and identify massive host specification‐linked parasite diversification events. The results of this study transform our understanding of the origins and evolution of parasitism in the most basal metazoan parasites known.

1999; Yoder & Nuismer, 2010) and scenarios such as the acquisition of new hosts or new niches as well as within-host competition foster parasite speciation (Morand, Krasnov, & Littlewood, 2015).
Despite the central importance of species interactions to the diversification of life, the knowledge about the processes by which species undergo reciprocal evolutionary change through natural selection is patchy (Carmona, Fitzpatrick, & Johnson, 2015). An explicit framework to understand the role of species diversification in coevolution is still required, and authors repeatedly caution against a direct link of the two (Althoff, Segraves, & Johnson, 2014;Brockhurst & Koskella, 2013;Poisont, 2015). However, investigating the extent of codivergence or parallel speciation of parasites and their hosts under strict analytical conditions (timing, sympatry, etc.) allows for insights into the joint evolutionary history of two organisms. This provides a basis for studies uncovering underlying processes that govern the nature of their interactions (Charleston & Perkins, 2006) and can help us to better understand the dynamics and history of biodiversity.
Myxozoans are a group of diverse cnidarians, representing almost one-fifth of presently known Cnidarian species (Zhang, 2011). Extremely reduced as an adaption to parasitism, myxozoans are composed of only a few cell types throughout their life cycle. Characteristic spores (Figure 1, centre) serve as transmission stages between invertebrate and vertebrate hosts, predominantly in aquatic habitats.
Myxozoans are strongly opportunistic and show a high degree of flexibility with regard to exploring hosts of different animal groups.
Annelids and bryozoans are known as definitive hosts with 55 life cycles elucidated to date (Supplementary file S1). Fish compose the largest number of known intermediate hosts, but amphibians, reptiles, birds and mammals are also exploited (Fiala, Barto sov a-Sojkov a, Okamura, & Hartikainen, 2015; Fiala, Barto sov a-Sojkov a, & Whipps, 2015). Some fish-parasitic myxozoans even evolved into hyperparasites of platyhelminth endo-and ectoparasites of fishes (Freeman & Shinn, 2011;Overstreet, 1976;Siau, Gasc, & Maillard, 1981), and a muscle-dwelling species was found in an octopus (Yokoyama & Masuda, 2001). The history of these host-parasite associations is largely unexplored, but their interactions likely represent important drivers of the evolution of parasitism at the base of the Metazoa.
Based on a limited number of life cycle discoveries, in 2007, it was first suggested that the invertebrate host type mirrored large-scale myxozoan phylogeny as well as 18S rRNA secondary structure (Holzer, Wootten, & Sommerville, 2007). Most recently, Kod adkov a, Barto sov a-Sojkov a, Holzer, and Fiala (2015) determined that cartilaginous fish represent ancestral states for a number of phylogenetic lineages and these basal parasite lineages and their cartilaginous hosts likely co-originated in the Silurian. These results prompted us to document the common evolutionary history of myxozoans and their vertebrate and invertebrate hosts by cophylogenetic analyses, study the patterns of cospeciation and diversification in different host groups and investigate temporal congruence of these evolutionary events, using molecular clock analysis. Cospeciation has rarely been observed in fish parasites (summarized in Vanhove et al., 2015), and studies analysing codivergence patterns of heteroxenous parasites in more than one of their host groups are still missing. Our analyses cover all myxozoan host-parasite associations for which sequence data are presently available and decipher the series of events that led to the distinct success of this enigmatic parasite group at the base of the Metazoa.

| Sequence data and phylogenetic analyses
A prerequisite for cophylogenetic studies is the reconstruction of reliable phylogenetic trees at both family and species levels. Largescale relationships among myxozoan taxa are relatively well established based on 18S rDNA sequences, and topologies are congruent with those based on the limited data available from other genes (recently reviewed by Fiala, Barto sov a-Sojkov a, Okamura, et al. (2015) and Fiala, Barto sov a-Sojkov a, and Whipps (2015)). We used a comprehensive collection of 18S rDNA sequences from 633 taxa (1,638 bp; Supplementary file S2), predominantly from vertebrate (563 taxa) and to a much lesser extent from invertebrate (110 taxa) hosts (with 40 species overlapping between the two data sets). For construction of invertebrate phylogenetic trees, 18S rDNA sequence data were used (23 taxa, 1,859 bp) while for vertebrates, mitochondrial 16S rRNA gene sequences were compiled (245 taxa, 1,572 bp), as the latter is the most comprehensively sampled gene in fishes.
For vertebrate hosts, a second data set was compiled focusing on a considerably larger data set (full mitogenomes, 105 taxa, 15,913 bp) and a higher taxonomic level, that is, host families. This data set allows for analyses of the early history of cophylogeny signatures in vertebrate hosts, independent from recent events.
Nucleotide sequences were aligned applying the MAFFT version 7.017 alignment (Katoh, Misawa, Kuma, & Miyata, 2002) implemented in GENEIOUS version 9 (Kearse et al., 2012), using the E-INS-i algorithm, with a gap opening penalty (-op) 1.5 (myxozoans)/1.5-4.0 (vertebrates and invertebrates) and gap extension penalty (-ep) 0.0. The alignments were edited, and highly variable sections were removed. Maximum parsimony (MP) analyses were performed in PAUP* version 4.b10 (Swofford, 2003), using a heuristic search with random taxa addition, the ACCTRAN option, TBR swapping algorithm, all characters treated as unordered and gaps treated as missing data. JModelTest (Posada, 2008) was used to select the best-fitting model of evolution, using the corrected Akaike information criterion. Maximum-likelihood (ML) heuristic searches were performed in RAXML version 7.0.3 (Stamatakis, 2006), using the GTR + Γ model of nucleotide substitution. Bootstrap searches included 500 replicates for both MP and ML analyses. Bayesian inference (BI) analyses were performed in MRBAYES version 3.0 (Ronquist & Huelsenbeck, 2003), using the GTR + Γ model of evolution. Posterior probabilities were estimated from 1 million generations with four simultaneous MCMC chains, sampled at intervals of 100 trees, with burn-in set to 10%. To ensure convergence and an effective sample size, results were verified with TRACER version 1.6 (Rambaut, Suchard, Xie, & Drummond, 2014). Problems associated with single-gene phylogenies of rDNA and mtDNA sequences (Jeffroy, Brinkmann, Delsuc, & Philippe, 2006;Rokas & Carroll, 2006) were disrupted by comparison of the resultant host phylogenetic trees with those of other studies employing large molecular data sets (Andrade et al., 2015;Betancur-R et al., 2017;Near et al., 2012;Nesnidal et al., 2013;Struck et al., 2011) and by exclusion of taxa with a strongly divergent position from these studies, for example, if a fish species did not cluster with other members of the same host family.

| Diversification estimates
Overall rates of diversification for myxozoans were estimated in R version 3.2.4 (R Core Team, 2013), using the APE package version 3.4 (Paradis, Claude, & Strimmer, 2004) and the GEIGER package version 2.0.6 (Pennell et al., 2014). Diversification estimates were calculated based on a pure birth process and on a speciation:extinction rate of 2:1. The constant rate test (Pybus & Harvey, 2000) was used to test for even rates of cladogenesis through time. The Monte Carlo CR (MCCR) test was used to account for incomplete taxon sampling. A lineages through time (LTT) plot was produced, and the relative cladogenesis (RC) statistic (Pennell et al., 2014) was used to detect unusually rapid diversification shifts leading to lineage-rich clades in the tree.

| Cophylogeny analyses
Two methodological approaches were used to compare host and parasite phylogenetic relationships: (i) event-based tree reconciliation and (ii) global Fit analyses. The reconciliation of tree topologies interprets co-evolution as a stochastic process with cospeciation (concomitant host and parasite speciation), duplication (parasite speciation on one host lineage), host switching (colonization of a new host) and lineage sorting (disappearance of a parasite lineage from a host) as discrete events. Partially resolved ML trees with polytomies for clades with bootstrap support <50% were used for cophylogeny analyses. In some phylogenetic trees where polytomies affected >35% of all taxa, fully resolved ML trees were used (indicated in results in Table 1). To eliminate the risk of maximizing cospeciation events by tree reconciliation (de Vienne et al., 2013), we performed a parameter-adaptive approach in CORE-PA version 0.5 (Merkle, Middendorf, & Wieseke, 2010), testing 10 4 cost sets on the quality function, using the simplex method and avoiding an a priori cost assignment. Global fit estimates were implemented in PARAFIT (Legendre, Desdevises, & Bazin, 2002; APE package version 3.4 in R), which uses a permutation test on a matrix of raw patristic distances to evaluate codivergence and thereby overcomes the need for well-resolved tree topologies. Distance-based cophylogenetic analyses are hence considered less biased than topology-based methods. Multiparasite host species as well as multihost parasite taxa are anomalous under a codivergence paradigm (Banks & Paterson, 2005). Most myxozoans are highly host specific with regard to their vertebrate hosts (Moln ar & Eszterbauer, 2015). Only a few marine histozoic species, that is, Kudoa spp. (Burger & Adlard, 2011;Gleeson, Bennett, & Adlard, 2010;Whipps & Kent, 2006) and Enteromyxum leei (Sitj a-Bobadilla & Palenzuela, 2012), are able to infect a broad range of hosts from different fish families. Myxozoans appear to be less host specific with regard to their invertebrate hosts, with, for example, Tubifex tubifex as a known host for >30 different taxa. It is unlikely that such host-parasite associations are the result of true cophylogeny but rather, for example, host switching to a more abundant alternative host (Lootvoet, Blanchet, Gevrey, Buisson, & Tudesque, 2013;Poulin, 2011). As the mechanisms leading to the generation of such associations are unknown, independent from the data set, species with >8 hosts or parasites were excluded from analyses. This affected only a small percentage of the data, namely 0.8% of myxozoan species and 2.9% of hosts (detailed in Supplementary file S3).
All myxozoans and their invertebrate hosts were analysed together. Due to the large number of taxa sequenced from their vertebrate hosts and the respective size of trees, the data set was split into subsets, representing the major clades in the parasite tree (Figure 1) but excluding malacosporeans, as uncertainty exists with regard to "true" fish hosts in this group, that is, hosts that can form spores, which are infective to an invertebrate host (Patra, Hartigan, Morris, Kod adkov a, & Holzer, 2017). Analyses were also performed after reduction in the large subclades "Biliary tract I," "Histozoic I" and "Histozoic III'; Figure 1) and of these subclades themselves, allowing assessment of recent and ongoing cospeciation in closely related myxozoans and their hosts. Results from diversification and cophylogeny analyses were considered statistically significant at a = 0.05.

| Molecular clock analyses
To provide meaningful information on the timing and divergence of the major clades of myxozoans in relation to that of their hosts and other metazoans, concatenated alignments of six protein-coding genes, that is, aldolase (200 aa), triosephosphate isomerase (217 aa), phosphofructokinase (175 aa), methionine adenosyltransferase (348 aa), elongation factor 1 alpha (418 aa) and ATP synthase beta chain (430 aa) were used. Using an existing data set (Erwin et al., 2011), we mined published and new genomic and transcriptomic parasite and host data sets resulting in the final alignment of 2,444 aa, including 10 myxozoans + 128 other metazoan taxa (Supplementary file S4). Divergence times were estimated using BEAST version 2.4.7 (Bouckaert et al., 2014;Drummond, Suchard, Xie, & Rambaut, 2012).
The ML tree topology computed using RAXML under the WAG+ Γ model of evolution was used as a user-specified starting tree. The best model of evolution for this analysis was selected by PROTTEST 3.4.2 (Darriba, Taboada, Doallo, & Posada, 2011), using the Akaike information criterion. The same model was used for molecular clock analysis. The BEAST input file was constructed using BEAUTI (within BEAST package). The lognormal relaxed molecular clock model, which accounts for independent rates of heterogeneity from lineage to lineage (Thorne & Kishino, 2002), and the Yule speciation prior set were used to calculate divergence times and corresponding credibility intervals. Fourteen fossil calibration points (Supplementary file S5; adapted from Dohrmann & W€ orheide, 2017) were used. MCMC analyses were run for 10,000,000 generations, and after 10% burnin sampled every 1,000 generations. The convergence of chains was confirmed using TRACER version 1.6 (Rambaut et al., 2014) and chronograms generated in TREEANNOTATOR version 2.4.6 (within BEAST package). Polypodium hydriforme is a sister taxon of the Myxozoa (Chang et al., 2015;Okamura, Gruhl, & Reft, 2015) with high nucleotide substitution rate and assumed long-branch attraction (LBA) to the Myxozoa (Evans, Lindner, Raikova, Collins, & Cartwright, 2008).
To estimate reciprocally independent divergence times of the Myxozoa and P. hydriforme and ascertain the influence of these taxa on the calculated origin of the Cnidaria, we performed BEAST analyses excluding either or both of the taxa from the data set. To confirm whether phylogenetic placement of the two long-branching taxa as sister lineages is correct, the position of each long-branching taxon was examined separately within a large metazoan phylogenomic data set (51,940 aa; Chang et al., 2015; data re-analysed).

| Myxozoa
In accordance with published records on smaller data sets (Barto sov a, Fiala, & Hyp sa, 2009;Fiala, 2006 stricto (s. str.; referred to as sphaerosporids) and two large clades, previously termed "freshwater" and "marine" clades (Fiala, 2006). All four lineages are highly supported (bootstrap values 95%-100%) in both MP and ML analyses and maximum posterior probabilities in BI (not shown). Mapping of host and habitat characteristics to the parasite tree ( Figure 1) provided new insights into the evolution and classification of the Myxozoa. While the taxa from the "marine" clade rarely inhabit freshwater environments (<2%), almost a quarter of the species clustering in the "freshwater" clade were collected in marine hosts and habitats. Twelve years after the original definition of these clades by Fiala (2006), we show that they are less defined by their environment (freshwater vs. marine) but rather by their invertebrate host type, oligochaete-infecting myxozoans (OIM) vs.

| Myxozoa-vertebrate hosts
Significant cophylogenetic signal was detected between myxozoans and their vertebrate hosts, in all major clades; however, the outcome of the analyses varied depending on the data set (Table 1). For all data sets but sphaerosporids, tree reconciliation methods ascribed a higher number of cospeciation events to higher-level taxa (host families and above) while PARAFIT detected more associations in lower level taxa (host species; Table 1) The PIM data set (161 parasites) contains two extremely diverse subclades of closely related taxa infecting the biliary tract and the muscle (biliary tract I and histozoic I in Figure 1). Their inclusion resulted in significant global fit (p = .001), but tree reconciliation by CORE-PA rejected a cophylogeny scenario of myxozoans and their fish hosts unless these two subclades were collapsed, resulting in a database of 76 vs. 161 parasite taxa for the PIM clade (Table 1, Supplementary file S6_Figure 2, PIM reduced data set). CORE-PA is sensitive to analyses of large phylogenetic trees where host switches become more and more likely while, simultaneously, cospeciation events become less likely (Merkle et al., 2010). Independent analysis of the two species-rich subclades resulted in cophylogenetic signal only by tree reconciliation but no significant outcome by distancebased methods ( #P = number of parasite taxa, #H = number of host taxa used for analysis (in case data sets differ due to malpositioning of host taxa in phylogenetic trees these are stated "tree reconciliation/Global fit"); #C/D/HS/S = number of cospeciation/duplication/host switching/sorting events, cost and quality (CORE-PA) and ParaFitGlobal value as well as significance (yes/no for CORE-PA, probability for PARAFIT). Reconciliations of myxozoans and invertebrate hosts based on 18S rDNA sequences, of vertebrate hosts based on 16S rRNA (vertebrate host species) or full mitogenomes (vertebrate host families). Topology-based cophylogenetic analyses used ML trees with nodes supported by a bootstrap value under 50 collapsed and fully resolved ML trees if >35% of taxa affected by polytomy (marked † ). ‡ second-best solution in CORE-PA chosen as best result due to unrealistically high costs for certain events, which are not biologically plausible (Desdevises, Morand, Jousson, & Legendre, 2002), clade names as in Figure 1 with PIM = polychaete-infecting myxozoans and OIM = oligochaete-infecting myxozoans, § excludes data from cyprinids. Significant outcomes are listed in bold.
analysis of phylogenetic congruence on host family level (mitochondrial genome data) produced a significant outcome in both methods (  To investigate the effect of the high rate of nucleotide heterogeneity on the position and divergence time estimate for Polypodium hydriforme and the Myxozoa and on potential long-branch attraction (LBA) phenomena, we performed independent analyses. The analysis of our and Chang's data set (Chang et al., 2015) proved that, exclud-  to have a faster rate of molecular evolution in order to win the "arms race" against their hosts (e.g., Bromham, Cowman, & Lanfear, 2013;Paterson et al., 2010). The rate heterogeneity across genes within the Myxozoa can be as high as that between myxozoans and other organisms (Hartigan et al., 2016;this study), and mitochondrial gene order and organization is highly variable (Takeuchi et al., 2015;Yahalomi et al., 2017), indicating a considerably accelerated rate of molecular evolution, possibly the fastest known among eukaryotes.

| Diversification rate and lineages through time
This may well be explained by the extraordinary level of radiation that occurred within this group (Castro, Austin, & Dowton, 2002;Eo & DeWoody, 2010

Independent from exact dates, Cnidaria originated in the Early-to
Mid-Cryogenian, suggesting a "slow burn" rather than a "Cambrian explosion" scenario for the evolution of the Cnidaria, a view supported by palaeontological (Budd, 2008;Fedonkin, 2003;Seilacher, Bose, & Pfluger, 1998;Van Iten et al., 2013) and molecular data (Cartwright & Collins, 2007;Peterson, Cotton, Gehling, & Pisani, 2008  and cell-in-cell development occurs, at least in initial stages. However, the binucleate cell in P. hydriforme infects fish (Raikova, 1994), while in myxozoans, it infects annelids (summarized in Morris, 2010) in contrast to first stages in fish that are multicellular (El-Matbouli, Hoffmann, & Mandok, 1995;Morris & Adams, 2008). The binucleate cell in P. hydriforme is a larval, diploid stage (Raikova, 1994), while in myxozoans, it is represented by a haploid binucleate cell that merges to form a zygote in the annelid host (Morris, 2012). Their homology is hence unlikely. Cell-in-cell development has been interpreted as another common feature of early P. hydriforme and most of myxozoan development, but it also occurs convergently in parasitic protists of marine invertebrates, the Paramyxida (Ward et al., 2016).
Recent genomic data further support an independent origin of Myxozoa and P. hydriforme: the genome size of the myxozoan Kudoa iwatai amounts to only 4% of that of P. hydriforme (22.5 Mb vs. 561 Mb), and rigorous reductions depleted genes related to development, cell differentiation and cell-cell communication, only in the Myxozoa (Chang et al., 2015). The overlap of exclusive orthologous groups of genes (OGs) between K. iwatai and P. hydriforme is lower than between K. iwatai and free-living Nematostella or Hydra (Chang et al., 2015), indicating a somewhat closer relationship to these freeliving taxa than to the parasitic sister taxon. Finally, P. hydriforme is a monotypic species that appears to be an evolutionary blind lineage, partially parasitic in old aciperseriform fishes, a living fossil that failed to evolve and radiate into modern fish lineages, while the Myxozoa, parasitic throughout their life cycle, successfully established in different invertebrate hosts, and thrived and diversified after multiple entries into fishes. The analysis of presently available developmental and genome-based features in the light of a dated origin and hostparasite phylogenetic convergence strongly suggests that the Myxozoa and P. hydriforme represent not one but two independent routes to endoparasitism in the Cnidaria and that the "Endocnidozoa" the Myxozoa as sister groups to the Medusozoa, it is likely that P. hydriforme and the Myxozoa represent independent lineages drawn together by long-branch attraction of their highly divergent genes and genomes. The lack of other, closely related taxa in combination with fossils (Wiens, 2005) presently prevents resolving this relationship.
Within the Myxozoa, phylogenetic clustering in the tree based on six molecular genes was in accordance with 18S rDNA tree reconstructions, with sphaerosporids consistently positioned basal to the PIM and OIM lineages, while, in the past, they had sometimes been reconciliated as sister to the PIM clade (Barto sov a et al., 2013;Karlsbakk & Køie, 2009). Basal to the known annelid hosts of myxozoans two groups can be found, Haplodrili (= Archiannelida, five families) and Sipuncula (Andrade et al., 2015;Dunn et al., 2008;Struck et al., 2011). Based on congruent myxozoan and invertebrate host trees, these offer themselves as ancestral hosts of sphaerosporids, while modern taxa may well have adapted more recent annelid lineages. The most basal sphaerosporid clade (Barto sov a et al., 2013) is represented by isolates exclusively from fishes in marine habitats, which also home these evolutionary old annelids (Struck et al., 2007). According to the present reconstruction of the history of myxozoan host acquisitions, sipunculids are excellent candidates for invertebrate hosts of basal sphaerosporids and are worthy of indepth study including 18S rDNA screening, especially because infected specimens were found to harbour more than one spore phenotype (Ikeda, 1912). In contrast to previous reports suggesting a different host for sphaerosporids (Barto sov a et al., 2013;Holzer et al., 2007), based on the present analysis, we are able to pinpoint specific taxa.
Findings of myxozoans in fish-parasitic flatworms (Freeman & Shinn, 2011;Overstreet, 1976;Siau et al., 1981) and a free-living mollusc (Yokoyama & Masuda, 2001) describe spores exhibiting a bilateral symmetry, in contrast to the triradial symmetry of spores from annelid hosts, and represent generic morphotypes known from fishes (Kudoa, Fabespora and Myxidium). Further supported by their close phylogenetic relationship to histozoic taxa from marine fishes (Freeman & Shinn, 2011), it is likely that they were acquired from fish hosts by blood and tissue feeding helminths and settled as hyperparasites within them.  (Canning, Okamura, & Curry, 1996) provide exciting perspectives regarding additional host discoveries that may further improve our understanding of the evolutionary history of the Myxozoa.
Another actual debate is that of a marine vs. freshwater origin of myxozoans. Cnidarians live predominantly in marine environments, and Kent et al. (2001) suggested that myxozoans first became endoparasitic in old marine annelid worms. However, this interpretation predates the rediscovery of malacosporeans, the oldest clade of myxozoans, in bryozoan hosts. Malacosporeans infections have only been detected in Phylactolaemata, the radix group of bryozoans that occurs exclusively in freshwater habitats (reviewed in Taylor & Waeschenbach, 2015). This, together with the fact that P. hydriforme occurs in freshwater fish, prompted some authors to suggest a freshwater origin of the Myxozoa. However, the last common ancestor of today's Phylactolaemata first evolved in marine environments and only secondarily occupied freshwaters habitats (Koletic, Novosel, Rajevic, & Franjevic, 2015).   (Tedersoo, Bahram, & Dickie, 2014) explain myxozoan hyperspeciation in some clades, and for a better understanding of the ecological and evolutionary principles than impact on diversity in this parasite group much research is still required. However, the evolutionary history of their hosts is presently an essential predictor of myxozoan diversity. Thus, our study clearly highlights the importance of host phylogenetic information for explaining parasite speciation. From an evolutionary perspective and based on the early birth age of the Myxozoa, it is possible that myxozoan parasites themselves contributed to promoting host diversification. Parasites have often been considered a threat to biodiversity because of their negative impact on host population persistence (de Castro & Bolker, 2005;Valenzuela-Sanchez et al., 2017), but they may actually play important roles in maintaining and promoting host biodiversity (Buckling & Rainey, 2002;Karvonen & Seehausen, 2012). In myxozoans, the oldest metazoan parasites on Earth, the evolutionary history of host interactions is key for their distinct success, leading to an extraordinary richness of modern taxa in all aquatic habitats and fish lineages, whose contribution to total levels and patterns of aquatic biodiversity on Earth is expected to be considerable.

| CONCLUSIONS
Our data suggest an origin of the Myxozoa in extinct marine Phylactolaemata and archiannelids, around 651 Ma. Myxozoans first evolved in invertebrates and their main lineages split before a two-host life cycle was acquired, with the members of three major clades of myxozoans strictly using specific invertebrate groups (Bryozoa, Polychaeta, Oligochaeta). Fish were conquered as second hosts multiple times, firstly when prehistoric cartilaginous fish became available (approx. 450 Ma). Chimaeras, sharks and rays still serve as hosts for the oldest lineages of the two major annelid-infecting clades. Myxozoans clearly show a common evolution with their vertebrate hosts, with lineages over time in (i) cartilaginous fishes, (ii) tetrapods and (iii) following the emergence pattern of modern ray-finned fishes. Due to their origin in invertebrates prior to the emergence of vertebrates, invertebrate and myxozoan phylogenies are highly congruent, while the co-evolutionary signatures of myxozoans and fish depend on the tested data set.
Despite a sister lineage relationship of Polypodium hydriforme and the Myxozoa, the origin of P. hydriforme in old acipenseriform fishes likely represents an independent conquest of fish as cnidarian hosts. We point out important differences in the development of the two parasite lineages and demonstrate that characteristics considered developmental homologies represent divergences. Myxozoans are parasitic throughout their life cycle and, in contrast to P. hydriforme thrived and diversified massively, especially after the acquisition of fish as secondary hosts, to whose diversification they may well have contributed, with massive bursts of speciation occurring only in highly diverse host lineages.

ACKNOWLEDG EMENTS
The study was supported by the European Commission's Research

CONFLI CT OF INTEREST
The authors declare no competing interests.