Contributions of biogeographical functions to species accumulation may change over time in refugial regions

Elevated biodiversity is the result of the cradle, museum or sink functions. The contributions of these three functions to species accumulation and their changes through time remain unknown for glacial refugia. Additionally, our understanding of the role these functions played during pre‐glacial periods is limited. We test for changes in contributions of functions through time leading to the current diversity patterns using a model refugium and taxon.


| INTRODUC TI ON
Glacial refugia are of particular interest to biogeographers, evolutionary biologists and ecologists who are trying to understand the heterogeneous distribution of species. Their often high species richness can be the result of three possible biogeographical functions: (a) high in situ speciation rates -the "cradle function", (b) high immigration rates -the "sink function" and (c) range contraction resulting in relict species -the "museum function" (Pulliam, 1988;Stebbins, 1974;Stenseth, 1984; Figure 1).
Nevertheless, the relative contribution of the three functions to species accumulation and their changes through time remains speculative for refugia. So too does the understanding of the role pre-glacial periods may have played in shaping extant diversity patterns of refugia.
Inferring the underlying processes requires a suitable model taxon that meets the following requirements -It should be widespread and species-rich to test the functional associations of species occurring in the refugia (i.e. a result of either the cradle, museum or sink function). Also, the taxon's evolution should extend over a considerable geological time to assess the change of these functions through time. A candidate taxon that fulfils the above-mentioned requirements is the aquatic snail genus Theodoxus Montfort, 1810.
Theodoxus may have originated in central Europe as early as the Palaeocene (Bandel, 2001;Kowalke, 2002). Its diversification has been suggested to be closely linked to both older drainage fragmentation and recent glaciations (Bandel, 2001;Bunje, 2005;Bunje & Lindberg, 2007). Today, the genus has a mostly temperate distribution across Europe, south-western Asia and northern Africa, where representatives inhabit freshwater to mesohaline aquatic systems (Bandel, 2001). Based on a model refugium (Anatolia) and a model taxon (Theodoxus), we test if there is indeed a change of contributions of functions through time leading to the current diversity pattern. Furthermore, using these models, we assess if pre-glacial periods may have played greater roles in biodiversity accumulation for regions considered glacial refugia than previously thought. To achieve this objective, our specific goals were to: (a) estimate ages and infer accumulation of lineages through time using a dated phylogeny based on both available fossil data and substitution rates, (b) estimate ancestral areas through biogeographical models subjected to the dated phylogeny and (c) use a set of operational criteria to define the contributions of the three functions through time (see Figure 1). As species assignment within Theodoxus remains controversially discussed and thus elevated species diversity in refugia is largely unsubstantiated, we used preliminary analyses to identify molecular operational taxonomic units (MOTUs) that could serve as phylogenetic species (through a combination of phylogenetic and barcoding approaches). In addition, we determine if this diversity is elevated in glacial refugia (through interpolated mapping and phylogenetically based calculations of species richness). accumulation through a temporal perspective to adequately explain present-day biodiversity patterns.

K E Y W O R D S
Anatolia, biodiversity accumulation, cradle function, museum function, sink function,

Theodoxus, Western Palaearctic
As a working hypothesis, we expect the sink and cradle functions to dominate in pre-glacial periods in Anatolia, given the long and complex limnological history and dispersal gateways with other parts of the WP, but to contribute relatively little to the accumulation of interspecific diversity in the region. In contrast, the sink and museum functions may dominate during the pronounced glacial cycles reflecting the retreat of species to southern WP refugia and/or range contractions of widespread species, and contribute more substantially to the interspecific diversity of the genus in Anatolia. By inferring the relative contributions of the three functions and their changes through time, we hope to improve our understanding of how these processes may have shaped biodiversity in similar glacial refugia worldwide. This study may also highlight the roles pre-glacial conditions have played in biodiversity accumulation and emphasize that contemporary patterns of diversity need to be considered in a temporal context.

| Sampling design and morphospecies identification
Three hundred and forty-eight specimens (343 Theodoxus, 4 Neritina, 1 Neritilia) were hand-collected from hard substrate in shallow waters or dredging off small boats in deeper water ( Figure 2). Specimens were preserved in 80% ethanol, before being photographed and utilized for DNA extraction. Emphasis was placed on collecting as many Theodoxus morphospecies from as many locations across the geographical range of the genus as possible. Collected specimens were identified based on regional species lists, original species descriptions and by consulting taxonomic specialists (see Appendix S1, Table S1.1). Most Theodoxus-based molecular studies to date either closely mirrored our own sample localities or lacked nuclear DNA (nDNA) data (Bunje, 2005(Bunje, , 2007Bunje & Lindberg, 2007;Fehér, Szabó, Bozsó, & Pénzes, 2009a;Fehér, Zettler, Bozsó, & Szabó, 2009b). Therefore, only an additional 34 (33 Theodoxus, 1 Neritina) carefully chosen published sequences, representing potentially unique species or populations, were included into analyses to avoid a maternal bias and oversaturation of our analyses with redundant data (Figure 2; Table S1.1). In all instances, special note was made of the location for mapping and biogeographical analyses (Table S1.1). In total, we were able to obtain 32 of the 37 commonly recognized Theodoxus morphospecies. The five missing species have either never been recorded alive and are thus impossible to incorporate in a phylogenetic framework (T. gloeri; Odabaşi & Arslan, 2015) or efforts to gain samples were unsuccessful (T. doriae, T. mesopotamicus, T. pallidus and T. schultzii).

| DNA extraction, amplification, sequencing, saturation and model selection
Foot tissue was used to extract total genomic DNA for each specimen via a DNeasy Blood and Tissue kit (QIAGEN). DNA sequences for two mitochondrial (mtDNA) gene fragments, cytochrome c oxidase subunit I (COI) and 16S rRNA (16S), and one nDNA intron fragment, ATP synthetase subunit α (ATPα), were then amplified (see Table S2 for primers and PCR conditions). Finally, purification of gene fragments and bidirectional Sanger sequencing were carried out by LGC Ltd (Berlin, Germany).
The three gene fragment datasets (COI, 16S and ATPα) showed no significant degree of saturation (p < 0.001). Thereafter, jMod-elTest 2.1.10 (Darriba, Taboada, Doallo, & Posada, 2012) and the Bayesian information criterion (BIC) (Schwarz, 1978) was used to determine the general best-fit model of sequence evolution for each dataset (Posada, Buckley, & Thorne, 2004). The HKY+I+Γ model was selected for COI and 16S, and the HKY+Γ model was selected for ATPα. F I G U R E 1 Operational criteria based on ancestral distributions used to identify and distinguish between the cradle, museum and sink functions. The functions are identifiable by: (a) the cradle function, having an in situ distributed ancestor and descendant with either strict in situ distribution or both in-and ex situ distribution, (b) the museum function, having an ancestor with both in-and ex situ distribution followed by a descendant with only in situ distribution and (c) the sink function, having an ex situ distributed ancestor with either a strict in situ distribution or with both in-and ex situ distribution [Colour figure can be viewed at wileyonlinelibrary.com]

| Determination of interspecific molecular diversity
Provided barcoding gaps are present between different monophyletic phylogroups containing a single morphospecies, a barcoding approach can be a useful tool to distinguish between different levels of diversity and define MOTUs (representing interspecific diversity) among paraphyletic species or phylogroups where complexes exist.
A COI barcoding approach has already been shown to be useful for defining phylogenetic species within other Neritidae (Chee & Mohd Nor, 2014;Frey & Vermeij, 2008).
To identify initial MOTUs of supported monophyletic phylogroups of individual Theodoxus morphospecies, we constructed two combined phylogenies of all genes (1,608 characters; 549 parsimony informative) based on unweighted maximum parsimony (MP) and model-informed Bayesian inference (BI). All analyses were conducted through the free online resource CIPRES Science Gateway (Miller, Pfeiffer, & Schwartz, 2010). The MP analysis was performed with PAUP 4.0b (Swofford, 2002), using the heuristic search option, with TBR branch swapping and 100 random taxon additions. In instances where multiple equally parsimonious trees were retrieved, only the best tree was saved during each replicate. The robustness of nodes were assessed by 10,000 bootstrap (BS) replicates and values >70% were considered supported (Felsenstein, 1985). The BI analysis was performed using MrBayes 3.2.6 (Ronquist et al., 2012) to determine the posterior probabilities (PP) of associations. Two parallel Markov Chain Monte Carlo (MCMC) simulations used five chains for 50,000,000 generations, saving one tree in every 2,000 generations. The first 50% of trees were discarded as burn-in as assessed by parameter convergence in Tracer 1.6 (Rambaut, Suchard, Xie, & Drummond, 2014) after standard deviation (SD) of split frequencies had reached stationarity.
To distinguish and substantiate between intra-and interspecific diversity we used MEGA 7.0.26 (Kumar, Stecher, & Tamura, 2016) to calculate uncorrected genetic p-distances for the barcoding gene COI. Mean values between and within initial supported MOTUs of singular monophyletic Theodoxus morphospecies were used as reference for intra-and interspecific diversity. These were then applied to paraphyletic assemblages or phylogroups where species complexes existed. In instances where mean values were between the thresholds of intra-and interspecific diversity attained for monophyletic morphospecies, the nearest supported node (PP ≥ 0.95 and/ or BS ≥ 70%) between thresholds encompassing a monophyletic entity was used to equate intraspecific diversity and define the remaining MOTUs.

| Diversity across the landscape
To infer if elevated interspecific diversity occurs in glacial refugia, a heat map indicating the spatial variation in diversity of Theodoxus MOTUs across the WP was produced in ArcGIS 10.4 (ESRI, 2016) following Neubauer et al. (2015) in order to correct for oversampling. The MOTUs, defined by molecular analyses and COI barcoding, were used to plot observed interspecific diversity. An Inverse Distance Weighting (IDW) interpolation surface was computed on the number of MOTUs per grid cell, using a cell size of 200 km 2 , a power of two and variable search radius with 12 nearest input sample points as settings. Where areas of concentrated diversity were found, phylogenetic diversity (PD; Faith, 1992) and sum of evolutionary distinctiveness ( sum Ed; Isaac, Turvey, Collen, Waterman, & Baillie, 2007;Redding, 2003) were calculated as indicators of the relative amounts of diversity thereof. For this, we used the R package Picante 1.6-2 (Kembel et al., 2010) for the R statistical environ-  Table S1.1 for more information [Colour figure can be viewed at wileyonlinelibrary.com] & Müller, 2010) from which most of the specimens were excluded, resulting in a species tree-like phylogeny of the 18 interspecific MOTUs) and a community matrix of presence-absence data were used.

| Molecular clock and ancestral distribution analyses
Spatiotemporal distributions of common ancestors of in situ MOTUs is needed to distinguish contributions of the cradle, museum and sink functions through time (Figure 1). To do this, BEAST 1.8.4 (Drummond, Suchard, Xie, & Rambaut, 2012) was used for constructing a multilocus dated molecular phylogeny. Input files were generated in BEAUti 1.8.4 (Drummond et al., 2012) Kowalke, 2002). For congruency, five independent MCMC simulations ran for 200,000,000 generations, sampling every 10,000 generations.
Validation of convergence and mixing was assessed in Tracer 1.6 to ensure that all effective sample size (ESS) values were >200. LogCombiner 1.8.4 (BEAST package) was used to combine log and tree files applying a 75% burn-in. Trees were subsequently summarized in TreeAnnotator 1.8.4 (BEAST package) with no additional burn-in.
To place the geographical perspective on phylogenetic divergence events, the R package BioGeoBEARS 0.2.1 (Matzke, 2013a) was used to estimate ancestral areas under six different biogeographical models for the entire dated phylogeny. This included the DEC, DIVA-like and BayArea-like models including the +J parameter for each that allows for founder-event speciation (for details on the differences between these models, see Matzke, 2013bMatzke, , 2014. Here, a simplified analysis with default settings in which only two areas were predefined, namely "Anatolia" and "non-Anatolia" was run, and the best-fit model was determined by using the Akaike information criterion (AIC) as implemented in BioGeoBEARS.
Improving the model by constraining the ancestral areas of lineages using fossil evidence is hampered by the paucity of records for Anatolia (but for examples, see Alçiçek, Wesselingh, & Alçiçek, 2015;Wesselingh, Alçiçek, & Magyar, 2008) and the general difficulties in assigning fossil Theodoxus species to lineages based on morphology alone.
To infer if the majority of diversity is attributed to specific points in time, two lineage-through-time (LTT) plot analyses were performed using the R packages Ape 4.1 (Paradis, Claude, & Strimmer, 2004) and Phytools 0.5-64 (Revell, 2012). In a first attempt, the complete BEAST MCC tree (excluding the outgroups) was subjected to the LTT plot analysis including a 95% confidence interval based on the posterior distribution that was generated in LogCombiner 1.8.4 (Drummond et al., 2012). A second LTT plot analysis was conducted using the same simplified tree with 18 specimens as representatives of each relevant MOTUs (see Section 2.4). For clarity, we defined clades and subclades based on dates. Clades represented lineages present at the end of the Miocene, while subclades represented lineages present at the end of the Pliocene. The uncorrected COI barcoding p-distances among these initial MOTUs of monophyletic morphospecies ranged from 13.81% (SD ± 0.33%; between MOTUs A and M) to 2.17% (SD ± 0.00%; between MOTUs I and J). On average, interspecific diversity was found to be 9.17% (SD ± 2.40%). In contrast, uncorrected p-distances within these MOTUs ranged from 0.00% (SD ± 0.00%; for MOTU J) to 2.30% (SD ± 1.32%; for MOTU N -see note in Table S1.3 for an explanation of the use of this group as a monophyletic morphospecies), with an average of 0.97% (SD ± 1.12%; Table S1.3). Given the barcoding gap, we considered mean uncorrected p-distances for phylogroups of >2.50% to represent interspecific and <2.00% intraspecific diversity. Where uncorrected p-distance values fell between the values of intra-and interspecific diversity, phylogenetic support was used to determine the level of diversity (see Section 2.3).

| Interspecific diversity
Through barcoding and phylogenetic support, the remaining Theodoxus spp. were found to occur across seven MOTUs (H, K, N, O, P, Q and R; Figure 3a). In total, 18 interspecific MOTUs were identified within Theodoxus.

| Diversity across the landscape
Inverse Distance Weighting interpolation mapping characterised five areas with some level of concentrated diversity across the WP: Anatolia, the Balkans, Central Europe, Iberia-Mauretania and Sarmatia (Figure 4). Anatolia is the most diverse region thereof considering both local concentration and the overall amount of diversity followed by the Balkans, Iberia-Mauretania, Central Europe and then Sarmatia (Figure 4).

| Molecular dating and ancestral distributions
Three clades and 12 subclades were defined based on geological periods (Figures 3a and 5). Molecular clock analyses, substantiated with MP and BI nodal support, suggest that Theodoxus likely  (Table S1.4). It shows that the most recent common ancestor   (Figures 4 and 5). Through these initial results, our biogeographical and temporal analyses revealed a complex pattern over time of the interplay between functions. Below follows a chronology of the accumulation of lineages, their ancestral distributions and thus the changing contributions of the cradle, museum and sink functions over time with the potential drivers thereof.

| Miocene-Pliocene transition
Two divergence events occurred within Theodoxus during the Miocene (Table 1; Figure 3a). The corresponding MRCAs had wide distributions that likely included both Anatolia and other parts of the WP (Figure 5; Table S1.4). Later, the MRCA of clade 3 had an Anatolian-specific distribution before it begun to diversify in the Early Pliocene (see C3; Figure 5). As such, Anatolia probably acted as museum during the Miocene-Pliocene transition (6.78-4.77 Ma; HPD 9.40-3.26 Ma). The museum function is the only function that can definitively be associated to Anatolian Theodoxus over this time period. The HPD of the two Miocene diversification events in Theodoxus aligns relatively well with the Messinian Salinity Crisis (MSC) (Krijgsman, Hilgen, Raffi, Sierro, & Wilson, 1999;Roveri et al., 2014). During the MSC, a positive freshwater budget in the Black Sea 5.5-5.33 Ma allowed freshwater to enter the Mediterranean Basin, known as the "Lago Mare" event (Gliozzi, Ceci, Grossi, & Ligios, 2007;Stoica, Krijgsman, Fortuin, & Gliozzi, 2016;van Baak et al., 2015). This event may have facilitated short-range dispersal between adjacent areas in the eastern Balkans and western Anatolia for fresh-and brackish water taxa (Gliozzi et al., 2007;Levy, Doadrio, & Almada, 2009). A number of fresh-and brackish water lakes existed during this time in Anatolia and may have acted as suitable environments for Theodoxus (e.g. Denizli, Söke, Karacasu, Çameli, Eşen, Baklan, Acigöl, Burdur, Tuz and Çankırı: Alçiçek, 2010;Gürbüz & Kazancı, 2015;Yavuz, Culha, Demirer, & Aydın, 2017). At the end of the MSC, in the Early Pliocene, Anatolia became more isolated due to the rising sea level in the Mediterranean (Karakitsios et al., 2017;Popov et al., 2006;de la Vara, van Baak, Marzocchi, Grothe, & Meijer, 2016). The web of ancient drainage basins within Anatolia would have, for a time, provided opportunities for the persistence of previously widespread species (i.e. relicts) before they diversified into independent lineages.

| Pliocene
During the Pliocene, the amount of lineages rapidly rose from 3 to 12 (Figure 3a,b). Of these further nine diversification events, six were strictly associated to Anatolia (Table 1; Figure 5). By the end of the Pliocene, 5 of the current 12 Theodoxus lineages were present in Anatolia ( Figure 5). All Pliocene Anatolian-associated descendants derived from Anatolian ancestors ( Figure 5). As a result, the entire build-up of Theodoxus diversity in Anatolia during this period was the result of in situ diversification (i.e. cradle function). The Anatolian Pliocene climate is considered to have been warmer and more humid than today (Alçiçek, 2010). Increased precipitation sustained a high diversity of freshwater environments until the earliest Pleistocene (Alçiçek, 2010). Additionally, ongoing tectonic activity at that time strongly affected the physiography of some of the basins (Brocard, Meijers, Willenbring, Kaymakci, & Whitney, 2015; Gürbüz & Kazancı, 2015;Özsayin et al., 2014). For example, basin fragmentation probably led to divergence differences in the mollusc assemblages of two adjacent basins of Denizli and Baklan (Wesselingh & Alçiçek, 2010). Similarly, vicariance as a result of drainage fragmentation and isolation over this time period was suggested as the main driver of F I G U R E 3 (a) Dated phylogeny of the genus Theodoxus constructed in BEAST based on COI, 16S and ATPα. Node labels denote divergence times in millions of years ago (Ma); node bars indicate the 95% credibility interval around these dates. Small squares at nodes indicate significant support of divergence events found with BEAST and other phylogenetic analyses (see Figures S2.1 and S2.2), as explained through the key. Where MOTUs (A-R) show conspecifics among a number of morphospecies, species names are given in order of their year of description. Morphospecies, incorporated from GenBank, where determination was potentially dubious are highlighted by an asterisk. Clades (C) and subclades (SC) are demarcated by dashed lines between MOTUs. (b) LTT plots indicating the build-up of lineages in Theodoxus over geological time. Dashed lines surrounding the solid LTT lines indicate the 95% confidence intervals. Where intra-and interspecific diversity diverge, interspecific diversity is highlighted in blue and intraspecific diversity in red. Transitions in geological ages are highlighted by narrow grey lines, while the grey bar marks the period of pronounced glacial cycles (last 900 kyr) [Colour figure can be viewed at wileyonlinelibrary.com] speciation in spined loaches and chub in Anatolia (Bohlen, Perdices, Doadrio, & Economidis, 2006;Durand et al., 2000). Although it remains difficult to differentiate the individual drivers of divergence events for each lineage, we hypothesize that in situ diversification in Anatolian Theodoxus was promoted by changes in topography and drainage evolution.

| Early Pleistocene to Middle Pleistocene
Six phylogenetic divergence events occurred during this period, which increased the total number of Theodoxus lineages from 12 to 18 (Figures 3a,b and 5). The number of lineages with a presence in Anatolia rose from five to eight ( Figure 5). Similar as for the Pliocene, Anatolian-associated descendants derived from MRCAs restricted to Anatolia ( Figure 5). Again, the entire build-up of Theodoxus interspecific diversity in Anatolia during this period was the result of in situ diversification (i.e. cradle function). Although the topography of Anatolia was still changing as a consequence of tectonic activity (Westaway et al., 2005), diversification may have also been driven by the result of early Quaternary glacial cycles. The glacial cycles that preceded the Middle Pleistocene were more frequent, yet less severe (Elderfield et al., 2012). They may not have been prominent enough to drive the retreat of species into refugia (i.e. sink function), but it is possible that the shifting climates would have driven habitat heterogeneity (Alçiçek, 2010;Alçiçek & Jiménez-Moreno, 2013), promoting further ecological speciation, such as in tree frogs (Gvoždík et al., 2010). Moreover, there is evidence suggesting that some of the Pliocene Anatolian lakes developed into wetlands and river systems by the end of the Gelasian (Alçiçek, 2010).

| Middle Pleistocene to present-day
Our results suggest a lack of interspecific diversification events during this period (Figure 3a,b). However, three MOTUs colonized Anatolia between the Middle Pleistocene and present-day increasing the Anatolian diversity from 8 to 11 lineages ( Figure 5). The additional MOTUs (K, Q and R) are all suggested to have arisen elsewhere in the WP ( Figure 5; Table S1.4). Their immigration into Anatolia is a clear indication of the sink function, which appears to be the only function contributing to the biodiversity build-up during this time.
Around 900 ka, glacial cycles became more pronounced leading to c. 100-kyr-long glacial cycles with massive fluctuations in temperature, precipitation and global sea levels (Elderfield et al., 2012). The low sea levels during glacial maxima (Deuser, 1972;Frigola et al., 2012;Lambeck, Rouby, Purcell, Sun, & Sambridge, 2014) caused the Black and Marmara seas to become cut-off from the Mediterranean Sea and develop brackish to freshwater conditions (Aloisi et al., 2015;Schrader, 1978;Verleye, Mertens, Louwye, & Arz, 2009). The lowered salinities would have facilitated the retreat of WP groups further south as well as promote exchange with the Balkan refugia.

| CON CLUS IONS
Our working hypothesis -sink and cradle functions dominating in pre-glacial times while contributing little to the accumulation of interspecific diversity versus sink and museum functions dominating in glacial times and contributing considerably to interspecific diversity -is partially rejected. Our results, indeed, suggest that these functions did not operate simultaneously in the accumulation of species richness in Anatolia. Specifically, the cradle and museum functions alternately acted during pre-glacial times where most interspecific diversity was accumulated, while during the glacial period only the sink function acted and contributed comparably little to diversity buildup. The fact that most of the extant lineages were present in Anatolia before the onset of the pronounced glacial cycles was somewhat unexpected given the available dispersal pathways to other parts of the WP during this time (e.g. the freshening of the Black Sea; Badertscher et al., 2011) and the possibility that pathways used to explain recent range expansion of a number of aquatic species out of Anatolia could be mirrored for immigration (Akın et al., 2010;Fehér et al., 2009a,b for examples of recent range expansions out of Anatolia). This minimal contribution of the sink and museum functions does not mitigate the role glacial cycles played. Perhaps, the Ice Ages were less effective in forcing temperate aquatic species to retreat into more opportune areas than previously thought. In turn, elevated diversity in refugia may rather be the result of early in situ speciation. These results highlight the importance of considering species accumulation through a temporal perspective to adequately explain present-day patterns of biodiversity. As this case study is based on a single model genus and region, we encourage future research to expand both taxonomically and spatially to test for the changing roles of cradle, museum and sink functions. High species diversity in refugia that has been mainly attributed to glacially forced retreats may need to be re-evaluated for regions that have a complex and long geological history.