Broad North Atlantic distribution of a meiobenthic annelid – against all odds

DNA barcoding and population genetic studies have revealed an unforeseen hidden diversity of cryptic species among microscopic marine benthos, otherwise exhibiting highly similar and simple morphologies. This has led to a paradigm shift, rejecting cosmopolitism of marine meiofauna until genetically proven and challenging the “Everything is Everywhere, but the environment selects” hypothesis that claims ubiquitous distribution of microscopic organisms. With phylogenetic and species delimitation analyses of worldwide genetic samples of the meiofaunal family Dinophilidae (Annelida) we here resolve three genera within the family and showcase an exceptionally broad, boreal, North Atlantic distribution of a single microscopic marine species with no obvious means of dispersal besides vicariance. With its endobenthic lifestyle, small size, limited migratory powers and lack of pelagic larvae, the broad distribution of Dinophilus vorticoides seems to constitute a “meiofaunal paradox”. This species feasts in the biofilm among sand grains, but also on macroalgae and ice within which it can likely survive long-distance rafting dispersal due to its varying lifecycle stages; eggs encapsulated in cocoons and dormant encystment stages. Though often neglected and possibly underestimated among marine microscopic species, dormancy may be a highly significant factor for explaining wide distribution patterns and a key to solving this meiofaunal paradox.

markers and employed methods. GMYC analyses indicated the same two minimum entity estimates in all gene tree combinations with the exception of four instances where estimates suggested more than two identifiable entities (see Table 1). All GMYC analyses were statistically significant (p ≤ 0.05). Analyses using mPTP (BEAST trees) estimated two phylogenetic entities throughout all datasets tested (Table 1). When using bPTP (BEAST trees), two of the multiple gene tree combinations (COI + CytB + 16S and COI + CytB + 16S + 18S) recovered more than two entities. Comparative PTP analyses in RAxML likewise resulted in two phylogenetic entities for all larger compiled datasets as well as COI. Yet, our RAxML analyses failed to recover the minimum two entities in five (out of 23) of the smaller data sets (see Table 1). The variation in minimum entities recovered is most likely due to differences in tree building algorithms and how missing data is treated. ABGD analyses of the single gene datasets (COI and CytB) again recovered two distinct phylogenetic entities (Table 1).
Although we base our conclusions on the species delimitation analyses listed above, we also calculated genetic distances with Mega v7.0 42 using the Kimura 2-parameter model with variation among sites modeled with a gamma distribution (shape parameter = 1). Positions containing missing data were eliminated. The number of base substitutions per site between sequences were highest for the sequenced COI fragment (645 bp) with comparable similarities among specimens of the 11 localities of D. vorticoides (97.2-100%) and among the single sampled population of D. taeniatus (97-99.7%), but with distances between these two sister clades being ten times higher (82.3-85.6% similarity) (Suppl. Table 2). Similarities were comparable for CytB (405 bp), ranging from 98.4-100% among D. vorticoides specimens, 96.6-98.5% among D. taeniatus specimens and 81.9-85.4% between these two sister species (Suppl. phylogenetic analyses. Our phylogenetic investigation of Dinophilidae was conducted employing both Bayesian probabilities and Maximum likelihood analyses for all single gene datasets as well as our final concatenated five gene dataset (18S rDNA, 28S rDNA, 16S rDNA, COI, CytB). All analyses recovered Dinophilidae monophyletic, containing three well-supported clades found to represent differing morphologies and life cycles, but with Dinophilus being recovered paraphyletic (Fig. 2). Marked clades C, D, and A in Fig. 2 Fig. 2). The internal resolution within the D. vorticoides clade was not congruent among the analyses due to the low level of diversity among sequences.
The hyaline taxa of Dimorphilus (Clade A, Fig. 2) are characterized by having dwarf males, fast life cycles enabling year-round reproduction, and females with transverse rows of dorsal body ciliation. None of these taxa were wild caught, but sampled either from marine station aquaria or stemming from laboratory cultures (Suppl. Table 1). Dimorphilus cf. gyrociliatus cultures from Japan and D. gyrociliatus cultures from Italy formed a clade (BPP = 1.0; MLB = 100), which was found to nest with D. cf. kincaidi (BPP = 0.99; MLB = 72); the three of them forming a sister clade to D. gyrociliatus, Xiamen, China (Fig. 2).

Morphological examinations and taxonomic implications. Specimens, both live and fixed, belonging
to the newly defined Dinophilus clade were examined using LM and SEM (Fig. 3) and all found to be strongly orange pigmented and sexually monomorphic. They have elongated cigar-shaped bodies with six indistinct trunk segments, a relatively broad mid-ventral ciliary band, and a prostomium with two pairs of anterior compound cilia and at least two incomplete transverse ciliary bands. Dinophilus vorticoides from West and South Greenland, Faroe Islands, Svalbard (Norway), Sweden, White Sea (Russia); D. taeniatus from Cornwall (United Kingdom); and Dinophilus sp. from Yucatán (México) show no obvious morphological differences and all possess two transverse ciliary bands per trunk on segments 1-5, with a dorsally incomplete transverse band on segment six. In contrast, the anterior prostomium and trunk of D. gardineri is densely ciliated with individual ciliary bands only distinguishable in the posterior trunk. These observations together with the phylogeny results in the following generic definitions.
Type species. Dimorphilus gyrociliatus (O. Schmidt, 1857). Type locality: Naples, Italy (sampled from algae growing on the piers and rocky shore off Santa Lucia).
Diagnosis. Females with hyaline body with dark red pigmented kidney-shaped eyespots. Prostomium with two dorsally incomplete transverse ciliary bands. Six trunk segments, each with single transverse ciliary band. Strongly dimorphic dwarf males with minute round bodies, no distinct body regions or segmentation, presence of anterior, ventral and posterior ciliation and a muscular copulatory organ. Female life cycle completed within three weeks, males within a week.
Etymology. From Greek "dimorphos" (from di-'twice' + morphē 'form'), to account for the dimorphic sexes of the genus containing dwarf males, and from Greek "philos" ('liking of ') in accordance with the similar ending of the type genus Dinophilus.
Remarks. Besides from D. gyrociliatus, only D. kincaidi is currently regarded a valid species within Dimorphilus based on morphology. Several species have been considered invalid due to poor descriptions and doubtful resemblance to dinophilid annelids 43 , often representing platyhelminths instead, e.g., D. sphaerocephalus Schmarda, 1861, D. borealis Diesing, 1862 and D. rostratus Schultz, 1902. However, species of Dimorphilus (especially D. gyrociliatus) have been reported from multiple disjunct localities around the world, some of which may represent genetically distinct species. For example, Dimorphilus apatris (Korschelt, 1882) and D. conklini (Nelson, 1907) were regarded junior synonyms to D. gyrociliatus due to morphological similarity 31  Emended diagnosis. Strongly orange pigmented, monomorphic, with kidney-shaped red pigmented eyes. Prostomium, buccal region and trunk with transverse ciliary bands; antero-dorsal trunk ciliation dense or restricted to double segmental bands; posterior segments with one-two transverse ciliary bands each. Life cycle long with prolonged encystment stage.

Remarks.
Aside from D. vorticoides, D. gardineri (New England, USA) shows morphological differences and must be regarded a valid species. This study also proves D. taeniatus (Plymouth, UK) to represent a genetically different entity (opposing the previous synonymy with D. vorticoides 40 ). Dinophilus jägersteni Jones and Ferguson, 1957 (US East Coast) shows some discrete differences in ciliation to D. gardineri and future genetic studies may  31 , but these decisions may also have to be reevaluated by incorporating molecular data. Dinophilus simplex Verril, 1892 is invalid and was rejected by Ruebush 44 .

Discussion
Molecular investigations of marine meiofauna generally argue against the ubiquity theorem (EiE) [20][21][22][23][24][25] , and the few meiofaunal annelids proposed to have broad amphi-Atlantic distributions have all been disproven cosmopolitans when investigated genetically 7,23 . Our recovered continuous boreal-Atlantic distribution of D. vorticoides is exceptional towards supporting the "Meiofauna paradox", allowing us to now revisit controversial hypotheses, both on evolutionary stasis 10 and EiE for marine meiofauna.
Broader distributions of meiofauna are now mainly attributed to contemporary events (e.g., rafting, current, drift), integrating processes of dispersal, over that of supercontinent populations. However, while some limno-terrestrial meiofauna seems to fit the EIE hypothesis 8 due to desiccant-tolerant stages capable of long-distance dispersal, desiccation tolerance in marine meiobenthos is very rare. Nonetheless, size and dormancy in various forms may still be a relevant trait for explaining wide distribution of marine meiobenthos, integrating processes of long-distance dispersal by means of megafauna, rafting, or even anthropogenic means of ballast water and sand 23,[45][46][47] . Dinophilids constitute a highly derived evolutionary lineage within Annelida [48][49][50][51] and show distinct and well supported clades and subclades (Fig. 2) contradicting evolutionary stasis. Moreover, D. vorticoides shows a mix of populations, especially among the Faroese and Greenland localities (see Fig. 1), whereby we attribute its broad distribution to a remarkable and previously inconceivable dispersal ability.
Dinophilidae, like other direct developing microscopic fauna, are considered to have low dispersal potential 10 and the low genetic diversity among D. vorticoides populations recovered herein is counterintuitive given their direct development (=lacking pelagic larval stages) and their limited migratory abilities, being small and moving by ciliary action. On the other hand, D. vorticoides and D. taeniatus have an encystment stage extending from July to October/November in Sweden and UK; more northern populations have delayed onset and excystment extended sometimes into late winter. Both cysts and eggs encapsulated in gelatinous cocoons 35,36 attach to substrates and are likely more tolerant to fluctuations in temperature and salinity than free moving adults 35 . Knowingly, all of Dinophilidae are described from interstitial sediments 31 , yet D. vorticoides, and likely other dinophilids, are not strictly bound by them 33,36 . Throughout our sampling, especially in the Faroe Islands, D. vorticoides was regularly found inhabiting filamentous algae and buoyant mollusk eggs, the later association not previously described. Seemingly, their willingness to selectively graze biofilms on drifting algae, buoyant invertebrate clutches, and potentially other marine debris (natural or anthropogenic) opens up their means of dispersal by ocean currents and winds 14 , especially in the egg or encystment phase of their life cycle. In the North Atlantic, floating algal clumps have long been associated as a type of microcosm or 'micro-island' that have the ability to transport and disperse macroscopic invertebrates, including larger annelids [45][46][47] . Although studies are limited and usually overlook micrometazoans, algal rafts are capable of traveling hundreds of kilometers, but the diversity of their faunal hitchhikers decreases with distance 45 . From a historical perspective, floating algae in European waters were responsible for seeding current algal communities throughout the Faroe Islands 46 and is a widely accepted means by which the rocky shore flora and fauna of Eastern North America was established post glaciation 45 . Collections of D. vorticoides on filamentous algae are quite telling towards the likelihood of such a transport mechanism, and unlike most intertidal invertebrates, their internal copulation with later cocoon deposition and an encysted summer/fall (in the far north also including early winter) stage would potentially increase their tolerance during transportation, aiding towards successful colonization events and erection of small populations within short intervals.
There are however obstacles still obstructing our understanding of the processes shaping the current and seemingly cosmopolitan distribution in D. vorticoides, or with the same notion, what is restricting D. taeniatus from dispersing into the colder, northern localities of D. vorticoides? Our phylogenetic delineations suggest an absence of reproductive barriers among D. vorticoides populations, corroborated by low genetic diversities; however, identifying evidence of dispersal is far more convoluted. In Guil's 9 review of micrometazoans, it was discussed that dispersal may not always translate into gene flow, suggesting that numerous ecological and organismal conditions need to be met prior. One of these a priori conditions is that a continuous wave of individuals would be reaching each of the locations, however, our collection sites of D. vorticoides are often inundated with sea ice, altering normal colonization pathways. Furthermore, the polar ice coverage in the latest glacial period engulfed several of the current locations of D. vorticoides, now spanning relatively young oceanic areas. Interestingly, Jägersten 35 observed D. vorticoides cysts encased in ice that continued to hatch upon thawing. Based on this indirect evidence of colonization along the retracting ice edge of the latest polar glacial coverage and the fact that our genetic distances within D. vorticoides has larger ranges than other meiofaunal investigations 7,23 , it would appear our distributions are heavily influenced by seasonal climatic events. Similarly, these events, including currents, may be restricting distributions of D. taeniatus, as northern Atlantic waters appear to not regularly mix with those surrounding the United Kingdom and/or D. taeniatus may not be able to tolerate the lower water temperatures in more northern waters. An interesting notion for future research is that genetically undetermined D. vorticoides/D. taeniatus populations have also been reported from more southern localities in e.g., Roscoff and Valencia, along the French and Spanish Atlantic coast, respectively 39 . It remains to be tested whether these populations represent the D. taeniatus species, which may hereby likewise have a broad but more southern distribution. While we are limited in our broad understanding of meiofaunal distribution in general, it appears that within Dinophilidae, especially D. vorticoides, represents a true cosmopolitan species, supporting the EiE (but the environment selects) hypothesis. While numerous questions still remain, our findings suggest that both historical and contemporary events are shaping distributions patters in meiofaunal annelids and likely to be group or even species specific. Our findings inevitably will provide a basis for future investigations whereby a more integrative approach focusing on connectivity may finally help in elucidating the meiofaunal paradox, sensu stricto. www.nature.com/scientificreports www.nature.com/scientificreports/ Using most conservative estimates, species delineations of single and combined gene datasets have identified two significant Dinophilus clades throughout the North Atlantic that are more closely related to each other than to other Dinophilus taxa from the eastern United States and México. We hereby reject the previously proposed synonymy of D. vorticoides and D. taeniatus 31,40 , and accept the validity of both species, now recognizing D. vorticoides as the prolific species of the North Atlantic.
Dinophilidae is now demarked by three well-supported clades with distinct reproductive modes that eventually lead to the formation of gelatinous cocoons (in Dinophilus and Dimorphilus) or clutches (in Trilobodrilus) that house yolky eggs. These cocoons or clutches are deposited in favorable habitats, including sediments, on algae, and along rocks or pilings, giving rise to directly developed free-swimming, but benthic, juveniles (Fig. 3). Systematically, these clades include Dimorphilus gen. nov., representing the smallest dinophilids (females ≤ 1.3 mm), displaying strong sexual dimorphism with completely hyaline dwarf males and females with a single ciliary band per segment and a rapid life cycle; Trilobodrilus (≤2.0 mm), being hyaline with a distinctive trilobed prostomium and limited ciliation on the trunk 30,31,52 ; and Dinophilus, being the largest dinophilids (≤3 mm), with easily recognizable orange-red pigmentation and two transverse ciliary bands on each trunk segment or a more random and occasionally denser distribution of ciliary tufts 34,35,40,44 . From a strictly morphological standpoint, Dimorphilus (e.g., D. gyrociliatus) shares several diagnostic morphological features of Dinophilus (e.g., D. taeniatus), so their previously proposed affinity had never come into question.
By incorporating years of meiofaunal collections, our phylogenetic analyses of Dinophilidae now provides clear support of separate evolutionary pathways between monomorphic and dimorphic taxa, promoting a common monomorphic ancestor for the clade Trilobodrilus -Dinophilus. Unfortunately, we cannot provide further insight whether the traits of the paedomorphic Dimorphilus are more closely representing the ancestor of all Dinophilidae or present a series of derived characters. Prior to this study, the most extensive molecular phylogenies within Dinophilidae were focusing on the systematics within Trilobodrilus 1,53 , and regrettably, these phylogenies only included a single Dinophilus species (syn. as Dimorphilus). While Kajihara et al. 53 did include Dinophilus sp. from Lizard Island (Australia) (previously published by Worsaae & Rouse 54 ), it had initially been incorrectly identified, thus preventing the discovery of a paraphyletic Dinophilus. This specimen was herein reexamined by means of LM and reidentified as Trilobodrilus sp.
In addition to published results, the composition of Dimorphilus, while fully supported (Fig. 2), is confusing since cultures once maintained by Bertil Åkesson of D. gyrociliatus collected in China, Italy, etc. have since been dispersed, making it unfeasible to trace if mixing of cultures has occurred, or which genetic sample stem from the type locality. While this was not the focus of the study, it appears that multiple independent lineages of D. gyrociliatus are being used throughout the literature (see David & Halanych 55 ), and our integration of sequences from both marine aquaria and laboratory cultures suggest that a larger and unidentified diversity is currently present within Dimorphilus.

Methods
In order to determine the extent of the distribution of D. taeniatus, samples were collected from coastal areas throughout the Faroe Islands, Greenland, México, Sweden, Norway, and the White Sea, Russia. Specimens were also obtained from the type locality of D. vorticoides in the Faroe Islands and for D. taeniatus from the United Kingdom in order to determine the validity and phylogenetic extent of these species. collection and examination. Interstitial members of Dinophilidae were extracted from fine sand and coral rubble 56 collected from the intertidal zone to 20 meters depth. Epibiont and epibenthic specimens were collected by hand from rocks, algae or pilings.
Light microscopy (LM) was used to examine newly acquired live material using an Olympus IX70 inverted microscope mounted with an Olympus DP73 digital camera. Detailed morphological examinations were made using a JEOL JSM-63335F field emission scanning electron microscope (SEM) at the Natural History Museum of Denmark, University of Copenhagen. Specimens were prepared for SEM following previously published protocols that included fixation in glutaraldehyde, postfixation in osmium tetroxide, and dehydration by an ascending ethanol series 30,54 . Prior to imaging, all specimens were critical-point dried, mounted on aluminum stubs, and sputter-coated with platinum/palladium.

Molecular laboratory methods.
Evolutionary relationships within Dinophilidae were examined using the ribosomal markers 18S rRNA, 28S rRNA and 16S rRNA, as well as the mitochondrial markers cytochrome c oxidase subunit I (COI) and cytochrome B (CytB). The combinations of conserved and fast evolving genes were selected to resolve inter-and intraspecific relationships among and between the dinophilid genera. Additionally, given that the three selected mitochondrial markers are fast evolving, they were also selected to resolve population level dynamics within the 'D. taeniatus/D. vorticoides clade' .
Total genomic DNA was obtained from individual dinophilid specimens using the Qiagen DNeasy Tissue & Blood Kit (Qiagen Inc., Valencia, CA, USA) following the manufactures protocol. DNA was extracted from at least three separate individuals from each of the collection localities of 'D. taeniatus/D. vorticoides' .
Samples were prepared for polymerase chain reactions (PCR), sequenced and aligned according to methods previously outlined 30 .
Generated sequences were deposited in GenBank ® and their accession numbers can be found in Suppl. www.nature.com/scientificreports www.nature.com/scientificreports/ Dataset assembly. To understand the relationships among the genera of Dinophilidae and identify the closest relative of the D. vorticoides/D. taeniatus clade, phylogenetic reconstructions for Dinophilidae were performed using gene data from only a single representative from each collection locality as well as any already deposited information available on GenBank. For species delineations within the 'D.vorticoides/D. taeniatus clade' , identical sequences were first identified using pairwise distances in Bioedit 57 and subsequently removed to avoid inclusion of redundant information.
Due to the fact that the position of Dinophilidae is still highly debated and unresolved, outgroup selection was based on recovered sister group relationships from recent phylogenomic investigations 48,49  Sequences were aligned using the MAFFT online platform 58 . Individual gene datasets were aligned under the L-INS-I interactive refinement method 59 . Datasets for 18S rRNA and 28S rRNA were aligned with the 'nwildcard' option selected, as this does not designate missing data as gaps. Alignments of COI and CytB were trivial, however, both datasets were aligned to check for directionality. Protein coding genes COI and CytB were checked for stop codons prior to phylogenetic analyses using Mesquite v.3.51 60 . Individual gene datasets were concatenated using Sequence Matrix 61 . phylogenetic analyses. Phylogenetic reconstructions were performed on individual gene datasets, as well as concatenated gene datasets, using both maximum likelihood (ML) and Bayesian methods.
Maximum likelihood analyses were performed in RAxML v.7.2.8 62 as implemented on the CIPRES Science Gateway 63 . Given that RAxML only implements general time reversible (GTR) models of sequence evolution for amino acids, a GTR model with corrections for discrete gamma distribution (GTR + Γ) was specified for individual gene and concatenated gene datasets. Non-parametric bootstrapping with 1,000 replicates was used to generate nodal support estimations 64 .
Bayesian analyses (BA) were performed using MrBayes v.3.2.6 65 as implemented on the CIPRES Science Gateway 63 . Prior to analyses, jModelTest 66 was used on each individual gene dataset (18S rRNA, 28S rRNA, 16 s rRNA, COI, CytB) to evaluate their optimal evolutionary model as estimated by the corrected Akaike information criterion (AICc). A GTR model with gamma distribution and a proportion of invariable sites (GTR + I + Γ) was shown to be the best estimate for 18S rRNA, COI, and CytB, while 28S rRNA and 16S rRNA were selected for a GTR + Γ model. Both individual and concatenated datasets were run with two independent analyses using four chains (three heated, one cold). Generation sampling was set to 30 million, sampling every 1000 generations. Burnin was set to 10 million generations. Majority-rule consensus trees (50%), posterior probabilities, and branch lengths were constructed using the remaining trees after burnin. Convergence of all MCMC runs were verified using TRACER v.1.6.0 67 .

Species delimitation.
We employed an integrative taxonomic approach, including morphological (LM and SEM) and DNA taxonomy to determine if D. taeniatus represented a single species throughout the sampled localities, or at the opposite extreme, if each population could be considered separate and an independently evolving entity. Three methods widely employed in DNA taxonomy were used 18 , including the generalized mixed Yule-coalescent (GMYC) model 68 , the Poisson tree process (PTP), including multi-rate (mPTP) 69 and Bayesian implementation (bPTP) 70 , and the Automated Barcode Gap Discovery (ABGD) 71 . Outgroups were removed prior to implementation of the beforementioned methods.
Methods of GMYC, mPTP, and bPTP utilized ultrametric trees generated using Bayesian Inference in BEAST (see below). Since BEAST generates ultrametric trees without smoothing (no data loss), all analyses incorporated BEAST trees to maintain consistency across methods used when generating tree topologies. Yet, to test the strength of our mPTP and bPTP delineations, these analyses were additionally run with methods of maximum likelihood, generating trees in RAxML v.7.2.8 62 (following methods listed above). Analyses using GMYC were performed in R v.3.5.1 (R Core Team, 2014) using the package SPLITS v.1.0-19 72 on phylogenies obtained from individual and combined datasets. PTP analyses were carried out on the mPTP online server (http://mptp.h-its. org) and the bPTP online server (http://species.h-its.org). No changeable settings are present on the mPTP online server, however, bPTP analyses were run using 10 4 MCMC generations with a burning of 0.1. Individual COI and CytB datasets were uploaded on the ABGD online platform (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb. html) and were analyzed using preset parameters 70 .
Ultrametric trees were generated using Bayesian Inference in BEAST v.1.8.4 73 as implemented on the CIPRES Science Gateway 63 . BEAUTi (Bayesian Evolutionary Analysis Utility) v.1.8.4 generated xml files for all BEAST runs. Independent BEAST runs were created for individual and well as combined and partitioned datasets. Tree priors for all analyses were selected under a Coalescent Process with constant population size. Nucleotide substitution models were estimated by the corrected Akaike information criterion (AICc) using jModelTest 66 . Both 18S rRNA and COI were selected for a generalized time reversible model with a proportion of invariable sites (GTR + I), 28S rRNA and 16S rRNA were selected for GTR with gamma distribution (GTR + Γ), and CytB was selected for Hasegawa, Kishino, and Yano model with gamma distribution (HKY + Γ). All datasets had independent MCMC analyses with 10 8 generations and trees were sampled every 10000 generations. TRACER v.1.6.0 67 was used to verify convergence of all MCMC runs. A maximum clade credibility (MCC) consensus tree was obtained for each BEAST dataset in TreeAnnotator v.1.8.4 after annotating the remaining 9001 trees after burnin.