Deceleration of morphological evolution in a cryptic species complex and its link to paleontological stasis

Morphological stasis or the absence of morphological change is a well‐known phenomenon in the paleontological record, yet it is poorly integrated with neontological evidence. Recent evidence suggests that cryptic species complexes may remain morphologically identical due to morphological stasis. Here, we describe a case of long‐term stasis in the Stygocapitella cryptic species complex (Parergodrilidae, Orbiniida, Annelida). Using phylogenetic methods and morphological data, we find that rates of morphological evolution in Stygocapitella are significantly slower than in closely related taxa (Nerillidae, Orbiniidae). Assessment of quantitative and qualitative morphology revealed the presence of four morphotypes with only subtle differences, whereas molecular data supports 10 reproductively isolated clades. Notably, estimates for the time of Stygocapitella species divergence range from ∼275 million years to ∼18 million years, including one case of two morphologically similar species that have diverged about 140 million years ago. These findings provide evidence for morphological deceleration and long‐term morphological stasis in Stygocapitella, and that speciation is not necessarily accompanied by morphological changes. The deceleration of morphological divergence in Stygocapitella can be potentially linked to niche conservatism and tracking, coupled with the fluctuating dynamics of the interstitial environment, or genetic constraints due to progenetic evolution. Finally, we conclude that failing to integrate speciation without morphological evolution in paleontology may bias estimates of rates of speciation and morphological evolution.

The occurrence of morphological stasis, defined as little or no morphological evolution over extended periods of time, remains a controversial topic in evolutionary biology (Futuyma 2010). Morphological variation is seen as a desired feature of any biological system and its absence is often interpreted as a potential failure to capture variation (Weiss 2011). Lineages with higher evolvability are expected to occupy broader range of habitats more quickly and efficiently, ultimately replacing less labile groups (Rabosky and Adams 2012) and hence, cases of low evolvability and morpho-logical stasis are expected to be exceptionally uncommon. Despite this, examples of morphological stasis are commonplace in the fossil record (Cheetham 1986;Futuyma 2005Futuyma , 2010Frame et al. 2007;Hunt et al. 2008;Hunt and Rabosky 2014;Voje et al. 2018), where series of invariant morphotypes occur at diverse time-scales in different organismal groups (Cheetham 1986;Hunt 2007).
A theory that aims to explain the occurrence of long periods of stasis is the punctuated equilibrium (Eldredge 1971;Eldredge and Gould 1972). In its essence, punctuated equilibrium suggests that species undergo long periods of morphological stasis, which are disrupted by rapid change during speciation (Eldredge and Gould 1972;Futuyma 2005). The postulation that morphological evolution occurs exclusively during speciation implies that adaptive and selective processes are insignificant during substantial parts of the evolutionary histories of species (e.g., Stanley 1975), thus challenging the accumulating evidence of the emerging "modern synthesis" (Futuyma 2005;Hunt and Rabosky 2014). Although the modern synthesis-punctuated equilibrium debate lasted for about two decades, over the years paleontological evidence was aligned with the major processes suggested by the modern synthesis: selection, drift, mutation, and gene flow (Hunt and Rabosky 2014). However, the conciliation of stasis with these processes was never thoroughly achieved. On one hand, researchers have argued that stasis could result either from biases in the paleontological evidence due to, for example, taphonomical biases (Kidwell and Holland 2002;Pennell et al. 2014), or either be due to the lack of statistical or sampling resolution (Frame et al. 2007;Voje 2016). On the other hand, competing views have focused on developing theoretical frameworks underlying the deceleration of morphological evolution, which include scenarios of stabilizing selection (Charlesworth et al. 1982;Hansen and Houle 2004;Futuyma 2010), niche conservatism and tracking (Futuyma 2010(Futuyma , 2015, fluctuating ecological conditions (Futuyma 1987(Futuyma , 2010(Futuyma , 2015Sheldon 1996;Smith et al. 2011), lack of new ecological interactions (Nordbotten and Stenseth 2016), constraints (Charlesworth et al. 1982;Maynard Smith et al. 1985;Wagner and Schwenk 2000;Hansen and Houle 2004;Futuyma 2010;Smith et al. 2011), recurrent bottlenecks (Futuyma 2010), physiological or behavioral adaptation (Lee and Frost 2002;Futuyma 2010;Lassance et al. 2019), and the influence of particular environments and environmental conditions (Westheide 1977;Futuyma 1987Futuyma , 2010Westheide and Rieger 1987;Giere 2009;Gueriau et al. 2016).
One suggested way of integrating components of punctuated equilibrium and modern synthesis can result from variation in rates of anagenetic and cladogenetic change (Futuyma 1987(Futuyma , 2015. If a cladogenetic event (speciation) occurs, and posterior anagenetic changes are slowed or nonexistent, daughter species will be similar in morphology. Although this idea is plausible, it cannot be tested in the paleontological record because paleontological species are diagnosed based on morphology alone and rely on the evolution of morphological differences (Jackson and Cheetham 1999;Aze et al. 2011;Strotz and Allen 2013). If daughter species are morphologically similar, these estimations will be biased. A solution to this problem is the combination of molecular phylogenetic tools coupled with morphological data in lineages displaying stasis (Bokma 2002(Bokma , 2008Mattila and Bokma 2008;Katz 2018). For example, Katz (2019) argues that integrating evidence from paleontology, paleobiology, and molec-ular phylogenetics is the key to understand the early-burst of the Angiosperms. Applying a holistic approach to the problem of stasis could potentially provide a link between paleontological evidence and evolutionary studies based on neontological evidence. Cryptic species complexes are a potential target for such an approach. Morphological stasis has been suggested to occur in extant cryptic species complexes-different species that are similar in morphology (Wada et al. 2013;Swift et al. 2016;Cerca et al. 2018;Struck et al. 2018). Cryptic species have been found widespread across the tree of life (Pfenninger and Schwenk 2007;Pérez-Ponce de León and Poulin 2016), and are being discovered at a faster pace. Matching morphological stasis as seen in cryptic species and paleontological stasis may provide the possibility to test for phylogenetic signatures of stasis in closely related extant taxa (Mattila and Bokma 2008) and offer the possibility to quantify gene flow and genetic divergence in these complexes (Sheldon 1996;Futuyma 2010;Hunt and Rabosky 2014). In cryptic species, morphological evolution can be potentially decelerated, if adaptive pressures focus on physiological, behavioral, or biochemical traits that have no bearing on morphology (Gómez et al. 2002;Lee and Frost 2002;Novo et al. 2010Novo et al. , 2012Struck et al. 2018;Lassance et al. 2019). Despite this potential, currently described cryptic species complexes are relatively young, contrasting with fossil record evidence that suggest long-lasting stasis.
A promising system to close this gap is the Stygocapitella cryptic species complex (Parergodrilidae, Orbiniida, Annelida), which until recently was thought to consist of a single species with a worldwide distribution (Schmidt and Westheide 2000;Struck et al. 2017). Recent evidence suggests that pronounced periods of stasis are intertwined with slight changes in morphology (Struck et al. 2017). However, this work focused on a limited set of species and populations with slight morphological differences, hence lacking resolution to describe morphological and genetic evolution. Here, we present an extended sampling of this genus, with special emphasis on the Northern hemisphere. Specifically, we determine genetic and morphological divergences in the Stygocapitella cryptic species complex and in closely related taxa.

PRESERVATION
Stygocapitella is an interstitial annelid, generally found around the high-water line of stable, sheltered gravel, or sandy beaches (Purschke 2018;Fig. 1). In each sampling location, we collected sediment samples by drawing transects from the high-water line to the foot of the dune. Every 1 m, we dug a 1-m deep hole and, in each of these, we collected approximately 375 cm 3 of sediment every 15-cm depth. Afterward, we extracted Stygocapitella specimens using the MgCl 2 method and a dissecting microscope green (Eastern Atlantic species C; EA C). Circles with several colors identify sympatry. (Westheide and Purschke 1988). For this study, a total of 33 sites across the Northern hemisphere were sampled ( Fig. 2; Tables S1 and S2). Specimens used for molecular biology were preserved in 95% ethanol and for morphological analyses in sucrose-picric acid-paraformaldehyde-glutaraldehyde (SPAFG) following Westheide and Purschke (1988) (Table S1, S3, and S4).

MOLECULAR BIOLOGY
Genomic DNA was extracted using phenol-chloroform or the EZNA Tissue DNA Kit (Omega Bio-Tek  Struck et al. 2002), and 18SF3 Stygo (CCTCGGGATTGGAATGAGTAC; Struck et al. 2017). After sequencing, all sequences were assembled and the ends automatically trimmed using Geneious (version 6.8.1). The quality of the assembly of each sequence was visually checked and each consensus sequence screened for contamination by doing NCBI megablast searches against the complete nonredundant database.

PHYLOGENETIC AND MOLECULAR CLOCK ANALYSES
Details about amplified markers as well as out-groups obtained from GenBank are provided in Table S1. The sequences of each gene were aligned with mafft version 7.310 using a maximum of 1000 iterations and the accurate localpair algorithm (Katoh and Standley 2013). ITS1 was aligned using the genafpair algorithm, which is optimized for gappy sequences. After alignment, both ends of the sequences were trimmed until the first position without missing data. These datasets were concatenated with FAS-conCAT version 1.1 (Kück and Meusemann 2010). A partitioned maximum likelihood (ML) analysis of the concatenated dataset was performed using IQ-tree version 1.5.5 (Nguyen et al. 2015) with an automatic determination of the best substitution model for each gene, 300 initial parsimony trees, 15 best trees retained during search, and 1000 bootstrap replicates. Similar settings were used for ML analyses of each gene independently. Molecular clock analyses were conducted using BEAST version 2.4.7 (Bouckaert et al. 2014). We used IQ-tree's Mod-elFinder to determine which substitution models best fit the dataset (Kalyaanamoorthy et al. 2017). After performing several runs using combinations of different prior models and genes, we selected CO1 and 18S for the final dating analysis because both genes are commonly used for molecular dating of annelids and showed best chain convergence and effective sampling sizes. BEAST was run with linked trees using a TN93 model with equal frequencies for 18S and HKY with estimated frequencies and a -distribution with four categories for CO1. As no fossil record is known for these taxa, we selected a relaxed, log-normal clock to account for rate heterogeneity across lineages and mean values of 0.0001425 and 0.0176 as substitution rates per million year for 18S and COI, respectively (Escalante and Ayala 1995;Pérez-Losada et al. 2004;Struck et al. 2017). We applied a birth-death model and constrained the in-group as monophyletic. A Markov Chain Monte Carlo was run for one billion generations, sampling every 100,000 th generation. We confirmed chain convergence using Tracer version 1.6 (Rambaut et al. 2007), with a 1% burn-in threshold. A Maximum Credibility Consensus Tree was obtained using TreeAnnotator (Bouckaert et al. 2014). Considering the potential biases and criticism of molecular clock approaches, we limit our interpretation of these results.
Polymorphic sites, haplotype diversity, sampling variation, and Tajima's D were calculated using DNAsp version 6.10.01 (Rozas et al. 2017) and the 16S marker, as it has the highest coverage of species and populations. For the isolation-by-distance test, we applied a mantel test in R using the "ade4" package (Dray and Dufour 2007). For this test, we selected the sequences for each of the four Atlantic species separately, and realigned each individual gene dataset using mafft as previously described. Pairwise F ST were obtained with DNAsp. The least distance between each pair of sites was calculated using google maps' "measure distance" function. After this, we modified the distance line by adding points until the minimum distance between two sites would not cross any landmass. The obtained F ST values and geographic distances between sites were used to compute a mantel test based on 9999 replicates using the "ade4" R package (Dray and Dufour 2007). For the EA C, we excluded the two American specimens to ensure that the test is not dominated by the long distance to these two specimens, as evidence shows that these might be a result of accidental human translocation (Radziejewska et al. 2006).

MORPHOLOGICAL FIXATION AND PRESERVATION
Morphological divergence was evaluated by scanning electron microscopy (SEM) and morphometrics. Due to the small size of the specimens, we were unable to genotype and obtain morphological data of the same individual (Westheide and Hass-Cordes 2001). Individuals used in morphological studies were identified indirectly based on the investigated site. For SEM, specimens were transferred to a buffered 1% OsO 4 solution for 1 h at ambient temperature and dehydrated in a graded ethanol series from 30% to 100%. Dehydrated specimens were critically point-dried with CO 2 , mounted on aluminum stubs, sputter coated with platinum, and examined with a Zeiss Auriga field emission SEM. In total, we investigated 73 specimens from 16 sites (ß4.5 specimens per from: Langebaan, Sarge Bay, Gnarabup Beach, Roche Harbor, Reuben State Park, 4 th of July Beach, Canoe Beach, Lubec, Lødingen, Henningsvaer, Schilksee,Île Callot, List, Bristol Channel, Gravesend, and Plymouth; Table S3). For morphometrics, quantitative measurements of sexually mature specimens were carried out based on pictures taken under a light microscope at an amplification of 10×. Images were assembled using Adobe Photoshop and morphometric traits measured in ImageJ. These included six body size characteristics: body length and width, prostomium length and width, and pygidium length and width; of a total of 123 individuals from 18 sites (ß7 specimens per site; Langebaan, Sarge Bay, Gnarabup Beach, Roche Harbor, Reuben State Park, 4 th of July Beach, Canoe Beach, Lubec, Lødingen, Henningsvaer, Schilksee,Île Callot, Glenancross, Keitum, Bristol Channel, Gravesend, Plymouth, and Ellenbogen; Table S4).

RATES OF MORPHOLOGICAL AND MOLECULAR EVOLUTION
First, we analyzed which morphological characters are variable and/or fixed within the Stygocapitella species complex. We did this to obtain character states for each lineage, and to assess the degree of variability among and across lineages. After this initial assessment, we analyzed (1) the number of segments with chaetae; and (2) the number and pattern of chaetal type in each individual chaetiger (Table S3). Only the first and the second chaetiger displayed variance in chaetal pattern, and hence the third and following segments were coded only once. Based on this information, we investigated rates of morphological evolution within Stygocapitella by means of an ancestral state reconstruction and by a regression of morphological data and pairwise genetic differences. We obtained ancestral state reconstructions by mapping morphological characters on a tree topology derived from the ML and Bayesian analyses above using Mesquite version 3.51. We applied both ML and parsimony reconstructions. In addition, we determined multidimensional morphological disparity (MMD) indices within and between species as described in Struck et al. (2017). In brief, this method relies on the decomposition of variance through a principal component analyses (PCA). Using both discrete morphological and morphometric data, we performed a PCA using the function "prcomp" included in the basic R package "stats" (R Core Team 2013). In both cases, the first four principal components were selected for the MMD calculations as they accounted for >99% of the variance. The MMD indices were plotted against the uncorrelated genetic distances of the 18S gene within and between Stygocapitella species, which we obtained using MEGA X (Kumar et al. 2018), by applying 500 bootstrap replications, the TN93 model, and a -distribution.
After this, we compared the rate of morphological evolution within Stygocapitella to other closely related groups within Orbiniida (Struck et al. 2015). To do so, we selected 12 species from Orbiniidae and 11 Nerillidae for which morphological and molecular datasets exist (Worsaae 2005;Bleidorn et al. 2009; Tables S5-S7). Importantly, both these taxa and Stygocapitella have a similar degree of genetic divergence (Struck et al., 2015). Conveniently, Orbiniidae comprises both infaunal and interstitial species, whereas Nerillidae consists exclusively of interstitial species such as Stygocapitella. The integration of the morphological records obtained herein and the records from the literature (Worsaae 2005;Bleidorn et al. 2009;Zrzavý et al. 2009;Struck et al. 2015) led to a morphological data matrix consisting of 32 species (nine from Stygocapitella, 11 from Nerillidae, and 12 from Orbiniidae) and 75 morphological characters (Tables  S5 and S6). We then conducted a PCA analysis of this dataset using the PCA option of the "FactoMineR" package (Lê et al. 2008). Based on first 18 principal components, which together explain 99.07% of the variation in the dataset, we determined MMD indices between species within Stygocapitella, Orbiniidae and Nerillidae separately, and tested if they were statistically different to each other using Tukey's HSD and pairwise Student's t-tests.
To contrast morphological and molecular evolution, we compiled a dataset of 18S sequences for each species (Table S7). Considering that 18S was the slowest evolving gene in the dataset, this gene is the ideal gene to analyze distantly related lineages. Based on this dataset, we reconstructed a ML tree using IQ-Tree (as described above) and obtained pairwise genetic distances between species. Pairwise MMD indices were plotted against the corresponding pairwise genetic distances using the "ggplot" function of the R package "ggplot2" (Wickham 2016) with a loess smoothed fit regression including confidence regions. Finally, we mapped the morphological characters on the ML tree using Mesquite version 3.51 with the ML reconstruction option. We counted the number of changes occurring at each branch, and plotted these on the ML tree to quantify the number of total changes per branch.

ECOLOGICAL DATA
Bioclimatic variables were downloaded from the world-climatic database using the "raster" R package (Hijmans 2014). Nineteen variables were downloaded with a 2.5 min of a degree resolution (21.62 m 2 at the equator) for each of the sampling sites. Because of the extensive sampling effort in the Eastern Atlantic, we subsampled sites in which only one species occurred. Ideally, we would include variables such as granularity, pH, and moisture content-but this was not possible at this stage given the requirement of specialized equipment that was not available at the sampled sites across the entire globe. Additionally, given the strong short-term fluctuations of these parameters (Giere 2009), a comprehensive dataset may only be obtained after repeated and long-term measurements over several years, accounting for events such as heavy rains and stormy weather. "Species" was treated as the dependent variable and each of the 19 climatic variables as an independent variable. Because species is a multinomial variable, we fitted a multinominal logistic regression using the R package "nnet" (Venables and Ripley 2002). After fitting each model, we performed a least squares means analysis using the R package "lsmeans" (Lenth 2013). We evaluated species-presence pairwise contrasts using F-ratios and its associated P-value.

NUMBER OF Stygocapitella SPECIES
We obtained sequences from four gene markers (i.e., 16S, COI, 18S, and ITS1) for 301 Stygocapitella subterranea, 12 Stygocapitella australis, and 18 Stygocapitella minuta specimens from 32 sites, as well as five orbiniid out-groups (Table S1). The concatenated dataset comprised 4,081 nucleotide positions. Both ML (Fig. S1) and Bayesian analyses (Fig. 3) of the concatenated molecular data resulted in the same topology. Ten well supported species were found with bootstrap values ࣙ95 and posterior probabilities ࣙ0.95 except for the Eastern Atlantic species B (EA B; bootstrap value = 87 and posterior probability = 0.81). All species, except S. subterranea, S. australis, and S. minuta, are new and formal description is pending. Relationships among species received high support with only two cases of bootstrap values <95. Single-gene phylogenies retrieved the same species as the concatenated datasets (Table S8; Figs. S1-S5). The highly conserved 18S gene was unable to unambiguously distinguish the most recent divergence between closely related species (Fig.  S2). In congruence with the phylogenetic analysis, the barcoding analyses (ABGD) found an intraspecific maximal distance of 0.031581 for COI and the recursive partitioning found eight species within S. subterranea.
The three species occurring in Europe, WA and EA C in North America, as well as EP D, EP C, and EP B in the Pacific Ocean occur in sympatry at a total of seven out of 32 sites (Figs. 1 and 2), often co-occurring within the same handful of sand. For example, three of the four Eastern Pacific species co-occur at the 4 th of July beach and all three Eastern Atlantic species are found at Musselburgh.

ANCESTRAL STATE RECONSTRUCTION
The eight discrete morphological characters present in Stygocapitella (out of a total of 75 assessed characters; Fig. 1) exhibited interindividual and interspecific variation. Ancestral character reconstructions of these discrete morphological characters suggest the occurrence of four distinct morphotypes based on the number of chaetigers and the number of specific chaetae in certain chaetigers ( Fig. 4A; Fig. S6). Morphotype #1 (red in Fig. 4) is restricted to S. minuta. Besides having only eight chaetigers, the first chaetiger of #1 has two whip-like chaetae and three bilimbate chaetae. All following chaetigers have four bilimbate chaetae. Although #1 is considerably distinct, the other three morphotypes are very similar to each other. All three (#2, #3, and #4) have 10 chaetigers and possess an additional chaetal type, forked chaetae. Morphotype #2, observed in S. australis, in specimens from 4 th of July beach and in one specimen from Plymouth in EA C consists of two whip-like, two forked, and one bilimbate chaeta in the first chaetiger and two bilimbate and two forked in the rest (green in Fig. 4; Figs. S6 and S7). In comparison to #2, #3 has one additional bilimbate chaeta in the first and second chaetigers ( Fig. 4; Fig. S6). This morphotype is observed in EA B, EA C, and WA as well as in some specimens of EA A and the Eastern Pacific. Finally, #4 has one additional bilimbate chaetae in the second chaetiger when compared with #3 and is present in the Eastern Pacific, in EA A and in one specimen from Gravesend in EA C (yellow in Fig. 4; Fig. S6). The ancestral state reconstructions unambiguously reveal that #1 is the ancestral condition for S. minuta; #2 for EP A and S. australis; #3 for EP C, EP D, and EA A; and #4 for EA B, EA C, and WA.
These results show that the Stygocapitella species complex is composed of several cryptic species as the majority of species are presently morphologically indistinguishable. Only S. minuta can be distinguished from the remaining species having a particularly distinct morphology. The remaining species display morphologies that are not species-exclusive, and differences between the three morphologies are minimal ( Fig. 4; Fig. S7-blue, green, and yellow morphotypes). The maximal difference is three additional bilimbate chaetae in a total of 44 chaetae for #4 in comparison to #2.

MORPHOLOGICAL DISPARITY WITHIN Stygocapitella
We obtained estimates of morphological disparity within the Stygocapitella cryptic complex and compared these to estimates from two closely related taxa. Mapping of the 75 morphological characters on the 18S ML tree shows that morphological evolution in Stygocapitella is slower when compared to nerillids and orbiniids (Fig. 4B). Indeed, in a total of 16 branches, only four are associated with morphological changes in Stygocapitella, while nerillids and orbiniids display changes at 19 out of 20 and 18 out of 22 branches, respectively. Within Stygocapitella, we found one branch showing three, two branches displaying two, and a single branch with one morphological change. This translates to an average of 0.5 changes per branch in Stygocapitella. In contrast, when considering branches with changes in nerillids, only six showed three or less changes, whereas 13 had between four to 12 changes, totaling to an average of 4.1 morphological changes per branch. Orbiniids displayed eight branches with one to three morphological changes and 10 branches with four to eight changes (Fig. 4B). The average number of changes per branch in orbiniids is 2.8. Thus, in Stygocapitella not only the percentage of branches with no changes is higher, but also the amount of change along a branch is considerably smaller considering both the average and the maximum number of changes.
PCA shows that Nerillidae, Orbiniidae, and Stygocapitella are clearly separated from each other (PC1 and PC2 together explain 64.4% of the variance). In addition to this, Nerillidae and Orbiniidae occupy a substantially larger morphospace area than Stygocapitella (Fig. 5A). MMD analysis shows that this difference is not restricted to the first principal components, but holds up for the first 18 principal components. MMD indices of Stygocapitella had a mean value of 1.14 with a standard deviation of 1.28 and a median of 0.41. For nerillids, the mean was 7.21 with a standard deviation of 2.76 and median of 7.86 and for orbiniids it was 7.46 with 3.29 and 7.96, respectively. Boxplots of the distributions of the MMD indices show that the quartiles of Stygocapitella do not overlap with the quartiles belonging to Nerillidae and Orbiniidae, whereas these two overlap completely (Fig. 5B). Accordingly, Tukey's honestly significant difference and Student's t-tests show that morphological disparity in Stygocapitella is significantly lower in comparison to both Nerillidae and Orbiniidae (P < 0.000001 in all cases), whereas there is no significant difference in morphological disparity between Nerill-idae and Orbiniidae (Tukey's HSD: P = 0.7343163; Student's t: Plotting of the MMD indices against pairwise genetic distances in Stygocapitella, nerillids, and orbiniids (Fig. 5C) indicates that lower MMD values in Stygocapitella are not an artefact of taxonomical ranks (Stygocapitella is a genus; Nerillidae and Orbiniidae are families). In Stygocapitella, MMD indices remain between 0 and 1 until a pairwise genetic distance of 0.025. These values increase slightly to ß4 only at higher molecular distances. This is in agreement with the mapping of morphological change along the 18S ML tree (Fig. 4B) where two branches (out of four having morphological changes) comprise five of the eight morphological changes that were reconstructed in Stygocapitella. In clear contrast, MMD indices in Nerillidae and Orbiniidae vary between 5 and 10 at relatively shallow genetic distances of <0.01 (Fig. 5C). Morphological disparity in these lineages remains at high levels with increasing genetic distances and only a few outliers display disparity values as low as Stygocapitella. These outliers are at genetic distances of 0.0375. Hence, independent of the genetic distance, morphological disparity in nerillids and orbiniids is on average 2-8 times higher than in Stygocapitella (Fig. 5C). Complementarily, we calculated MMD indices for discrete morphological traits and for the morphometric data in the Stygocapitella complex while taking the interindividual variation into account and plotted these against genetic distances in 18S (Fig. S7). Considering only MMD indices, morphotype #2, #3, and #4 cannot be separated from each other neither in discrete traits nor in morphometrics. Morphotype #1 is clearly set apart from the remaining three but only when considering discrete characters (Fig. S7A). When considering morphometric characters, morphotype #1 partially overlaps with the remaining three (Fig.  S7B). Hence, as above, the morphological similarity across different Stygocapitella species is high. Nonetheless, when considering both genetic distances and MMD together we find three clearly separated clusters. The only exceptions being morphotypes #3 and #4. Despite the very low morphological disparity, high genetic divergence can be observed between species with different morphotypes indicating that speciation in the Stygocapitella complex is not accompanied by morphological changes given the characters analyzed herein.
Finally, we plotted the results of the ancestral state reconstruction (Fig. S6) on a time-calibrated tree using the time estimates obtained from the molecular clock analysis (Fig. 3; Fig. S6). These results suggest that morphotype #3 is the ancestral condition of EA B, EA C, and WA and that these three species diverged about 18 million years (5-37 million years) ago (Fig 4A; Fig. S6). Morphotype #4 was reconstructed as the ancestral condition for a clade comprising EA A, EP C, and EP D, as well as for EA B, EA C, and WA (Fig. S6). The age of divergence for this clade (also including EP B, for which no morphological data were obtained) was estimated to be about 64 million years (33-104 million years) ago (Fig. 4A). Finally, morphotype #2 was reconstructed as the ancestral state for the whole radiation, except for S. minuta, and it was dated at ca. 140 million years (75-205 million years; Fig 4A;  Fig. S6). The age for Stygocapitella spp. was dated at about 275 million years (124-438 million years). These dates are congruent with previous estimates on a substantially smaller dataset based only on 18S (Struck et al. 2017), which calculated an age of 270 million years for the whole complex and 83 million years for the split of S. australis from S. subterranea. Hence, even though different morphotypes can be detected, long-term morphological stasis is evident. When considering the 95% confidence interval, morphotype #2 has been maintained for at least 75 million years, and #3 and #4 for at least 5 and 33 million years, respectively. To put the degree of these morphological variations into perspective, even when considering the absolute lowest estimated minimal divergence time for #2 of about 75 million years, only four morphological differences evolved in Stygocapitella, whereas almost the whole radiation of mammals took place during the same time interval.

SEPARATION
We focused on the Atlantic species for the analyses addressing potential ecological drivers of morphological stasis such as niche conservatism and/or the occurrence of fluctuating ecological dynamics (Sheldon 1996;Futuyma 2010;Lindholm 2014). However, analyses at the macro-climatic scale using 19 variables found no statistical differences in the ecological preferences within and among the species (Table S9), with P-values being generally >0.70. The lowest observed P-value was P = 0.22, for the comparison between EA B and EA C for annual mean temperature.
To investigate potential causes of stasis including niche tracking and reduction of standing genetic variation, we assessed the dispersal capacity of different Stygocapitella species. Phylogenetic reconstructions indicate that at least five trans-oceanic dispersal events (two across the Pacific Ocean and three across the Atlantic) are necessary to explain present-day distribution. This includes historical transitions, as well as modern translocations potentially due to human activity (Fig. 3, Fig. S1). Using Tajima's D, we found indications for reduction in genetic diversity in several populations (D < 0 or populations with no polymorphism; Table S10) including Bakka, Henningsvaer, Kristineberg, and Lødingen as part of EA A; Glenancross,Île Callot, Keitum, Morsum, and Nairn as part of EA B; and Bristol Channel, Ellenbogen, Hörnum, Musselburgh, and St. Efflam as part of EA C. Except for Henningsvaer (EA A), Keitum (EA B), Ellenbogen, Hörnum, Musselburgh, and St. Efflam (EA C), all other populations had significant P-values or had only a single haplotype. Finally, signatures of balancing selection or population contraction were detected for Cutty Sark (EA A), Ardtoe (EA B), Gravesend, List, and Plymouth (EA C) but none of these were significant (D > 0; Table S10).

EVOLUTION
To the best of our knowledge, this study describes the longest occurrence of stasis in cryptic species complexes known to date. It substantially exceeds results from other cryptic species complexes, which uncovered patterns of morphological stasis over a few millions of years. For instance, independent lineages of Cavernacmella snails have remained morphologically similar despite estimates of species divergence of about 3 million years ago (Wada et al. 2013) and Mastigias jellyfish have remained morphologically similar for 6 million years (Swift et al. 2016). In about 275 million years, the Stygocapitella complex evolved at least 10 reproductively isolated species with only with minor morphological differences. Reproductive isolation of these species can be indirectly inferred by complete congruence of independently evolving genetic markers, even when several species occur in sympatry (Coyne and Orr 2004), by species-specific ITS1 length (Fig. S8), and by the ABGD analysis.
Rates of morphological evolution are best described as a continuum with two ends represented by "acceleration" and "deceleration" (e.g., Stuck et al 2018). In such framework, cryptic species complexes represent cases of substantially decelerated morphological evolution. This is evident when comparing Stygocapitella to closely related taxa that demonstrate substantially faster rates of morphological evolution (Fig. 4b). On the other end of the spectrum, accelerated morphological evolution occurs in, for instance, adaptive radiations (Gillespie 2004;Simões et al. 2016) and in character displacement (Brown and Wilson 1956). Generally, increase in rates of morphological evolution seems to be connected with the availability of ecological niches (Gillespie 2004;Losos 2010). Between these two extremes, one finds cases where morphological divergence follows genetic divergence or where changes are not as pronounced as in adaptive radiations or conserved in cryptic species complexes. A continuum view will be most informative to future works in morphological evolution that should seek to understand the causes leading to shifts in the pace of morphological evolution (acceleration and deceleration), but also lineage-level morphological evolvability and constraints.
Variations in the rate of morphological change (i.e., acceleration or deceleration/stasis) are likely due to selective pressures, but also "drift in morphology." Selection can lead to the evolution of new morphologies, but it can also act to conserve a morphotype (Lynch 1990;Smith et al. 2011;Lidgard and Love 2018). Similar to molecular evolution, the absence of selection or a weak selective pressure on a particular morphotype or trait would lead to "neutral nonadaptive change" (Lynch 1990;Smith et al. 2011). Neutral nonadaptive change occurs when separate populations or species (reproductively isolated or not) accumulate differences in morphology by chance. This could result, for instance, due to random fluctuations in the mean value of a given trait or morphology or fixation of a different trait.

PUNCTUATED EQUILIBRIUM
Stasis in Stygocapitella is occasionally interrupted by pulses of quantitative and qualitative change, as predicted by punctuated equilibrium (Eldredge and Gould 1972). It is important to note that due to the lack of fossil data, and due to the few events of morphological change, we cannot infer the exact tempo and mode of morphological evolution in this group. Morphological changes could have occurred either along with speciation events as predicted by the original formulation of the punctuated equilibrium or they could have occurred in a series of subsequent changes as predicted by gradualism (Mattila and Bokma 2008;Landis and Schraiber 2017). In any case, this is in line with suggestions of punctuated equilibrium-like patterns described in biology, in molecular evolution (Pagel et al. 2006), molecular phylogenies (Bokma 2008), mammal body mass evolution (Mattila and Bokma 2008), and paleontological studies such as the armor development shifts between adaptive peaks in the three-spined stickleback that fits the punctuated equilibrium model (Hunt et al. 2008). Similar to Stygocapitella, the evolutionary patterns of other cryptic species complexes such as the Cavernacmella snails (Wada et al. 2013) and the Mastigias jelly fish (Swift et al. 2016) are best explained by stasis in combination with occasional morphological change in some lineages. In Cavernacmella, several species remain morphologically similar, whereas five became morphologically distinct. Morphological differences in these five lineages seem to be associated with the colonization of limestone outcrops (Wada et al. 2013). This colonization event led to the accumulation of differences in shell morphology in a few thousand years despite the evidence for 3 million years long stasis in C. minima lineages. This suggests that there are no constraints preventing the evolution of new shell shapes, and that shifts in the pace of morphological evolution could be linked to the ecology or habitat of Cavernacmella species, similar to adaptive radiations (Losos 2010). Similarly, in the Mastigias species complex, stasis in oceanic offshore species is repeatedly interrupted following independent colonization of marine coastal lakes. Both examples suggest that some degree of stabilizing selection on morphology occurs and when a group of individuals colonizes a new environment, there is a release of the selective pressure leading to the appearance of morphological differences. In Stygocapitella, we could not find macroecological differences within three morphologically similar species (EA B, EA C, and WA), and between them and the morphologically distinct species (EA A). However, our analyses do not include potentially relevant microecological parameters affecting the interstitial environment such as granularity, salinity, or moisture distribution (Giere, 2009).
Cryptic species strongly support the hypothesis that speciation is not necessarily coupled with morphological change (Rabosky and Adams 2012;Wada et al. 2013). In Stygocapitella, 12 out of 16 branches are not associated with morphological change in the assessed characters. In addition, the occurrence of speciation without morphological change adds a layer of complexity to the relationship among speciation, extinction, and morphological evolution (Alizon et al. 2008;Silvestro et al. 2018;Katz 2019). Failure to integrate speciation without morphological evolution in paleontology may bias estimates of rates of speciation and morphological evolution (see Alizon et al. 2008). Because in paleontology, a "species" is defined based on morphological differences (Futuyma 1987) and considering the evidence for the widespread occurrence of cryptic species (Pfenninger and Schwenk 2007; Pérez-Ponce de León and Poulin 2016), paleontological estimates of speciation and extinction could be wrong and simply reflect rates of morphological change through time. To mitigate this, we suggest that, when possible, researchers should obtain ratios of cryptic species in extant taxa when focusing on a given paleontological group. For instance, accounting for the presence of cryptic species in marine microplankton fossils demonstrates that the lumping of several cryptic species overestimates the amount of variation within morphospecies, and leads to an apparent slow-down of the rates of evolution (Alizon et al. 2008).
Although rates of anagenesis and cladogenesis have been estimated in the paleontological record (Jackson and Cheetham 1999;Aze et al. 2011;Strotz and Allen 2013), variation of these rates has received less attention. Stygocapitella has two nodes (in a total of eight) associated with more than one morphological change, a single node with a single morphological change, and five nodes without any change. This implies that at least five speciation events (i.e., cladogenesis) have occurred without any change in the characters assessed. Additionally, considering that S. australis and EP A share the same morphotype but branch off consecutively, the two subsequent morphological changes occurring in this clade had to be anagenetic. Futuyma (1987) has argued that if variation occurs in rates of anagenesis and cladogenesis, it could result in a pulsated pattern. This is seen in Stygocapitella, where we observe that anagenetic changes are lower than the rate of cladogenesis, which results in a pattern of "pulsated changes" among morphotypes.
A recent review has argued that the view on living fossils would benefit from shifting from pattern descriptions to the causes underlying consistency in morphology through time (Lidgard and Love 2018). This is similar to the recently proposed views for cryptic species . Considering that morphological deceleration occurs in both cryptic species and living fossils, and that it may be due to the same underlying forces, the integration of cryptic species and living fossils could be accomplished under a synthesis of morphological deceleration and stasis. Indeed, cases of long-term stasis in cryptic species could represent living fossils, given a paleontological record (Lidgard and Love 2018). Our results show that some Stygocapitella species have been morphologically identical in the assessed characters for at least 50 million years and if there was a fossil record for these lineages, it is very likely that extant species would have been considered living fossils.

MORPHOLOGICAL STASIS
Selection is, likely, the force that maintains highly similar/identical morphologies in Stygocapitella species. This may be due to the special properties of the interstitial realm (the space between sand grains ;Noodt 1974;Westheide 1977Westheide , 1987Giere 2009). The interstitial environment is characterized by limited and three-dimensionally structured space, and by constant changes in chemistry due to tides, wave action, seasonal changes, weather conditions, salinity, temperature, and input of organic matter. Although these factors vary extensively on short time scales (daily and seasonally), it has been shown that abiotic conditions within the interstitial realm have not changed for millions of years (Noodt 1974;Westheide 1977;Westheide and Rieger 1987;Giere 2009). This is in agreement with the "plus ça change, plus c'est la même chose" model by Sheldon (1996), which states that taxa inhabiting environments with severe short-term abiotic fluctuations, yet stable in the long term, can display stasis. In this model, species are expected to be efficient niche trackers, ultimately resulting in niche conservatism and stasis. Fluctuating conditions have been shown to occur in areas where the three European Stygocapitella species occur (i.e., Sylt, Schilksee, and Tromsø; Schmidt 1968Schmidt , 1969Schmidt , 1970Schmidt , 1972aSchmidt , 1972b. Furthermore, given their mitochondrial divergence (Figs. S4 and S5), Stygocapitella species seem to be efficient niche trackers, capable of finding new areas where their habitat occurs. This is shown by their broad distribution ranges, for example, in several populations all individuals have the same haplotype, despite the long distance between sites. Although this is a potential scenario, our attempt to test for reduction in genetic diversity, as a result of recurrent bottlenecks and founder effects (as suggested for meiofauna), did not reveal any general pattern (Andrade et al. 2011;Derycke et al. 2013). If this was the case, absence of standing genetic variation would diminish the potential for selection to act, and potentially drive morphological stasis (Futuyma 2010), yet this must be confirmed using genomelevel data. Importantly, Nerillidae, which are also interstitial, do not exhibit signs of stasis, and hence not every taxa inhabiting the interstitial environment is under stasis.

Conclusions
The Stygocapitella cryptic species complex is characterized by decelerated rates of morphological evolution and by long periods of morphological stasis. In about 275 million years, the Stygocapitella complex evolved at least 10 reproductively isolated species, but only four distinct morphotypes with few differences. This highlights that morphological evolution should be represented as a continuum (from accelerated to decelerated) and suggests that the Stygocapitella cryptic species complex is one of the most extreme examples for morphological deceleration. Even though we cannot provide conclusive evidence about the causes of stasis in Stygocapitella, stasis is likely maintained by niche conservatism coupled with the ability to track favorable habitats. On the other hand, the comparison with other interstitial species shows that the interstitial realm per se is not causing stasis. The increasing numbers of publications describing the occurrence of cryptic species suggest that speciation without morphological changes might be commonplace, and that the presence of morphologically similar species can bias paleontological rates of speciation, extinction, anagenesis, and cladogenesis.

AUTHOR CONTRIBUTIONS
JC, THS and GP have conceived this study. JC, CM, GP and THS have collected specimens in the field and obtained morphological data. JC, DSt, DSi, JW have generated genetic data. JC, DD and THS have analyzed the data. JC wrote the manuscript, with contributions from the remaining authors.

ACKNOWLEDGMENTS
The author are grateful to G. Paulay, C. E. Mills, B. Holthuis, and T. Miller for field site suggestions in the United States, and to T. Worsfold, A. Mackie, H. Reiss, and L. Jørgensen for laboratory space in the United Kingdom and Norway. The authors thank L. Thorbek for her assistance in sequencing, N. Budaeva for primer suggestions, and I. Modesto for fieldwork support in Northern Norway. The authors are grateful to P. Wagner and to three anonymous reviewers whose comments substantially improved the manuscript. The authors acknowledge the use of Norwegian national e-infrastructure for high-performance computing and storage via the projects NN9408K and NS9408K, respectively. JC is grateful to J. M. Branco ("É um lindo sonho para viver; Quando toda a gente assim quiser"), and to C. Booth and J. Lewis for late-night discussions about punctuated equilibrium. Fieldwork was partly funded by the Ragen Award from Friday Harbor Laboratories and a Den Grevelige Hjelmstjerne-Rosencroneske Stiftelse ved UiOslo (JC) and by the EU Assemble program (THS). A Forbio travel-grant led JC to Osnabrück to take Light Microscopy and SEM photographs with GP and CM. The authors thank K. Voje and L. Bachmann for reading and commenting a previous version of this manuscript. This is NHM Evolutionary Genomics lab contribution number 14.

DATA ARCHIVING
The raw genetic data are found in NCBI GenBank. The references are provided as part of the Supplementary Material. The morphological raw data are also provided as part of the supplementary data. Tree files deposited as "25302 in TreeBase.