Life in a rock pool: Radiation and population genetics of myxozoan parasites in hosts inhabiting restricted spaces

Introduction Intertidal rock pools where fish and invertebrates are in constant close contact due to limited space and water level fluctuations represent ideal conditions to promote life cycles in parasites using these two alternate hosts and to study speciation processes that could contribute to understanding the roles of parasitic species in such ecosystems. Material and methods Gall bladder and liver samples from five clinid fish species (Blenniiformes: Clinidae) were morphologically and molecularly examined to determine the diversity, prevalence, distribution and host specificity of Ceratomyxa parasites (Cnidaria: Myxozoa) in intertidal habitats along the coast of South Africa. Phylogenetic relationships of clinid ceratomyxids based on the SSU rDNA, LSU rDNA and ITS regions were assessed additionally to the investigation of population genetic structure of Ceratomyxa cottoidii and subsequent comparison with the data known from type fish host Clinus cottoides. Results and discussion Seven Ceratomyxa species including previously described Ceratomyxa dehoopi and C. cottoidii were recognized in clinids. They represent a diverse group of rapidly evolving, closely related species with a remarkably high prevalence in their hosts, little host specificity and frequent concurrent infections, most probably as a result of parasite radiation after multiple speciation events triggered by limited host dispersal within restricted spaces. C. cottoidii represents the most common clinid parasite with a population structure characterized by young expanding populations in the south west and south east coast and by older populations in equilibrium on the west coast of its distribution. Parasite and fish host population structures show overlapping patterns and are very likely affected by similar oceanographic barriers possibly due to reduced host dispersal enhancing parasite community differentiation. While fish host specificity had little impact on parasite population structure, the habitat preference of the alternate invertebrate host as well as tidal water exchange may be additional crucial variables affecting the dispersal and associated population structure of C. cottoidii.


Introduction
Intertidal rock pools where fish and invertebrates are in constant close contact due to limited space and water level fluctuations represent ideal conditions to promote life cycles in parasites using these two alternate hosts and to study speciation processes that could contribute to understanding the roles of parasitic species in such ecosystems.

Material and methods
Gall bladder and liver samples from five clinid fish species (Blenniiformes: Clinidae) were morphologically and molecularly examined to determine the diversity, prevalence, distribution and host specificity of Ceratomyxa parasites (Cnidaria: Myxozoa) in intertidal habitats along the coast of South Africa. Phylogenetic relationships of clinid ceratomyxids based on the SSU rDNA, LSU rDNA and ITS regions were assessed additionally to the investigation of population genetic structure of Ceratomyxa cottoidii and subsequent comparison with the data known from type fish host Clinus cottoides.

Results and discussion
Seven Ceratomyxa species including previously described Ceratomyxa dehoopi and C. cottoidii were recognized in clinids. They represent a diverse group of rapidly evolving, closely related species with a remarkably high prevalence in their hosts, little host specificity and frequent concurrent infections, most probably as a result of parasite radiation after multiple speciation events triggered by limited host dispersal within restricted spaces. C. cottoidii represents the most common clinid parasite with a population structure characterized by young expanding populations in the south west and south east coast and by older populations in equilibrium on the west coast of its distribution. Parasite and fish host population structures show overlapping patterns and are very likely affected by similar oceanographic barriers possibly due to reduced host dispersal enhancing parasite community differentiation. While fish host specificity had little PLOS  Introduction 31,38] display a high degree of morphological spore plasticity. This raises questions whether the observed variation is determined by deformations of their presumably thin-walled shell valves, a feature previously documented for some ceratomyxids [39], or alternatively is due to the concurrent infections of several Ceratomyxa species of different morphology, as mixed myxozoan infections [17], including those of Ceratomyxa spp. [40], are a common phenomenon.
Myxozoans are generally recognized as parasites with an accelerated rate of molecular evolution [41,42] explained by an extraordinary level of radiation [43,44]. Such rapid molecular evolution and associated long-branch attraction (LBA) can greatly obscure the results of phylogenetic studies of the Myxozoa [41]. Substitution rates significantly differ among myxozoan lineages. For example, the Ceratomyxa clade, placed within the marine myxosporean lineage [45,46] includes several long-branching taxa sometimes clustering together, most probably due to LBA [46]. Ceratomyxa spp. from the same fish host family (and sometimes from the same geographic location) tend to cluster together, however, this is not always the rule [27,29].
Though population structure of both, host and parasite, is most relevant to co-evolutionary processes [47][48][49] and myxozoan evolutionary history has been significantly linked to hostparasite co-evolution [50,51], unfortunately, to date, investigations determining the population structure of myxozoan parasites in relation to that of their hosts are still missing. For this purpose, we selected the host-parasite model system of Clinus cottoides, rendering the detailed knowledge of its biology, distributional range, and population structure [10][11][12], and Ceratomyxa cottoidii for which information on its fish host species and distribution was collected during the previous [20] and this studies. Moreover, an ecological system with geographically limited host dispersion as seen in C. cottoides is ideal for study.
Our aims were to i) examine diversity, prevalence, distribution and host species spectra of ceratomyxids from South African clinids, ii) reconstruct the phylogenetic relationships among clinid ceratomyxids and to other members of the Ceratomyxa clade, and iii) reveal the patterns of population structuring of C. cottoidii and to investigate if they are linked to geographical or host-driven isolation of parasite populations. Based on previous knowledge, we hypothesized that ceratomyxids species in South African clinids are strictly host-specific, closely related and occur in mixed infections. Finally, we hypothesized that C. cottoidii as a specialist parasite follows the population structure of its fish host C. cottoides.

Sample collection, study area and time schedule
In total, 143 fish specimens of the family Clinidae (Perciformes) belonging to

Processing of samples
Fish livers and gall bladders (in total 152 samples; S1 Table) were examined directly in the field or fixed in 70% ethanol and examined subsequently in the laboratory for the presence of myxosporean infections by light microscopy (LM) on Leica DM750 and Olympus BX51 microscopes. Plasmodia and myxospore morphology were documented with an Olympus DP70 digital camera using differential interference contrast. Samples were stored in 400 μL of TNES urea buffer (10 mM Tris-HCl with pH 8, 125 mM NaCl, 10 mM EDTA, 0.5% SDS and 4 M urea) or in 90% ethanol for subsequent DNA extraction.

DNA extraction, PCR amplification, cloning and sequencing
Total DNA was extracted from all collected samples using a standard phenol-chloroform protocol, after an overnight digestion with proteinase K (50 μg mL -1 ; Serva, Germany), at 55˚C. DNA was re-suspended in 50-100 μL −1 of DNAse-free water and left to dissolve overnight at 4˚C.
Different regions of molecular markers within the nuclear ribosomal operon (SSU rDNA, LSU rDNA, ITS region) were amplified in this study with primer combinations as listed in Table 1. At the beginning we aimed at obtaining the sequence data of all three markers for Ceratomyxa parasites using the general Ceratomyxa (SSU rDNA) or general clinid Ceratomyxa (LSU rDNA, longer ITS region) primers from a random selection of several samples to compare the phylogenetic clustering patterns. Later on, we selected the longer ITS region as the most suitable marker for species detection based on the highest amount of data informativeness and best tree resolution and continued in PCRs of additional samples only for this region. Using these approaches, seven Ceratomyxa spp. were uncovered (C. cottoidii, C. dehoopi, further undescribed species Ceratomyxa sp. 1, Ceratomyxa sp. 2, Ceratomyxa sp. 3, Ceratomyxa sp. 4 and Ceratomyxa sp. 5). However, amplification of all three loci from samples chosen for amplification of all markers was not successful in all cases. This situation was further complicated by identification of concurrent infections of multiple Ceratomyxa spp. from clones placed in different clinid Ceratomyxa subclades (representing species) originating from single PCR products. To overcome these obstacles, which would bias the exact assessment of diversity and prevalence of Ceratomyxa spp., we designed seven sets of species-(lineage) specific primers targeting the shorter ITS region of each Ceratomyxa species (Table 1). This region was nested within the previously amplified longer ITS region and provides enough interspecific sequence variation and sufficient number of sequences for primer design. We subsequently PCR-screened all 152 samples.
A table with detailed information on sample localities and fish species, DNA extracted, numbers of PCR products and clones sequenced for each locus of a particular Ceratomyxa species is supplemented (S1 Table).

Phylogenetic analyses
Maximum parsimony (MP) analysis was performed in PAUP Ã v4.b10 [56], 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. Maximum likelihood (ML) analysis was performed in RAxML v7Á0Á3 [57] using the GTR + Γ model of nucleotide substitution. Bootstraps were based on 1000 replicates for both analyses. Bayesian inference (BI) analysis was performed in MrBayes v3.0 [58], using the GTR + Γ + I model of evolution. Analyses were initiated with random starting trees, four simultaneous MCMC chains sampled at intervals of 100 trees and posterior probabilities estimated from 1 million generations for the Ceratomyxa-SSU rDNA dataset (burn-in 100,000), 3 million generations for the clinid-Ceratomyxa-SSU rDNA dataset (burn-in 600,000), 3.5 million generations for the clinid-Ceratomyxa-LSU rDNA dataset (burn-in 350,000) and the clinid-Ceratomyxa-ITS dataset (burn-in 710,000). Suitable burn-in levels were chosen in Tracer v1.4.1 [59].

Distances
The intraspecific and interspecific divergences were determined from proportional distances (in %) which were calculated in Geneious from the datasets previously used for the phylogenetic analyses. These datasets were additionally adjusted by excluding the very short sequences and by trimming the 5' and 3' ends of the remaining sequences to achieve their same length. The intraspecific cut-offs for the studied markers in the group of long branching ceratomyxids from clinid hosts were as follows: SSU rDNA-7.5%, LSU rDNA-10%, ITS region-20%.

Haplotype networks, AMOVA and population genetic statistics
The ITS region (Ceratomyxa cottoidii-ITS haplotype dataset, 76 sequences) was selected for haplotype networks, analysis of molecular variance (AMOVA) and population genetic statistics, as it provides enough sequence variability to distinguish the population structure patterns in myxozoans [60].
To examine the evolutionary relationships among haplotypes in populations of C. cottoidii from different localities and clinid fish species statistical parsimony networks (TCS) based on pairwise differences were constructed using PopART v1.7 [61]. Sampling localities were set as traits for the haplotype network which grouped the sequences (number of sequences = n) with regard to geography (Cape Columbine n = 6, Jacobs Bay n = 2, Granger Bay n = 1, Mouille Point n = 4, Kommetjie n = 1, Kalk Bay n = 24, False Bay n = 1, De Hoop n = 14, Jongensfontein n = 2, Herolds Bay n = 19, Jeffreys Bay n = 2). Fish host species were set as traits for the haplotype network which grouped the sequences with regard to host specificity (C. cottoides n = 61, C. superciliosus n = 9, C. acuminatus n = 1, C. brevicristatus n = 1, M. dorsalis n = 4).
Population structure was estimated by AMOVA in PopART. Individuals represented by ITS region sequences were grouped by sampling locality to populations and then nested within groups delimited by putative oceanographic barriers and host-driven isolation. We performed AMOVAs at three different levels to quantify how much variation is partitioned: (i) among individuals within populations, (ii) among populations within groups, and (iii) among groups of populations. A priori groupings were as follows: i) "Geographic groups" formed according to the geographic separation of the South African coastline [11]: "West coast group" (n = 14): localities Cape Columbine, Jacobs Bay, Granger Bay, Mouille Point, Kommetjie; "South west coast group" (n = 25): Kalk Bay, False Bay; "South east coast group" (n = 37): De Hoop, Jongensfontein, Herolds Bay, Jeffreys Bay, ii) "Barrier groups" formed in accordance with the main barriers to gene flow recognized for the populations of the type fish host Clinus cottoides (HB1, HB2 and HB3 in Fig 1) [11,12] with an exception that the most eastern locality, Jeffreys Bay, was not analyzed as a separate group due to low sample size (n = 2), so it was merged into a single group with the south coast localities To characterize the diversity of populations and their demographic history we performed several population genetic statistics for C. cottoidii groups that showed significant variation among the groups tested in the previous AMOVA (i.e. above defined geographic and barrier groups). DNA polymorphism statistics (number of haplotypes, number of segregating sites, nucleotide genetic diversity, haplotype diversity) and statistical tests aiming at detection of population size changes (Tajima's D, Fu & Li's F and D and Ramos-Onsins and Rozas R2) were calculated using DNAsp v6.10.04 [62]. The significance of all tests was determined by 10,000 coalescent simulations implemented in DNAsp. Highly significant negative values of Fs and small positive values of R2 indicate population growth whereas positive Fs values can result from balancing selection, population bottlenecks or the presence of population structure [63,64]. Significantly negative values of Tajima's D indicate recent range expansion (excess of low frequency polymorphisms), whereas significantly positive values are a signature of bottleneck or selection. A non-significant D is consistent with a population at drift-mutation equilibrium [64,65].

Ceratomyxa diversity as observed by light microscopy
Myxospores (Fig 2A-2E) and/or plasmodia (Fig 2F) of the genus Ceratomyxa were found in 40% of examined samples with the majority represented by infected gall bladders (58/152), and with few of them (3/152) found in bile ducts (liver). High morphological plasticity of Ceratomyxa spores was commonly observed in our samples with spore morphotypes similar to previously described C. cottoidii (Fig 2A) and C. dehoopi ( Fig 2B) typical for their slender spores, however, oval round spores suggestive of novel species were also frequently seen (Fig 2C and  2D). Spore shape and size (length and thickness) varied even within single fresh samples.
Aberrant Ceratomyxa spores with three polar capsules and three valves were occasionally observed ( Fig 2E). Concurrent infections of Ceratomyxa with other myxozoan genera, Sphaeromyxa and Kudoa were also detected.

Phylogenetic positioning of clinid ceratomyxids within the Ceratomyxa clade
The Ceratomyxa clade is an assemblage of diverse groups that evolved initially at different substitution rates and some of them have undergone subsequent radiation events [46]. This pattern is evident in ceratomyxids from South African clinids that clustered into a single wellsupported long-branched lineage ("Ceratomyxa group from clinid fish") which further split into several short-branched species (Fig 3). The whole lineage grouped with other extremely long-branching taxa with various scenarios in tree topology, most probably due to LBA artifacts. In ML and BI analyses (Fig 3), clinid ceratomyxids grouped with Ceratomyxa aegyptiaca from soleid fish from a Mediterranean coastal lagoon in northern Africa and Ceratomyxa longipes from a gadid fish sampled in the North Sea. The whole group then clustered either with a group including Ceratomyxa ayami and Ceratomyxa sp. M0304 (ML analysis, Fig 3) or formed a polytomy with two groups (Ceratomyxa ayami and Ceratomyxa sp. M0304; C. anko, C. appendiculata, C. pantherini, C. vikrami and Pseudalataspora kovalevae, BI analysis, BI posterior probability = 0.81). Another alternative in tree topology was observed in MP analysis, where clinid ceratomyxids grouped with C. longipes and then with the group formed by C. anko, C. appendiculata, C. pantherini, C. vikrami and Pseudalataspora kovalevae (MP bootstrap support <50%). Ceratomyxa thalassomae (EU045332)

Phylogenetic relationships within the clinid Ceratomyxa group
In this study, we amplified almost 400 clinid Ceratomyxa sequences belonging to different molecular markers of the nuclear ribosomal operon (SSU rDNA: n = 104; LSU rDNA: n = 118; Longer and shorter ITS region: n = 170, S1 Table) that were all used for subsequent phylogenetic analyses. After performing the analyses, unique sequence data representing each Ceratomyxa species from different fish hosts and localities were selected for GenBank submission (74 sequences, S2 Table). All other sequences are available in the S1 Dataset containing the untrimmed and trimmed fasta datasets of all molecular markers analyzed in this study. Blasting of newly obtained sequences against GenBank entries has shown that our data represent the first publicly available sequences of Ceratomyxa spp. from clinids. Phylogenetic analyses of SSU rDNA (S1 Fig) and ITS data (Fig 4) recognized seven well-supported novel lineages of Ceratomyxa, likely species, from South African clinids, whereas only five lineages were recovered in the LSU rDNA-based tree (S2 Fig) due to missing data from the two remaining lineages (Ceratomyxa sp. 2 and Ceratomyxa sp. 3). Based on the consensus of biological data (type host species and type locality) and, wherever possible, morphological data two of the lineages are considered to represent the previously described species C. cottoidii and C. dehoopi. The remainder of the sequences likely represents novel species (see paragraph "PCR screening, prevalence, distribution, co-infections and morphological remarks"). In the ITS-based tree, C. cottoidii clustered with a well-supported clade of Ceratomyxa sp. 1, Ceratomyxa sp. 2 and Ceratomyxa sp. 3 to create a common sister group to C. dehoopi (Fig 4). However, positioning of C. cottoidii and C. dehoopi varied depending on the marker used (S1 and S2 Figs). Another stable clade in all single gene analyses was formed by Ceratomyxa sp. 4 and Ceratomyxa sp. 5. Some members, i.e. Ceratomyxa dehoopi and Ceratomyxa sp. 3, created long branches to their sister taxa within the clinid Ceratomyxa group (Figs 3 and 4, S1 and S2 Figs).

Species distances
SSU rDNA sequence variability among clinid ceratomyxids ranged from 7.8 to 15.7% whereas intraspecific divergence ranged from 0.4% in Ceratomyxa sp. 5 to the significantly higher level (7%) in C. cottoidii (Fig 5A). A similar trend was observed for the LSU rRNA gene with interspecific variation of 18.2-28.3% and intraspecific variation ranging from 4.7 to 9% (Fig 5B). The highest variation was calculated for the ITS region showing 20.9-33.6% interspecific divergence and 1.7-18.4% intraspecific distance (Fig 5C).

PCR screening, prevalence, distribution, co-infections and morphological remarks
Species-specific PCR screening revealed 100% of samples positive at least for one Ceratomyxa species compared to only 40% parasite prevalence determined based on morphological observations (see paragraph "Ceratomyxa diversity as observed by light microscopy"). By speciesspecific PCR, C. cottoidii was recognized to be a true generalist with a high prevalence in all fish host species including its type host C. cottoides, while general clinid Ceratomyxa PCR revealed its presence in four out of five fish host species (Table 2). It is widely distributed from the west (Cape Columbine) to the east coast (Sea View) of South Africa with almost 100% prevalence at all studied localities. The only exception is its lower prevalence in Jacobs Bay (west coast) (Fig 6) where it was missing in all C. superciliosus and present in 80% of M. dorsalis samples. However, C. superciliosus was a common host of C. cottoidii at other localities ( Table 2). On the other side, C. dehoopi is a parasite with a narrower host range (three out of five clinid species) and with the highest prevalence in C. superciliosus (  Nodes with minimum support in all analyses are without labels distributional range spans from Cape Columbine (west coast) to Herolds Bay (south coast). It was not detected at any of the two east coast localities (Fig 6). The remainder of Ceratomyxa spp. is represented by generalists, present in 4-5 examined fish species and distributed from the west to the east coast of South Africa with variable prevalence at each locality ( Table 2, Fig  6). Importantly, the application of species-specific PCR was especially crucial for the assessment of the fish host species spectra of C. dehoopi, Ceratomyxa sp. 1 and Ceratomyxa sp. 5 for which general clinid Ceratomyxa PCR amplification was suggestive of their strict host specificity only in one fish species (Table 2). Co-infections of single fish host samples by several Ceratomyxa species were extremely common and occurred in 71% of all samples. One sample from C. superciliosus from Cape Columbine even contained all Ceratomyxa spp. which was confirmed by variable shapes and dimensions of myxospores observed in this sample. On the other hand, many samples with single spore morphotype observed by LM showed multiple species infections. Single infections (29% of all samples) were mainly represented by C. cottoidii (42 samples) and only two samples exclusively contained C. dehoopi or Ceratomyxa sp. 5.
We were not able to draw any conclusions regarding the measurements of C. cottoidii spores from the original species description (Reed et al. 2007) and from our samples with single infections, for which LM measurements were taken, as the comparison would be biased by potential spore shrinkage of our ethanol-fixed material (details on spore measurements in S3 Table). Unfortunately, the lack of any samples with single infections of C. dehoopi and Ceratomyxa sp. 1-5 accompanied by spore measurements prevented analysis of spore morphometry for species comparison/descriptions. Therefore, Ceratomyxa sp. 1-5 currently represent phylogenetic lineages and remain undescribed species.

Haplotype networks of Ceratomyxa cottoidii populations
The haplotype network of C. cottoidii ITS region sequences with traits defined as geographic localities was represented by 39 haplotypes and showed a considerable complexity of haplotype grouping (Fig 7A). Isolates of C. cottoidii fell into six main haplotype groups. Haplotype group A, separated from the rest of the network by 11 mutations (nucleotide changes), was represented by samples from Cape Columbine, the most western locality of our sampling range. Haplotype groups B and C were mainly created by samples from Kalk Bay (south west coast) mixed with samples from the western localities. Haplotype group D was mainly represented by samples from De Hoop and Herolds Bay (south east coast) mixed with samples from the south-western locality Kalk Bay. Haplotype group E contained exclusively samples from the south-eastern coast. Two haplotypes from the south-eastern coast locality Jeffreys Bay formed haplotype group F separated from haplotype group E by two mutations. Haplotype groups B-E showed a radial branching pattern with the main haplotype surrounded by satellite low frequency haplotypes which is typical for young expanding populations. Due to small sample size (4 and 2 sequences, respectively), existence of a similar burst pattern cannot be excluded for haplotype groups A and F (Fig 7A).
The separation of haplotypes according to the fish host species was not evident in our analyses. Haplotypes from different fish species mixed with each other and the only two unmixed groups of parasite haplotypes were biased by their exclusive origin from a single host species

Population genetic structure of Ceratomyxa cottoidii
Approximately half (55.4%) of the variation in C. cottoidii was attributable to highly significant differences within populations whereas about 42.6% was associated with significant differences among geographic groups. On the other hand, very low variation was observed among populations (2.6%). The second AMOVA showed slightly higher variation among the barrier groups (50.0%) with highly significant covariance components. The within population variation was comparable with the one from the previous analysis. The lowest level of among group variation was found for AMOVA of host groups (26.0%), the analysis reaching highly significant levels only within and among populations (  dividing C. cottoidii populations into groups delimited by major oceanographic barriers, especially if followed by the separation according to host gene flow barriers, explained significant amount of variation, however, additional factors except of fish host species, seem to impact population structuring of C. cottoidii. As revealed by AMOVA and haplotype networks, calculations of population genetics statistics for C. cottoidii using the dataset of parasite sequences from different fish species and their subsequent comparison with the population data from the type host C. cottoides is not a problem due to the lack of significant structuring of C. cottoidii populations according to fish host species.

Sequence diversity indices and demographic history of Ceratomyxa cottoidii populations
The differences between C. cottoidii populations based on geographic separation observed in haplotype networks and AMOVA were also evident in their population genetics statistics ( Table 4). Results of DNA polymorphism statistical analysis showed generally high values of haplotype diversity with small differences between the analyzed groups. From the total of 39 haplotypes, the highest number of haplotypes was found in the south west (n = 20) and south east coast (n = 20) for the geography-related groups and in B2 (n = 22) for barrier groups. The nucleotide diversity was the lowest for the south east coast (resp. B3 group). Tajima's test showed significant negative values only for south east and B3 populations whereas Fu and Li's tests revealed significant negative values for the south west, south east coast, B2 and B3 populations. The R2 test was significant only for the south east coast and B3 group populations ( Table 4). The results indicate that the south west, south east coast, B2 and B3 populations are expanding with the evidence of recent expansion in south east coast and B3 populations. In contrast, equilibrium and no population growth are suggested for the west and B1 populations at the edge of the distributional range of C. cottoidii.

Discussion
This is the first study to determine myxozoan diversity, distribution, prevalence and host specificity in a variety of clinid fish species, typical inhabitants of intertidal environments in South Africa. Moreover, the evolutionary history of the rapidly evolving clinid ceratomyxids has been assessed within this group as well as within the Ceratomyxa clade. Lastly, the most common clinid parasite in South Africa, Ceratomyxa cottoidii, has been subjected to population genetics analyses and its population structure has been compared with the data known for its type fish host, Clinus cottoides. This is a novel approach in myxozoan research as none of the Table 3. AMOVA partitioning of genetic variance for shorter ITS region within populations and groups of Ceratomyxa cottoidii. A priori groupings followed the separation of C. cottoidii sequences according to the i) geography of South African coastline (geographic groups), ii) main barriers to gene flow recognized for the populations of the fish host Clinus cottoides (barrier groups) and iii) fish host species (host groups). The covariance components are indicated as ÃÃ for highly significant (P < 0.01) and
We have revealed a high degree of species diversity represented by seven closely related Ceratomyxa spp. from the intertidal fishes Clinidae along the coast of South Africa. Such parasite radiation is probably the result of multiple speciation events triggered by limited host dispersal within this endemic area. Future application of a holistic approach in species descriptions based on detailed LM examination, proper photographic documentation and in-depth molecular identification is a good direction in accounting myxozoan species diversity, which currently seems to be underestimated [30,[73][74][75][76].
PCR screening with species-specific primers was crucial for the discovery of the remarkably high prevalence of ceratomyxids in clinid fish, with all examined samples infected with at least one parasite species and with some ceratomyxids even reaching 100% prevalence in certain hosts and at certain sites. Our findings underpin the rock pools as ideal habitats for promoting myxozoan life cycles and for the establishment of a diverse assemblage of organisms as reported previously for other animal groups [14][15][16]. Specific PCR assays were also essential for the exact assessment of Ceratomyxa spp. diversity in clinids, as frequent mixed infections, indeed common in myxozoans [17,40], complicated such evaluation. Present evidence of multiple Ceratomyxa spp. infections indicates that the extreme spore variability observed in our samples, most likely resulted from concurrent infections of closely related species of different morphology. This may also be a case for Ceratomyxa sp. from M. dorsalis [37]. Moreover, high rates of mixed Ceratomyxa spp. infections encountered in such a restricted organ space are suggestive of a low or even lacking within-host competition in this niche. Nevertheless, a temporal and spatial separation of the infections during the host's development may still occur similarly as reported for two competing marine myxozoan species [77].
In this study, we have unexpectedly revealed that Ceratomyxa spp. from South African clinids are not strictly host-specific parasites as a general rule for other ceratomyxids [25][26][27][28][29][30]. This may be caused by the close evolutionary relationships of South African clinids [8] and associated likely similar immune system traits of these fish as well as by generally low host immune response to coelozoic myxozoan parasites [78,79]. We assume that several clinid fish species are true hosts for individual Ceratomyxa spp. as spores have been observed in the site of sporogony (gall bladder) at least in several samples in which the parasite has been molecularly detected.
The phylogenetic relationships of clinid Ceratomyxa spp. with other ceratomyxids were provided in this study. However, their accurate assessment was difficult due to their unstable clustering and weak nodal supports within the Ceratomyxa clade most probably due to LBA artifact. Nevertheless, increased taxon sampling with some closely related species, possibly from other blenniiform fishes, would facilitate to break the long-branch [46,80] and thus improve inferring the exact positioning of this fast-evolving parasite group.
There is no universal criterion regarding what constitutes a sufficient level of sequence variation to represent distinct species in the genus Ceratomyxa. The intraspecific threshold values used for clinid ceratomyxids in present study are higher than commonly reported in other myxozoan clades [17] but they have been carefully evaluated with regard to the fast evolution of the studied molecular markers in this long-branched parasite group and we recommend using these cut-offs for this specific group of ceratomyxids in future studies. For example, relatively low values of maximum SSU rDNA intraspecific divergences were reported for ceratomyxids from Indian (up to 0.8%) [81], Australian (up to 1.3%) [26,27,30] and North Atlantic fishes (up to 1.6%) [46], whereas wide ranges up to 7% in C. cottoidii were encountered for South African clinid ceratomyxids (present study). Such large intraspecific variability of C. cottoidii may be a result of fast evolution of its genes concordant with increased substitution rates and associated tree branch lengths [82]. Therefore, settings of the intraspecific cut-offs and the following assessment of intraspecific and interspecific variation of a molecular marker have to be taken into account in relation to the substitution rate in a particular Ceratomyxa group. As evolutionary rates also vary for other myxozoan lineages, investigations into the sequence variation for a particular phylogenetic group are crucial for the assessment of the species concept in the Myxozoa. Moreover, it would be desirable to assess the intragenomic variation of ribosomal sequences for each myxozoan species from one individual (spore), however, low DNA yields associated with single spore extractions are limiting factors for such approach. Information about the intragenomic variation would be of interest also in the present population genetics analysis of C. cottoidii, however, distinguishing of multiple haplotype infections from the real intragenomic variability was not possible due to the fact that DNA was extracted from hundreds of spores from each infected fish individual. Moreover, spores in each host specimen may have been a genetic mixture that developed after an infection from a mixed pool of invertebrate alternate hosts harboring different C. cottoidii haplotypes. It is highly probable that intragenomic variation exists in C. cottoidii but this would only have a significant effect on "within individual level of variation" that was not assessed in our analyses. An alternative how to get around the intragenomic variability issue would be the use of mitochondrial (maternally inherited) genes, for example cytochrome c oxidase I or cytochrome b, but this approach was unfortunately not feasible in our case as such genes are difficult to amplify in myxosporeans in general and lead to several misinterpretations due to occurrence of pseudogenes (numts) [73].
C. cottoidii was shown to represent the most common parasite of clinids with a widespread distribution along the South African coastline. Investigation of a potential host-driven isolation of C. cottoidii populations, impelled by the newly discovered generalist life history of this parasite, showed no evidence for separation of parasite haplotypes by fish hosts. However, a certain bias may exist as the largest haplotype groups were characterized by a single host species. On the other hand, population structuring based on geography-related groupings (geographic and barrier groups) showed important overlapping patterns with fish host population structure [11,12] consistent with the hypothesis that reduced host dispersal, as typical for C. cottoides [10], enhances parasite community differentiation. C. cottoidii population structure is formed by young expanding populations as well as by older populations in equilibrium. The older populations are represented by haplotypes from B1 group, which includes the most western localities, most probably at the edge of the distributional range of C. cottoidii. In these peripheral populations, no haplotype mixture is evident in the associated haplotype group A, likely representing the populations under recent speciation, which further supports the clustering of Cape Columbine sample (nr. 1462) as a separate lineage within (Fig 4) or at the base of the Ceratomyxa clade (S1 and S2 Figs). The young expanding populations of C. cottoidii are mainly represented by haplotypes from groups B2 and B3 covering western, south-western and south-eastern South African localities. Though some of the western localities are present in geographically mixed haplotype groups B and C (west + south west coast) and D (south west + south east coast), the west coast haplotypes never mixed with south east coast haplotypes. This pattern of population mixing does imply existence of some gene flow barriers on one side along with still ongoing migration between south west coast and south east coast locations on the other side, which is most probably facilitated by Agulhas ring eddies between Cape of Good Hope and Cape Agulhas in the zone of contact of Agulhas and Benguela current (Fig 1). Haplotype network suggests that some additional barrier to gene flow may exist within the south east coast populations (separation of haplotype group F), however, analysis of this hypothesis by splitting the south east group into south coast (De Hoop, Jongensfontein, Herolds Bay) and east coast (Jeffreys Bay) groups was not feasible due to low sample size from the eastern locality. Therefore, a more efficient fish taxon sampling, especially from the peripheral eastern and western localities, and following haplotype sequencing are necessary to support our conclusions which would additionally provide more data for detailed gene flow and migration analyses.
Our data suggest that oceanographic barriers around Cape Point, Cape Agulhas and east of Jeffreys Bay (east of Algoa Bay) region play important roles in the distribution and population structuring of the parasite as they do for the fish host Clinus cottoides [6,11,83]. This may be caused by closely linked host-parasite co-evolution in these space-limited habitats as genetic structure and co-divergence of host-parasite populations is higher in parasites infecting hosts with limited dispersal abilities [84]. However, only 43-50% of molecular variation in C. cottoidii sequences was ascribed to differences among groups associated with geography (Table 3) while a higher percentage of variation among groups (localities) was encountered for C. cottoides populations exhibiting significant oceanographic separation [11]. Such difference in variation may be linked to the complex myxozoan life cycle involving a vertebrate (intermediate) and an invertebrate (definitive) host, which adds another variable to the dispersal and associated population structure of C. cottoidii, a myxozoan with a presumably two-host life cycle [34]. C. cottoides is most commonly associated with mid-shore areas [2,12], and while this fish host shows strong site fidelity, its ceratomyxid parasites may have a high dispersal capacity by distributing spores to other rock pools via tidal water exchange or wave splashes. An invertebrate host can be present in these rock pools or alternatively, it can inhabit other terrestrial or planktonic habitats and thus contributing further to the gene flow in C. cottoidii populations. The high infection prevalence in fish, however, suggests that the invertebrate host habitat overlaps with that of C. cottoides. Excellent candidates that are common in the intertidal rocky shore habitats are sedentary polychaetes (Annelida: Polychaeta) from the families Cirratulidae, Spionidae, Orbiniidae, Arenicolidae, Flabelligeridae, Sabellariidae, Terebellidae, Sabellidae, Serpulidae and Spirorbidae [52,85]. A systematic approach specifically targeting the species present in rock pools in order to find myxozoan parasites within these is needed.

Conclusions
We show that ceratomyxid species from South African clinids are a diverse group of fast evolving closely related parasites with high prevalence in their fish hosts, little host specificity and frequent concurrent infections, most probably as a result of radiation and no competition within the space limited host niche. Ceratomyxa cottoidii shows overlapping population structure with its type fish host, C. cottoides, however, data on the definitive host is required to unravel the complex network.
Several genetic studies of marine organisms, including C. cottoides, have shown that the sections of South African coastline coined as marine protected areas, representing hotspots of species richness and endemism and including a high diversity of habitats [86], require more protection [11,83]. As both parasite diversity and distribution are closely linked to that of the host, any change in the conservation status and distribution of the host directly impacts that of the parasite and vice versa [87]. As Ceratomyxa parasites of clinid fishes are well represented in the South African marine fauna (present study), they may significantly impact their fish host populations. Investigations into the aquatic parasite biodiversity and distribution are of highest priority, as global climate change can shift the balance in healthy parasite-rich ecosystems where parasites represent one of the most susceptible groups to environmental change [87]. Besides the implications of global climate change for fish populations, a more detailed research of myxozoan life cycles as well as the roles of these parasites in food webs and trophic transfers and their impact on the health of fish hosts are desirable not only from a South African [23] but also from a global perspective.
Supporting information S1 Table. List of clinid fish samples molecularly examined on the presence of Ceratomyxa spp. A detailed information on the host species and organ, locality, result of species-specific ITS-based PCR screening (+/-) and number of SSU rDNA, LSU rDNA and ITS PCR amplicons and/or clone sequences are available for each Ceratomyxa species along with the Gen-Bank accession numbers (corresponding to those present in the phylogenetic trees). gb = gall bladder.