Cyphoderia ampulla (Cyphoderiidae: Rhizaria), a tale of freshwater sailors: The causes and consequences of ecological transitions through the salinity barrier in a family of benthic protists.

The salinity barrier that separates marine and freshwater biomes is probably the most important division in biodiversity on Earth. Those organisms that successfully performed this transition had access to new ecosystems while undergoing changes in selective pressure, which often led to major shifts in diversification rates. While these transitions have been extensively investigated in animals, the tempo, mode, and outcome of crossing the salinity barrier have been scarcely studied in other eukaryotes. Here, we reconstructed the evolutionary history of the species complex Cyphoderia ampulla (Euglyphida: Cercozoa: Rhizaria) based on DNA sequences from the nuclear SSU rRNA gene and the mitochondrial cytochrome oxidase subunit I gene, obtained from publicly available environmental DNA data (GeneBank, EukBank) and isolated organisms. A tree calibrated with euglyphid fossils showed that four independent transitions towards freshwater systems occurred from the mid-Miocene onwards, coincident with important fluctuations in sea level. Ancestral trait reconstructions indicated that the whole family Cyphoderiidae had a marine origin and suggest that ancestors of the freshwater forms were euryhaline and lived in environments with fluctuating salinity. Diversification rates did not show any obvious increase concomitant with ecological transitions, but morphometric analyses indicated that species increased in size and homogenized their morphology after colonizing the new environments. This suggests adaptation to changes in selective pressure exerted by life in freshwater sediments.

INTRODUCTION

The major types of ecosystems on Earth shape the distributions and evolutionary histories of the organisms inhabiting them. Indeed, crossing the boundary between two ecosystems characterized by sharp differences in climatic (temperature, humidity, seasonality) or physiological parameters (pH, salinity) demands profound adaptations (Simpson, 1984). This explains why ecological transition events have only occurred on a few occasions in the evolutionary history of taxa. Despite its apparent permeability, the interface between marine and freshwater is one of these barriers that only relatively few groups have been able to overcome (Hutchinson, 1975). Crossing this salinity barrier implies adapting to new ionic concentrations and fundamental changes in the mechanism of osmotic pressure regulation. Accordingly, certain ancient lineages such as Echinodermata, Ctenophora (Gray, 1988) or polycystine radiolarians (Boltovskoy, 2003) seem to have never managed to achieve this transition. Other groups have been more successful; the first adaptations in plants (Mansouri et al., 2019) and animals (Hong et al., 2019; Zhang et al., 2017) can be witnessed at transcriptomic level before populations became separated and speciation occurred.

The selective pressures exerted on the organisms that overcame the salinity barrier led to population diversification and, eventually, to ecologically driven speciation in the new free niches (Rautsaw et al., 2021). These speciation events occurring during ecological transitions, together with massive extinctions, are considered as the main drivers for the disruption of diversification rates across the Tree of Life (Schluter, 2000). However, the outcome of these transitions does not follow a single rule and differs deeply among taxa. Selective pressures can vary sharply in the newly colonized environment, and act differently on the organisms depending on their own ecology and life history (Pallarés et al., 2017). For example, the effect of marine‐fresh water transitions on the rate of cladogenesis and morphological evolution (i.e., changes in phenotypic traits) varied depending on taxa, as exemplified in several bony fish families: needlefishes (Belonidae) evolved towards dwarfism (Kolmann et al., 2020) when colonizing freshwater; pufferfishes (Tetraodontidae), in contrast, increased body size heterogeneity (Santini et al., 2013) and underwent lower rates of cladogenesis; therapontid grunters (Therapontidae) increased both cladogenesis rate and morphological diversification in parallel to dietary diversification (Davis et al., 2012, 2016). Thus, there is no single scenario describing the effect of crossing the salinity barrier on the rate of diversification in fishes. However, while this question is still open to debate in relatively well‐known groups such as vertebrates, it has never been investigated in other groups of eukaryotes. There is therefore a clear a necessity to build hypotheses that may be generalized across all groups of organisms, or at least across the eukaryotes.

Eukaryotic diversity is composed to a large extent of unicellular, microbial‐sized organisms collectively called protists. This paraphyletic assemblage includes all eukaryotes with the exception of fungi, animals and plants (Taylor, 2003). Many major protist clades have been able to cross the salinity barrier, a phenomenon that is considered infrequent in plants, animals and microbes (Logares et al., 2009). The phenotypic changes brought by this transition have remained understudied because many of these protist groups have a relatively conserved morphology, with few variable traits, which make it difficult to analyse patterns of morphological diversification. One exception is testate amoebae: they possess a self‐secreted shell or test, whose morphology varies between species (Kosakyan et al., 2016). Their shape and composition have been thought to have an adaptive value to environmental factors (González‐Miguéns et al., 2021). Amongst them, family Cyphoderiidae (Rhizaria) is one of the few clades that managed to cross the salinity barrier (Heger et al., 2010). These organisms are typically found in the benthic fraction of freshwater ecosystems among sediments and aquatic plants (Todorov et al., 2009), as well as in the supralittoral zone of sandy marine ecosystems (Golemansky, 1970). The family comprises six genera with marked interspecific morphological differences (Heger et al., 2011), which makes it an excellent model for studying transitions through the salinity barrier in protists and in eukaryotes in general.

Cyphoderia is the largest genus of the family with 14 species, eleven of which have been reported from freshwater ecosystems. Cyphoderia ampulla (Ehrenberg, 1840) is the most widespread species. This freshwater species possesses tests that are typically round in cross‐section with a characteristic aperture forming an angle like a “bota bag” (Figure S1). Cyphoderia ampulla has a worldwide distribution in freshwater sediments or amongst wet peat mosses (Todorov et al., 2009). However, molecular studies showed that, despite a homogeneous shell outline and ecology, this taxon should be considered as a paraphyletic species complex (Heger et al., 2010, 2011). This homogeneity in shape and ecology contrasts with the variable morphology and ecology reported among other genera within Cyphoderiidae, and among species within the genus Cyphoderia (Todorov et al., 2009).

Here, we reconstruct the evolutionary history of Cyphoderiidae, focusing on the paraphyletic group C. ampulla, but also including its closest relatives within Cyphoderia. We based our inferences on molecular sequence data (nuclear SSU rRNA gene and mitochondrial cytochrome oxidase subunit I gene [COI]), ecological data (inferred from the habitat sampled) and morphological data obtained from the bibliography and from our own isolates from the Iberian Peninsula. In addition, we integrated environmental sequence data obtained from the EukBank public repository (Berney et al., 2017). Based on this information, we aimed at: (1) finding the original habitat of Cyphoderiidae and the direction of the ecological transitions, (2) inferring whether these ecological transitions were synchronized in geological time or occurred independently, (3) evaluating whether these transitions generated a significant change in the tempo of molecular diversification, and (4) testing if salinity barrier transitions were correlated with changes in the shell morphology, as it could be expected if shells have an adaptive value.

MATERIALS AND METHODS

Data acquisition

Sampling and specimen preparation

We sampled different freshwater sites from the Iberian Peninsula, a well‐known diversity hotspot, where we expected to have a high probability of increasing the number of lineages and the taxonomic representation of Cyphoderia ampulla. The localities sampled are detailed in Figure 1 and Table S1. The most active collected cells, those showing a high mobility under the microscope, were isolated for molecular analyses and deposited individually in Eppendorf tubes containing 100 μl of buffer for DNA extraction (Chomczynski & Sacchi, 1987). DNA was extracted from single cells as described in Duckert et al. (2018). The remaining cells were deposited on stubs for scanning electron microscopy (SEM) analyses.

FIGURE 1

Maps showing the localities from where Cyphoderiidae COI and SSU rRNA sequences were retrieved. Triangles and diamonds represent freshwater and marine environments, respectively. The colours of the symbols represent the different molecular clades

DNA extraction and amplification and microscopical observation are described in “DNA extraction and amplification” and “Microscopical observation” of Appendix 1.

Data analyses

Phylogenetic analyses

We built the COI database by combining our single cell sequences with publicly available sequences obtained from Heger et al. (2011) (Table S1), ending up with a total of 52 sequences from Cyphoderia and three from Pseudocorythion acutum, which were used to root the tree. The SSU rRNA gene database included sequences from Heger et al. (2010) and Wylezich et al. (2002), as well as environmental clones from Bass and Cavalier‐Smith (2004) (Table S1). To obtain the most complete picture of the global diversity of Cyphoderiidae, we included environmental sequences retrieved from the EukBank public repository, which contains eukaryotic metabarcodes based on the V4 region of the SSU rRNA gene, generated from about 13,000 samples distributed worldwide, involving 122 projects from over 90 research teams (Berney et al., 2017). Based on this database, we perform a blastn analyses (Altschul et al., 1990) to extract the amplicon sequence variants (ASVs) belonging to Cyphoderiidae, obtaining 64 sequences that were added to our initial data set. The tree was rooted with 26 sequences of Euglyphina (Kosakyan et al., 2016). In all, we analysed 146 SSU rRNA gene sequences; sensitivity analyses were performed with and without environmental sequences to check if the topology of the phylogenetic tree was affected (phylogenetic trees and nexus files are available at TreeBASE permalink http://purl.org/phylo/treebase/phylows/study/TB2:S29006).

Tree topologies and node supports were evaluated with both Bayesian inferences (BI) and maximum likelihood (ML) for each marker/gene individually. ML analyses were conducted using iq‐tree (Nguyen et al., 2015) and Bayesian inference (BI) analyses were conducted using mrbayes 3.2.7a (Ronquist et al., 2012). Phylogenetic analyses are described in Tree analyses of Appendix 1.

Divergence time estimation

Lineage divergence times were estimated in beast v.1.8.4 (Drummond et al., 2012), as implemented in the CIPRES Science Gateway version 3 (Miller et al., 2010) and using the SSU rRNA gene data set. We implemented the general time reversible (GTR) substitution model with gamma across‐site rate variation and an invariant portion (I + G). We calibrated lineage divergence times using two external points from sister groups belonging to Euglyphida, based on fossil evidence: (1) Trinema spp. Penard, 1890: Fossils from this genus consist in Trinema shells found in Borneo and dated from the late Miocene (Schiller, 1997), and a second fossil resembling a modern Trinema lineare from the French Pliocene (Boeuf & Gilbert, 1997). We decided to constrain the age of the most recent common ancestor (MRCA) of Trinema (the node grouping T. lineare and Trinema enchelys) with this second fossil, using a lognormal prior to reflect the uncertainty in the fossil calibration, an offset = 5 million years ago (Ma) and a standard deviation = 1. (2) Scutiglypha node: several fossils have been identified as belonging to genus Scutiglypha (represented as Euglypha). These fossils date back to the Eocene (Barber et al., 2013; Schiller, 1999) and Miocene (Foissner & Schiller, 2001; Schiller, 1997) periods, respectively. They are characterized by strongly scutiform scales without spines, which correspond with the morphology of both Scutiglypha cristata and Scutiglypha tuberculata but not with that of Scutiglypha acanthophora. Thus, we considered S. acanthophora and S. cristata as the crown group, and constrained the stem node of this group (S. acanthophora, S. cristata and S. tuberculata) with a lognormal prior with offset = 40 Ma and standard deviation = 1.

Although some Cyphoderia fossils exist (Waggoner, 1996), the preservation state of this material does not allow comparisons with modern taxa morphology. Therefore, we decided not to use these fossils as calibration points in our analyses. Sensitivity analyses were carried out to evaluate the effect of tree priors on the molecular clock. We performed several runs, with the parameters described above, and using alternative tree (Yule process, birth‐death and birth‐death incomplete sampling) and molecular clock (strict, uncorrelated, relaxed, random and fixed) priors (phylogenetic trees and nexus files are available at TreeBASE permalink http://purl.org/phylo/treebase/phylows/study/TB2:S29003). Next, marginal likelihood values were calculated for each model, using path sampling (PS) and stepping‐stone (SS) sampling (Baele et al., 2013). We then compared the model with the lowest marginal likelihood estimated, as by ps and ss (Table S2) against the remaining models in pairwise Bayes Factors comparison (Kass & Raftery, 1995) (Table S3). Finally, we also extracted the time of the most recent common ancestor (tMRCA) for the nodes of interest (freshwater clades 1 and 2, see Section 3) from the beast posterior distribution of each run using treestat (Vargas et al., 2014, 2020), in order to capture the uncertainty in divergence time estimates for those nodes.

Our sensitivity analysis indicated that the most adequate priors for our database was birth‐death (Gernhard, 2008) as tree prior and the uncorrelated lognormal distribution, for the relaxed molecular clock. We repeated three times the analysis under this model, each with two MCMC chains for 20 × 107 generations, sampling every 1000th generation, and then used tracer version 1.7.1 to monitor convergence and adequate mixing (effective sample size [ESS] >200) (Rambaut et al., 2018). We combined the three beast runs with Log Combiner, resampling states every 20,000th, with a burnin of 25%. A maximum clade credibility (MCC) tree was constructed in treeannotator version 1.8.4 (Drummond et al., 2012).

Additionally, due to the difficulties in establishing correspondences between fossils and nodes in the phylogenetic tree, we performed a second sensitivity analysis with the parameters described above but with different node positions for the calibration points, to see how these affect the estimation of nodal ages (more information in “beast sensitivity analysis” in Appendix 1).

Ancestral state reconstructions (ASR)

To reconstruct the ancestral habitat state, we used the combined MCC tree based on the SSU data set and the three beast runs, as described above. This tree was pruned to remove the outgroups, using the function drop.tip in the R package ape version 5.5 (Paradis et al., 2004). We used the phytools R package (Revell, 2012) to infer the habitat history of Cyphoderiidae. We estimated the posterior probability of the discrete characters marine and freshwater for each ancestral node, using an empirical Bayesian inference method, that is, stochastic character mapping (Huelsenbeck et al., 2003). We used the function make.simmap (Bollback, 2006) to generate 500 stochastic character maps from our data set, under two different models: all rates different (ARD) model, which allows different rates of transition between all states, and equal rates (ER), which enforces the same rate of transition between states. We then summarized these sets of stochastic maps with the function densityMap, plotting the posterior probability of being in each state across all the edges and nodes of the tree.

For the COI data set, we used the Bayesian tree, as ML and BI recovered identical topologies. We first transformed this nonultrametric consensus tree into an ultrametric one with the function chronos in ape R package (Paradis et al., 2004) and then resolved the basal polytomy in the mrbayes tree by outcome assigning a branch length of 0.00001. To reconstruct ancestral states on this tree, we followed the same procedure described above.

Diversification analyses

Our aim here was to perform an exploratory analysis of the diversification dynamics of Cyphoderiidae considering all the data we collected. To define these OTUs, we used the criterion proposed by Heger et al. (2011) to delimit species, which is based on support from molecular markers (both SSU rRNA gene and COI gene grouped conspecific sequences unambiguously together), as well as on morphological and ecological consistency. Missing species in a phylogeny can mislead estimates of diversification rates (Höhna, 2014; Sanmartín & Meseguer, 2016); although we cannot correct our analyses for possible missing species, as the number of lineages considered as independent taxonomic units or OTUs in our analysis is considerably higher than any of the estimates given in the literature or based on Cyphoderiidae systematics. However, given the size of the database explored, we can reasonably assume that we have a good representation of the groups diversity. We thus performed our diversification analyses assuming that we have retrieved all the diversity in Cyphoderiidae, that is., assuming a sampling fraction of 1.

We used the MCC tree described above, without outgroups, and further removed the intraspecific variability by pruning the tree to leave only one single terminal per species (74 tips). We then used the function ltt.plot in the package ape version 5.5 (Paradis et al., 2004) to plot the number of lineage diversification events through time. These lineage‐through‐time (LTT) plots were drawn for each of the individual beast runs, the combined MCC tree, and two subtrees representing the trajectories of the freshwater lineages, and the marine lineages, separately. For each tree, we also estimated the gamma statistics (Pybus & Harvey, 2000), and assessed its significance by comparing the observed value against a distribution of values from 1000 random trees with 74 tips, generated under a Yule model with the empirical speciation rate estimated in ape.

We used a reversible jump Metropolis‐coupled Markov Chain Monte Carlo (MCMCMC) approach implemented in the Bayesian analysis of macroevolutionary mixtures (bamm version 2.5); Rabosky, 2014; Rabosky, Donnellan, et al., 2014), following the bamm tutorials (Mitchell & Rabosky, 2017), to model the dynamics of speciation and extinction over time and across branches in the combined MCC tree above (with every tip representing a single species). The prior block parameter values were chosen using bammtools package (Rabosky, Grundler, et al., 2014), with a Poisson rate prior of 1.0. We ran four chains with 10 × 106 generations, sampling every 1000th, the final posterior distribution was estimated from 1000 trees. MCMCMC operators were used under default settings. bamm has been criticized recently because of its biased modelling of the probability of rate shifts occurring in extinct lineages (i.e., a lineage undergoes a rate shift and then goes extinct), leading in some cases to incorrect likelihood estimation (Laudanno et al., 2021; Moore et al., 2016). We used here an alternative lineage‐specific birth‐death shift model (lsbds), which uses the incremental‐time, state‐dependent differential equations developed by Maddison et al. (2015) to model shifts in diversification rates dynamically (Höhna et al., 2019), and which is implemented in the Bayesian software revbayes version 1.1.1 (Höhna et al., 2014; Höhna, Landis, et al., 2016). We followed the tutorial (https://revbayes.github.io/tutorials/divrate/branch_specific.html), using four rate categories (i.e., allowing for three rate shifts in the speciation rate), ρ = 1 (global sampling fraction), and a chain length of 10 × 106 generations. We then obtained the branch‐specific diversification rate estimates for each branch of the tree with a rejection‐free stochastic rate mapping algorithm and dividing the phylogeny into 500 time slices (Freyman & Höhna, 2019). We also performed two additional analyses to explore the effect of assuming varying levels of missing species in the phylogeny: using a global sampling fraction of ρ = .5 and ρ = .25.

Both bamm and lsbds model changes in diversification rates over time that are lineage‐dependent. To explore the possibility of a rate shift that affected all lineages in the phylogeny simultaneously (e.g., a mass extinction event), we ran an episodic birth‐death model to estimate the timing and magnitude of discrete shifts in diversification and relative extinction rates. We used the compound Poisson process of mass extinction (comet) model (May et al., 2016) implemented in the R package tess (Höhna, May, et al., 2016), with a log‐normal empiricalHyperPriorForm and sampling probability = 1, a chain length of 20 × 106 iterations, sampling every 100th generation, and a minimum number of effective sample size of 1000 to stop the run. comet uses independent compound Poisson process priors to estimate the posterior probability of the timing and magnitude of rate shifts in speciation and extinction; mass extinction events that deplete the standing diversity by a certain fraction (ρ) can also be estimated by integrating out speciation and extinction rate shifts via Bayes factor comparisons (May et al., 2016). Additionally, we ran an episodic birth death model that is implemented in revbayes and uses auto correlated Brownian rates across time slices, and a horseshoe estimator for sparse signals (Magee et al., 2020). We followed the tutorial (https://revbayes.github.io/tutorials/divrate/ebd.html#Hoehna2015a) with ρ = 1 and a chain length of 10 × 106 generations. We used the R package revgadgets version 1.0 (Tribble et al., 2021) to visualize the rates. For both comet and revbayes, we performed two additional analyses changing the global sampling fraction to ρ = .5 and ρ = .25, considering that we only sampled the 50% and 25% of the total diversity of Cyphoderiidae, respectively.

To explore if diversification rate shifts across branches detected by bamm and lsbds (results) could be explained by transitions between the marine and fresh‐water habitat, we performed state‐dependent extinction and speciation diversification rate analyses. We first ran binary state speciation and extinction (BiSSE) (Maddison et al., 2007) analysis using the R package diversitree (FitzJohn, 2012) and two states: (1) marine and (2) freshwater. We performed empirical Bayesian estimation of the parameters using the MCMC option. As our tree only included 74 tips, we used an exponential prior with a rate of 1/(2r) (where r is the character independent diversification rate) as recommended in the diversitree tutorial; the step size (argument w) was set as in Revell and Betancur (2017). We performed two runs of 30,000 steps, with a final burnin of 20% for all runs. We also repeated the BiSSE analysis in revbayes (Höhna, Landis, et al., 2016), which implements a fully Bayesian approach, following the procedures in (López‐Estrada et al., 2021) and (https://revbayes.github.io/tutorials/morph/morph_more.html) (Figure 3c). We ran the analysis with a log‐standard deviation of = 0.587405 and a chain length of 40,000 generations; we used 500 time slices to estimate the timing and number of state changes along branches using stochastic character mapping. The BiSSE model has been criticized because it assumes a null model for statistical testing, the constant‐rate birth death process (i.e., heterogeneity in diversification rates occurs virtually in any phylogeny; Rabosky & Goldberg, 2015, 2017). To test for the existence of diversification rate shifts that are unrelated to state transitions in the focal trait (marine‐freshwater habitats), we compared four different SSE models (explained in Appendix 1) in the R package HiSSE (Beaulieu & O’Meara, 2016) using the Akaike information criterion (AIC).

FIGURE 3

(a) Time‐calibrated maximum clade credibility (MCC) tree based on the SSU rRNA gene data set of family Cyphoderiidae. Branches are proportional to absolute time, with the x‐axis representing geological epoch boundaries. Branch colours and pie charts at nodes indicate marginal probabilities for ancestral state reconstructions and stochastic mappings for marine (red) and freshwater (blue) habitats. Lineage through time plots isare shown for different MCC trees: the thin green lines represent each of the independent beast runs for the entire phylogeny; graygrey line represents the combined MCC tree; red line represents the marine lineages only, and the blue line, the freshwater lineages only. (b) Marginal densities for the age of the most recent common ancestor (tMRCA) for the freshwater clade 1 (orange) and the freshwater clade 2 (purple). A schematic illustration of sea level fluctuations over geological epochs (Eberli, 2000) is shown on the top. The three green lines represent the marginal densities of the tMRCA of each independent beast run. (c) Posterior distribution of speciation and transition rates for marine and freshwater habitats, and transition rates between obtained them, with the BiSSE model in diversitree and revbayes

Morphometric analyses

We used five continuous morphometric traits for the morphometric analyses, as described in Todorov et al. (2009): (1) Length of the shell, (2) width of the shell, (3) large pseudostome axis, (4) length‐large pseudostome axis ratio, and (5) length‐width ratio. We combined measurements from Todorov et al. (2009) with our new data from the Iberian Peninsula, in total 719 cells (the data are available in Morphobank permalink http://morphobank.org/permalink/?P4082). We followed the same groups created by Todorov et al. (2009), as the monophyly of each of these lineages was supported by both SSU rRNA gene and COI gene sequences (Heger et al., 2011). Our cells were ascribed to these groups based on their COI and SSU rRNA gene phylogenetic position.

We carried out a principal component analysis (PCA) to observe the relationships between the different phylogroups (lineages recovered monophyletic with COI and SSU), using the stats package, with the prcomp function with the five continuous morphometric traits described above. Finally, we plotted the results using the package ggplot2 (Wickham, 2016) and used the function factoextra (Kassambara & Mundt, 2020) to see the contribution in the analysis of each variable used.

Hypervolume

We calculated the overlap between each phylogroup using the R package hypervolume (Blonder et al., 2014), allowing the quantification of the phylogroups morphospace as determined by the measured of five continuous morphometric traits: shell length, shell width, large pseudostome axis, length‐large pseudostome axis ratio and length‐width ratio. We used the first three principal components from the previously described PCA, to obtain a limited number of linearly uncorrelated variables. We used a Silverman bandwidth estimator (Silverman, 2018) with a quantile threshold of 0.05 to infer the Kernel density estimations (KDEs). Then, we created a multidimensional space for each phylogroup, calculating its size and the extent of overlap with other phylogroups, using the function hypervolume_overlap_statistics and hypervolume_distance, respectively (Blonder et al., 2014).

Combining phylomorphospace and hypervolume

We tested how far ecological transitions are correlated with changes in the morphology by mapping the morphological diversification history of the freshwater and marine clades, respectively. As not all measured specimens yielded a DNA barcode, but all morphology groups were recovered at least with two monophyletic sequences (both COI and SSU), we consider here each lineage‐morphology group as an independent evolutionary unit. We therefore pruned the SSU rRNA MCC, leaving only one tip per morphology group to link the molecular and the morphometric data. Then, we took the centroids of the hypervolume described above (Table S4) and mapped the morphological diversification history of each lineage by plotting the tree tips in a two‐dimensional space. We reconstructed the ancestral states along morphospace axes (hypervolume centroid), using the function phylomorphospace in the phytools package (Revell, 2012; Sidlauskas, 2008). The hyperphylomorphospace resulting from this combination of approaches allows the evaluation of both lineage morphological similarity (distance between each lineage) and the trajectory of the phenotype (directionality of tree branches). The procedure followed is described in Appendix 1.

Morphological similarity and equivalence: Testing the possible morphological adaptive value

Finally, we wanted to test how far the shell morphology of the different species was a response to selective pressures exerted by the new environments found after the transition across the salinity barrier. Assuming that there is a selective pressure of the environment (freshwater) on shell morphology, we would expect to see a homogenization of the shells. If this was the case, morphologies would be more similar between freshwater clades than with their phylogenetic sister marine clades. Thus, to statistically test the morphological patterns observed in Cyphoderia, we calculated two indexes: (1) “morpho‐niche” equivalence: whether the morphologies of resembling species are effectively indistinguishable; and (2) “morpho‐niche” similarity, whether the morphology model in one species predict the morphology of a second species better than expected by chance (Warren et al., 2008).

We first defined a background morphology, a “potential morphology” that can represent each morphology group, by drawing a hypervolume (as described above) with five continuous morphometric traits (see morphometric analyses). As the hypervolume generates random points, inferred by a Silverman bandwidth estimator (described above), we considered these estimated points as the potential morphological background, which were used to calibrate our variables (aka the native area in a niche equivalence test). To do this calibration, we performed the PCA‐env (Broennimann et al., 2012) with the function dudi.pca, using the random points of the hypervolume and the measures of the organism. Then, we performed niche similarity and equivalence tests (Warren et al., 2008), in the R package ecospat (Broennimann et al., 2018), using the functions ecospat.niche.similarity.test and ecospat.niche.equivalency.test, respectively, with 1000 replications in both tests. We also plotted the Niche overlap between each freshwater phylogroup with the function ecospat.plot.niche.dyn.

RESULTS

Molecular phylogenetic and divergence time estimation analyses

The nomenclature of the clades (freshwater clades 1 and 2) and the taxonomic identifications of the species were given after Heger et al. (2011). The COI tree included 52 sequences (Figure 2; Table S1). Our phylogenetic analyses recovered the monophyly of Cyphoderia with a Bayesian posterior probability (PP) of 1 and a maximum likelihood bootstrap (ML) of 100. The freshwater clades 1, 2 and the Cyphoderia major (Penard, 1890) clade were recovered with a PP = 1/ML = 91, PP = 1/ML = 96 and PP = 1/ML = 95, respectively. The freshwater clade 1 grouped six mitochondrial lineages (1–6 orange colour), the Cyphoderia major clade, two lineages (7 and 8 green colour), and the freshwater clade 2, four lineages (9–12 violet colour), all supported by PP = 1 and ML = 95. Cyphoderia compressa (Golemansky, 1979) was recovered as sister to freshwater clade 1 with a PP = 0.94 and ML = 66; the addition of sequences from our new isolates is not congruent with the topology in Heger et al. (2011), where C. compressa is sister to freshwater clade 2 (PP = 0.54 and ML = 60 values). Cyphoderia compressa appears as paraphyletic in the COI tree, which strongly suggests that both mitochondrial lineages are distinct species (and potentially a novel clade not recovered by SSU).

FIGURE 2

Bayesian phylogenetic tree based on COI sequences of family Cyphoderiidae. Posterior probability values (Bayesian analysis) and bootstrap values (maximum likelihood analysis) are represented at each node. Pie charts at nodes represent marginal probabilities for ancestral state reconstructions obtained in phytools (red colour for marine environment, blue for freshwater). Branch colours represent the molecular clades. GenBank accession codes are given on the tips. Species names correspond to the taxonomic assignation of the organisms based on morphology for each mitochondrial lineage; indices correspond to illustrations. Symbols represent the original habitat of the organism. The images on the right are scanning electron microscopy photographs showing the shell morphology at the same size scale, together with details of the scaling patterns of each phylogroup. The map represents the localities sampled for each clade

For the SSU rRNA, we generated 146 sequences, in which 120 represent Cyphoderiidae (Figure S2; Table S1). In our phylogenetic analyses, we recovered the monophyly of Cyphoderiidae and Cyphoderia (freshwater 1, 2, marine 1, 2, C. compressa, C. littoralis and C. major clades) with a PP = 1/ML = 100 and PP = 1/ML = 100, respectively. The freshwater clade 1, 2 and Cyphoderia major clade were recovered with a PP = 1/ML = 100, PP = 1/ML = 100 and PP = 1/ML = 100, respectively. Sequences from the EukBank database grouped into a new freshwater clade (3) with a PP = 1/ML = 100, coming from a single geographic location (Lake Sanabria, Spain). The freshwater clade 1 is composed of four nuclear lineages, Cyphoderia major clade by one, and freshwater clade 2 by two all supported by PP = 1 and ML <85. Cyphoderia compressa is the sister group to freshwater clade 1 with a PP =0.99 and ML =98; recovering the same topology as in our COI tree. EukBank sequences also highlighted the existence of new clades, only present in marine environments (marine clades 1, 2, 3; Figure 3a and Figure S2).

The maximum clade credibility (MCC) tree obtained from beast (Figure 3a, Figures S3–S5) recovered the same clades as the nonclock trees obtained with iq‐tree and mrbayes (Figure S2). The most recent common ancestor (MRCA) of Cyphoderiidae was estimated to have been living during the Mesozoic, around 121 million years ago (Ma) in the Early Cretaceous (95% highest posterior density [HPD] = 57.8–223.2 Ma). This node is where the largest discrepancy between sensitivity analyses was found (Figures S3–S5). More fossil data, if possible, from Cyphoderiidae themselves, will be needed to resolve more accurately the age of this clade. However, the nodes including both freshwater clades 1 and 2 (and Cyphoderia ampulla) were more homogeneous among sensitivity analyses (Figure S4). They support the existence of a MRCA in both clades during the Miocene. The MCC tree indicates that the MRCA of the freshwater clade 1 probably occurred around 14.2 Ma (95% HPD = 5.7–28.4 Ma), and for the freshwater clade 2, at 6.3 Ma (95% HPD = 1.6–14.6 Ma) obtained similar results with other models (Figure S6).

Ancestral state reconstruction and diversification analyses

Ancestral state reconstruction results obtained with phytools for the COI and SSU trees were congruent in recovering that the most probable habitat for the common ancestor of family Cyphoderiidae and genus Cyphoderia was marine (Figures 2 and 3a and Figure S7). The larger SSU tree recovered four branches in which the ecological transitions occurred from marine to freshwater, and we did not detect any transition in the opposite direction (Figure 3a). The lineages of these four freshwater clades (Cyphoderia major clade, freshwater clades 1, 2 and 3), are named in the following analyses as “freshwater lineages”, for short.

The LTT plots of the combined MCC tree (Figure 3a) show a nonsignificant positive gamma (γ) value (γ = 0.5469, p = .5845) indicating that the speciation rate was constant over time; while the random Yule trees showed a significant negative gamma (γ = –0.1490, p = .033). The LTTs for the marine (γ = 0.2625, p = .793) and freshwater (γ = 0.7702, p = .4412) lineages gave also a nonsignificant positive γ (Pybus & Harvey, 2000).

bamm and revbayes branch‐specific diversification rate analyses gave congruent results (Figure 4a,b), showing an increase in the diversification rates in the clade formed by genus Cyphoderia and the marine clade 3. When experimentally decreasing the sampling fraction (ρ = .5 and ρ = .25), the overall diversification rate increased but branch diversification patterns were similar to those obtained with ρ = 1 in each of the two analyses (Figure 4 and Figure S8).

FIGURE 4

Diversification analyses: (a) revbayes: Estimated branch‐specific speciation rates in the pruned SSU rRNA tree, each tip representing an independent taxonomic unit; circles at nodes are marginal probabilities for ancestral character states with BiSSE (red: marine; blue: freshwater; branch colorscolours present marginal stochastic character mapping); the net‐diversification graph shows results from the episodic birth‐death model implemented in revbayes (autocorrelated Brownian model). (b) bamm: Estimated branch‐specific speciation rates in the pruned SSU rRNA tree; the net‐diversification graph shows results from the episodic birth‐death model in comet (compound poissonPoisson prior model). (c) Estimates of speciation and extinction rates from the episodic birth‐death model implemented in revbayes. (d) Estimates of speciation and extinction rates, and relative extinction rates, from the episodic birth‐death model implemented in comet; the extinction plot depicts also the probability over time of mass extinction events; significance was assessed using Bayes Factor comparisons (BF >2, positive support Kass & Raftery, 1995).

The episodic birth‐death model showed differences between the revbayes and comet analyses (Figure 4 and Figure S8). revbayes recovered an upward trend in the net diversification rate, which is mainly explained by an initially high extinction rate (between 120 and110 Ma), and a gradual, slow increase in the speciation rate (Figure 4c). comet inferred more constant net‐diversification rates over time ending with a rapid decrease around 5 Ma (Figure 4b), caused by a significant downward shift in the speciation rate (Figure 4d, BF > 2). There were no significant shifts in the extinction rate; however, the comet analysis supported a mass extinction event between 120 and 110 Ma in the Bayes factor comparisons (BF ~2), in agreement with the revbayes signal (Figure 4c). Experimentally decreasing the sampling fraction (ρ = .5 and ρ = .25) recovered patterns that were very similar to those obtained with the full sampling by the two analyses, showing an overall increase in the net diversification rate (Figure S8).

BiSSE with Diversitree showed that both speciation and extinction rates are higher in freshwater lineages (λ = 0.07702517; µ = 0.03678705) than in marine lineages (λ = 0.0458489; µ = 0.01513946) (Figure 3c). BiSSE with revbayes found similar results: (λ = 0.06372; µ = 0.02799) in freshwater and (λ = 0.0466; µ = 0.01666) in marine lineages (Figure 3c). The speciation rate was higher in the freshwater lineages than in the marine lineages in 92.31% of the post‐burn posterior distribution, suggesting that speciation rate increases when transitions occur but that these differences are not significant (Figure 3c). “Null” was the best model to explain the observed diversification patterns, according to AIC comparisons (Table S5). This indicates that the “habitat” trait marine‐freshwater does not explain the diversification rate heterogeneity among clades in the phylogeny.

Morphometric and hyperphylomorphospace analyses

We delimited 10 morphogroups, or independent taxonomic units (Heger et al., 2010; Todorov et al., 2009), for the morphometrics analysis. These morphogroups were defined, based on the molecular results, as: Cyphoderia littoralis (Golemansky, 1973) (dark red in Figure 5), Cyphoderia compressa (light red in Figure 5), Cyphoderia major (green in Figure 5), freshwater clade 1 (orange in Figure 5) and 2 (violet in Figure 5).

FIGURE 5

(a) Bayesian phylogenetic tree based on the SSU rRNA data set, pruned from Figure 2, leaving one tip per species; branch colours represent the habitat (red, marine; blue, freshwater), (b) Principal component analysis of five continuous morphological traits. (c) Hyperphylomorphospace of the phylogeny (a) with the combinations of the three principal components (b); symbols represent the hypervolume centroids of each phylogroup, and the lines, the Kernel density of the random points generated in the hypervolume

The first two axes in the PCA (PC1 = 65.96%; PC2 = 25.87%) explained 91.83% of the total variation of the morphological data set (Figure 5b). The variables with the highest contribution were shell width and length (Figure S9A). Both PCA and hypervolume plots showed that C. littoralis, C. compressa and C. major groups have specific morphospaces that do not overlap (Table S6). In turn, freshwater clades 1 and 2 do overlap, with their centroids close to each other (Table S6). This morphological resemblance contrasts with their phylogenetic position, indicating that these two groups are distantly related (Figure 5c).

The first two axes (axis1 = 63.24%; axis2 = 24.88%) of the PCA‐env analysis, used to calculate morphological similarity and equivalence, explains 88.12% of the morphological data set variation (Figure S9B). The Schoener's D metric presents similar values with the hypervolume overlap (Tables S6 and S7), showing overlap only between freshwater clades 1 and 2 (Table S7). Cyphoderia trochus palustris (Penard, 1902), Cyphoderia ampulla 2 and C. ampulla 5 (freshwater clade 1) present a morphology more similar than expected by chance with C. ampulla 11 (freshwater clade 2) p‐value > .05) (Warren et al., 2008), and equivalent between C. ampulla 2 and C. ampulla 11 p‐value > .005 (Figure 6). This also contrasts with their respective phylogenetic positions (Figure 5c).

FIGURE 6

Top: composite represents the morphology overlaps between the phylogroups inferred from COI and SSU results (names in the diagonal); the x‐ and y‐axes correspond to PC1 and PC2, respectively, and the green colour represents the overlap between each group. Bottom: morphology overlap comparisons between each morphogroups in base to the p‐values obtained in the morphology equivalency/similarity (representing only A to B groups) tests

DISCUSSION

Marine origins of family Cyphoderiidae and Cyphoderia, and their parallel preadaptation to freshwater

Ancestral state reconstruction analyses inferred the marine habitat as the most likely original habitat for both genus Cyphoderia and family Cyphoderiidae. This result is supported by ASR analyses on both mitochondrial (COI) (Figure 2 and Figure S7) and nuclear (SSU rRNA gene) (Figure 3a) markers, and by state‐dependent models (Figure 4). In line with this finding, our survey of the EukBank environmental SSU rRNA V4 sequence data set revealed that the majority of Cyphoderiidae amplicon sequence variants (ASVs) stored there come from marine coastal sediments (Figure 1 and Figure S2). Furthermore, all sequences from deep environmental clade come from these environments (Figures 2 and 3 and Figure S7). This stands in contrast with the fact that most species described in Cyphoderia (14) are from freshwater environments, while only three species (Cyphoderia littoralis, Cyphoderia compressa, and Cyphoderia truncata [Schulze, 1875]) have been described from the marine environment. This knowledge gap demands more research on testate amoebae from marine ecosystems, a surprising fact given that current hypotheses based on fossils and phylogenomic reconstructions support a marine origin for the group (Porter et al., 2003; Porter & Knoll, 2000).

Most marine Cyphoderiidae specimens have been discovered in shallow areas near the coast (Figure 1; Bass & Cavalier‐Smith, 2004). Some observations describe their presence in the supralittoral zone of sandy beaches, and especially abundant several metres inland from the coastline (Golemansky, 2007). Interestingly, in our study, the Cyphoderia compressa clade, a euryhaline group, was reconstructed as sister to freshwater clade 1, while marine clade 1 (sequences retrieved from sandy beaches, probably also euryhaline) appeared as sister to the freshwater species Cyphoderia major (Figures 2 and 3a). This suggests that euryhaline environments possibly hosted the ancestors of each of these freshwater cyphoderiid clades.

The euryhaline environment is exposed to quick changes in salinity (Golemansky, 2007), and fluctuating osmotic pressures. These fluctuations imply physiological stresses to which marine organisms need to be adapted. Euryhaline environments possibly acted as a transition environment in the preadaptation of marine organisms to freshwater ecosystems, as these transitions, and vice versa, from freshwater to the marine habitat, are known to require more than one transitional step (Coates, 1993; Lee et al., 2012). Therefore, many species that crossed the salinity barrier seem to have descended from euryhaline ancestors that lived in naturally fluctuating environments; examples of this evolutionary pathway can be found in bony fishes, such as freshwater terapontid grunters (Davis et al., 2012), killifishes (Whitehead, 2010), or in the three‐spined stickleback (Raeymaekers et al., 2005). Other freshwater animals exhibit similar patterns (Lee & Bell, 1999), as well as some clades of eukaryotes, such as cryptomonads (Shalchian‐Tabrizi et al., 2008), dinoflagellates (Logares et al., 2007), diatoms (Alverson, 2014; Alverson et al., 2007), and even plants (Dittami et al., 2017). These dynamic euryhaline ecosystems have sometimes facilitated the reverse colonization of marine systems by freshwater organisms (for example, in some groups of insects such as the marine caddisflies, Riek, 1977), but this event is rare (Lee & Bell, 1999). For example, our analysis, based on a large data set of sequences retrieved from the EukBank repository, did not reconstruct any instance of recolonization of the marine environment by freshwater Cyphoderiidae (Figure 4). This asymmetry in the salinity barrier transition is common in eukaryotes (Betancur et al., 2015), and has been explained by a decreased probability of metabolic adaptation to the higher osmotic pressures that characterize the marine environment, or by competitive exclusion of freshwater species by marine resident communities (Bloom & Lovejoy, 2012; Whitehead, 2010). In Cyphoderia, one explanation is that the lack of salinity fluctuations in freshwater ecosystems prevented the colonization of euryhaline zones, due to either exclusive competition with marine clades or the loss of key adaptations to high salinity environments, or both.

From a common marine ancestor, our phylogenetic reconstruction indicates that there have been at least four independent events of transition from marine to freshwater environments in Cyphoderiidae: in the Cyphoderia major clade and in the freshwater clades 1, 2 and 3 (Figure 2a). This suggests that those key innovations that facilitated adaptation to freshwater environments in the family evolved independently from distantly related and, presumably, euryhaline ancestors. This potential case of parallel preadaptive evolution (Anker et al., 2006) provides strong evidence for a key adaptation in the evolutionary history of Cyphoderiidae (Futuyma, 1998), which could facilitate further adaptive radiation (Mayr, 1960). Similar parallel ecoevolutionary pathways have been found in other eukaryotic groups, such as fishes (Whitehead, 2010) and crustaceans (Mamos et al., 2016), as well as in some groups of protists like diatoms (Alverson et al., 2007) and dinoflagellates (Logares et al., 2007).

Sea level fluctuations favour ecological transitions through the salinity barrier

Ecological transitions through the salinity barrier can occur at any time, but they are more likely to happen when marine incursions (transgressions) expand the contact zones between freshwater systems and the sea, often isolating landlocked populations in zones with fluctuating salinity (euryhaline) (Schultz & McCormick, 2012). Fluctuations in the sea level may drive and structure species distributions, connectivity between populations, and can also drive speciation‐extinction patterns. As described above, the creation of euryhaline zones during marine transgressions provided an opportunity for marine organisms to colonize freshwater environments. However, these saltwater incursions could also have triggered extinction events in other, already‐adapted freshwater lineages, opening ecological opportunities for new colonizers from the marine and euryhaline environments to diversify via competition release; the latter has been often posited as a driver of adaptive radiation (Losos, 2010; Yoder et al., 2010) in marine fishes (Betancur et al., 2012). For this reason, the glacial‐interglacial cycles of the Pleistocene have been considered as key periods in the diversification of continental organisms. Eustatic changes in the sea level (Lambeck et al., 2002) during this epoch probably facilitated the ecological transition from marine to freshwater ecosystems in crustaceans (Dooh et al., 2006), and, in general, led to accelerated diversification rates in freshwater species (Hewitt, 2000).

Miocene ecosystems underwent even higher sea level fluctuations than the Pleistocene (Kominz et al., 2008; Westerhold et al., 2005), and this geological epoch was probably one of accelerated diversification rates (Miller et al., 2005). Indeed, our divergence time estimation analysis, and the sensitivity analysis exploring the uncertainty in lineage divergence times, suggest that ecological transitions through the salinity barrier in Cyphoderiidae date back to the mid‐late Miocene (Figure 3a, Figure S6), a period of massive sea level decreases. We postulate here that the cyclic fluctuations of sea level during the mid‐late Miocene generated new niches for marine Cyphoderiidae in periods of marine transgression (rise in sea level), allowing them to colonize freshwater ecosystems (via the euryhaline transition habitat) and favoured by the extinction of freshwater competitors. Subsequent speciation would have taken place in periods of marine regression (decrease in sea level) via competitive exclusion, thus generating parallel ecological transitions cyclically. The correspondence between Miocene sea‐level fluctuations, transitions through the salinity barrier, and diversification have been extensively studied in fishes (Bloom & Lovejoy, 2017; Lovejoy et al., 2006; Tsigenopoulos et al., 2003). This pattern has also been recurrently observed in other eukaryotes, such as amphipods (Mamos et al., 2016; Yang et al., 2013), molluscs (Vermeij & Wesselingh, 2002), plants (Triest & Sierens, 2014) or diatoms (Hayashi et al., 2018). Again, cyphoderiids seem to follow evolutionary patterns that had been already described in macroscopic organisms, suggesting similarities in ecological transition pathways.

It has been demonstrated that diversification rates increase substantially with every transition through the salinity barrier. This event has been mostly described in clades of fishes (Bloom et al., 2013) and diatoms (Alverson, 2014). The explanation invoked is the lower connectivity of freshwater systems in comparison to the sea, along with the higher variability in abiotic factors. In our analysis, both state‐dependent models (BiSSE in diversitree and revbayes) and the gamma test results support higher rates of cladogenesis in freshwater than marine lineages (Figure 3). However, a better fit to the null model (Table S5) suggests that the heterogeneity in diversification rates observed in Cyphoderiidae cannot be explained by transitions towards the freshwater habitat. Instead, the branch‐specific diversification rate models (bamm, revbayes BDS) suggest that speciation rates may be clade‐dependent rather than state‐dependent: for example, the group formed by Cyphoderia and Marine clade 3 exhibit higher diversification rates than other marine or freshwater clades in Cyphoderiidae, such as Freshwater clade 3 or Marine clade 4 (Figure 4).

The most recent common ancestor of genus Cyphoderia was dated at c. 50 Ma, during the Eocene, and the episodic birth‐death models in comet and revbayes (Figure 4) show a gradual increase in the diversification rates from the Early Cenozoic onwards (Figure 4). Silica is a limiting factor in aquatic systems, and its increased bioavailability during the Eocene has been linked to diversification events in organisms that use this resource to build mineral structures (Lazarus et al., 2009; Rabosky & Sorhannus, 2009). It is possible that the higher diversification rate exhibited by genus Cyphoderia and the Marine clade 3 was favoured by the ecological opportunity offered by the increase in silica bioavailability during this epoch. Competition between cyphoderiids and other benthic groups (e.g., diatoms) for the water dissolved silica might have caused an acceleration of speciation events in these clades; similar hypotheses have been posited for planktonic groups that construct siliceous structures, such as radiolarians or diatoms (Falkowski et al., 2004; Katz et al., 2004; Sims et al., 2006). Furthermore, we observed a decrease in speciation rates in Cyphoderiidae starting during the Pliocene that is temporally congruent with similar patterns in diatoms and radiolarians (Katz et al., 2004; Lazarus et al., 2009). Though other explanations cannot be discarded (lack of resolution in the selected molecular markers, incomplete taxon sampling, etc), the observed similarity in diversification patterns between such phylogenetically distant groups as testae ameba, diatoms, and radiolarians, suggest a clade‐wide effect of abiotic factors (long‐term environmental changes, geological events), acting at the biome level (Ezard et al., 2011). Further studies will help determine the drivers behind these patterns.

Concomitant changes in morphology and adaptive value of the shell

Mineral elements have played a major role in the shaping of evolutionary patterns in many groups of microeukaryotes, such as coccolithophorids, diatoms and foraminiferans (Katz et al., 2004). When comparing freshwater and marine Cyphoderiidae, Golemansky (2007) noticed the larger size of shells in the freshwater forms, and hypothesized that the small and compressed shells of marine species were an adaptation to life between sand grains. Our study confirms that all surveyed freshwater species are larger than the marine isolates (Figures 5 and 6). All transitions inferred in Cyphoderiidae occurred from the marine environment towards the freshwater habitat (Figure 3). We observed that, alongside a general increase in the shell size, there was a homogenization of the shape after the transition, that is, freshwater forms became ballooned or expanded in comparison with marine forms (Figure 5, PC1 axis). Interestingly, other testate amoebae living in marine supralittoral systems (e.g., Arcellinida) are smaller than their freshwater counterparts (Golemansky & Todorov, 2004).

We thus hypothesize, reversing Golemansky's hypothesis, that the increase in size in freshwater species compared with marine lineages of Cyphoderiidae was an adaptive process acting on the freshwater lineages, and working directly or indirectly through the shell morphology. Adaptive selection in response to similar selective pressures may also explain the polyphyletic nature of the complex Cyphoderia ampulla sensu lato (specimens named C. ampulla in Figure 2; and those in Freshwater clades 1 and 2 in Figure 3). Our phylomorphospace and principal component analysis show that Freshwater clades 1 and 2 are morphologically more similar to each other than to other phylogenetically closer marine species (i.e., C. compressa clade, Figures 5 and 6). However, while size and shape are conserved between specimens with C. ampulla morphology, the shape, size and distribution of the siliceous scales that form their shells are highly variable (Husnot, 1943). This morphological diversification of the shell scales was probably correlated with molecular diversification. Heger et al. (2011) noticed that Cyphoderia ampulla clades could be easily discriminated based on the overlapping or nonoverlapping of the scales, a finding that remains valid after adding new isolate sequences (Figure 2). In Euglyphida, in general, scale size has been shown to be highly variable and probably species‐specific (Chatelain et al., 2013; Wylezich et al., 2002). Altogether, this suggests that scale shape and size are less prone to selective pressures than shell shape and size, and should therefore constitute better traits for clade delimitation in Cyphoderiidae systematics.

CONCLUSIONS

Cyphoderiidae is one of the few groups of benthic marine protists that successfully managed the transition from marine to freshwater ecosystems. These transitions are inferred to have occurred during the mid‐late Miocene, a period of frequent fluctuations in the sea level. The ancestors of present‐day freshwater species were most probably marine coastal organisms that colonized freshwater systems with the euryhaline environment as a transitional step. This probably took place during marine transgressions, associated to ecological release, which was in turn favoured by the extinction of already‐adapted freshwater clades, and the ecological opportunity offered by the increasing levels of water‐dissolved silica. Finally, our study confirms that microscopic organisms such as benthic protists, may exhibit similar evolutionary patterns to those observed in macroscopic organisms, supporting the idea that global patterns apply to all life on Earth, no matter the organisms’ size.

AUTHOR CONTRIBUTIONS

Rubén González‐Miguéns and Enrique Lara designed the study; Rubén González‐Miguéns, Carmen Soler‐Zamora, Andrés Blanco‐Rotea, María Isabel Carrasco‐Braganza, David de Salvador‐Velasco, Antonio Guillén‐Oterino, and Daniel Tenorio‐Rodríguez carried out the fieldwork; Rubén González‐Miguéns and Carmen Soler‐Zamora performed the laboratory work; Cédric Berney analysed the EukBank data set; Carmen Soler‐Zamora prepared the images; Rubén González‐Miguéns and Isabel Sanmartín performed the diversification analyses and subsequent interpretation of the results; Rubén González‐Miguéns performed the rest of the analyses, interpreted and discussed the data together with Enrique Lara; Rubén González‐Miguéns, Carmen Soler‐Zamora, Fernando Useros, Sandra Nogal‐Prata, Cédric Berney, David Velázquez, Thierry J. Heger, Isabel Sanmartín and Enrique Lara wrote the manuscript.

CONFLICT OF INTEREST

The authors have no conflicts of interest.

BENEFIT‐SHARING STATEMENT

Benefits generated: There are no benefits outlined in the Nagoya protocol associated with this study to report. All data produced in this research are deposited in public databases as described above.

Supplementary Material

Click here for additional data file.