Phylogeny, ancestral ranges and reclassification of sand dollars

Classification of the Class Echinoidea is under significant revision in light of emerging molecular phylogenetic evidence. In particular, the sister-group relationships within the superorder Luminacea (Echinoidea: Irregularia) have been considerably updated. However, the placement of many families remains largely unresolved due to a series of incongruent evidence obtained from morphological, paleontological, and genetic data for the majority of extant representatives. In this study, we investigated the phylogenetic relationships of 25 taxa, belonging to eleven luminacean families. We proposed three new superfamilies: Astriclypeoidea, Mellitoidea, and Taiwanasteroidea (including Dendrasteridae, Taiwanasteridae, Scutellidae, and Echinarachniidae), instead of the currently recognized superfamily Scutelloidea Gray, 1825. In light of the new data obtained from ten additional species, the historical biogeography reconstructed shows that the tropical western Pacific and eastern Indian Oceans are the cradle for early sand dollar diversification. Hothouse conditions during the late Cretaceous and early Paleogene were coupled with diversification events of major clades of sand dollars. We also demonstrate that Taiwan fauna can play a key role in terms of understanding the major Cenozoic migration and dispersal events in the evolutionary history of Luminacea.

The irregular echinoids of the order "Clypeasteroida" sensu A. Agassiz 1872-1874 or "sand dollars" sensu lato, constitute a morphologically well-defined group of burrowing sea urchins living on sandy bottoms from the intertidal to the deep-sea 1-3 .Because of their characteristic flat body, they are commonly known as sand dollars or sea biscuits and were previously subdivided into two suborders "Scutellina" and "Clypeasterina", respectively.However, the monophyly of this traditionally defined morpho-group has been challenged since the advent of molecular studies [4][5][6][7][8][9] .Members of the former suborder "Scutellina" were resolved as sister group of irregular echinoids of the suborders Cassiduloida + Echinolampadoida 7 , and together, they were sister to the "Clypeasterina" 7,8 .Based on molecular evidence, Mongiardino Koch et al. 8 proposed an elevation of the suborder "Scutellina" to "Scutelloida" and restriction of usage of the order "Clypeasteroid" solely to member of the former suborder "Clypeasterina".To date, four modern irregular echinoid clades, Clypeasteroida, Scutelloida, Cassiduloida, and Echinolampadoida are recognized; these taxa together constitute the clade Luminacea 10 .Luminacea first appeared in the Middle to Late Jurassic 7 .As bioturbators, they provide significant services to the marine ecosystem.Dead Luminacea also contribute an important part of shallow marine sediments (e.g., in form of sand dollar mass deposits 11 ) and have indirectly influenced the ecosystem functions responding to global changes, in modern time as well as in the past.

Discussion
Based on the most dense taxon sampling to date, comprising nine out of ten families and 18 out of 29 extant genera of sand dollars sensu lato, our phylogenetic results reject the sister-group relationship between Clypeasteroida and Scutelloida, which is congruent with previous DNA-based studies 4,5,7,8 .The monophyly of the currently used superfamily Scutelloidea Gray, 1825 17 (including the families Astriclypeidae, Dendrasteridae, Mellitidae, Scutellidae) as well as the Scutelloida 7,8,10 in other usage should be re-examined.Among the Scutelloida, the sistergroup relationship between Laganiformes and Scutelliformes was often suggested by previous studies, yet the phylogenetic inferences were either based on fewer representative genera 4,5,7,8 or unbalanced character-sampling of a molecular vs. morphology data in combined analyses 7,13 .This relationship is rejected by the present study.We argue that the attributed morphological similarity in these two groups could be the result of convergent evolution.
The phylogenetic position of the family Taiwanasteridae, represented by Sinaechinocyamus mai, has been reviewed recently 9,13,21,33,34 .It is closely related to Dendraster, Echinarachnius, and Scaphechinus (Fig. 2).Among them, the eccentric apical disc is an adaptation to the feeding strategy of Dendraster and its feeding posture 35 .Other than that, the four groups are similar-non-lunulate, with similar petals (Figs.S3, S4).Sinaechinocyamus is a relatively long ranging genus with good fossil records since the late Miocene (~ 8 Ma) in Taiwan 36 .Furthermore, www.nature.com/scientificreports/new occurrences of Sinaechinocyamus have been found in coastal regions in Japan (Fig. S5), extending the geographic distribution outside of Taiwan Strait.
The new superfamily Taiwanasteroidea is proposed here to include families of Dendrasteridae, Echinarachniidae, Scutellidae, and Taiwanasteridae.Descriptions of Taiwanasteroidea together with other two new superfamilies Astriclypeoidea and Mellitoidea are given at the end of this section.
With the new phylogeny of the Luminacea reconstructed, we demonstrate important evidence on solving the mystery of lunule origins.Seilacher 25 stated that lunule evolved at least six times in sand dollars.Combining molecular and fossil evidence, Mongiardino Koch and Thompson 7 proposed a lunulate clade consisting of Astriclypeus and Mellita (as well as other fossil forms).In our study, lunule belong to two distinct clades (Astriclypeoidea and Mellitoidea) separated by the non-lunulate clade (Taiwanasteroidea) (Figs. 2 and 3).This indicates independent origins for the lunule in the Luminacea.The formation of lunule in Mellitidae consists of plates with a festooned arrangement and that is significantly different from the ones in Astriclypeidae with a cross-linked arrangement (see ref. 25 ).In a broader sense, the macroevolutionary trend of Echinoidea suggested on the basis of molecular evidence 9 similarly reflects the "Dollo's law of irreversibility" 37 , a hypothetical scenario showing: (A) innovation of lantern in early echinoids during the late Paleozoic; (B) loss of lantern at adult stage among early irregular echinoids during the Mesozoic; and (C) independent re-innovation of the lantern with modifications 33 in adult clypeasteroids and scutelliformes during the Cenozoic.
With new biogeographic analyses presented here (Fig. 4A-C), hypotheses 12,18,25,38 on the origin(s) of the modern groups of sand dollars can be tested.Smith 12 suggested that modern sand dollars have three biogeographical radiation hotspots to explain their present-day pattern of diversification.First, the diverse Arachnoidinae (including the genera Arachnoides and Fellaster) within the Clypeasteroida originated in the Australian region.Second, the ancestral Scutelloida emerged in the Caribbean Sea/ Mediterranean Sea.Third, the ancestral Rotulidae arose in tropical West Africa.Based on the abundant, solid fossil evidence, Seilacher 38 further noted that Taiwan could be another radiation hotspot where modern Astriclypeidae emerged.There are 11 geographic provinces defined in this study (Fig. 4A).As such, migration can be well illustrated with the biogeography network visualization software StrainHub 39 .The tropical eastern Indian and western Pacific Ocean (EIWP) have the highest source/hub ratio (SHR).Tropical eastern Pacific (EP) and northwestern Pacific (NWP), and tropical western Atlantic and Caribbean Sea (WA) have moderate SHR (Fig. 4B).Other nodes have low SHR.In addition, we reconstruct the polarity and frequency of migration events between geographic provinces.EIWP to NWP and EP to WA have a high frequency of migration.EIWP to EP has moderate frequency of migration.The other edges show a low frequency of migration.
Results show that EIWP is one of the key biodiversity and migration centers for shallow marine echinoids (Fig. 4B).In addition to Australia that has been identified previously 12,25 , the region around Taiwan is an important migration juncture for key echinoids mentioned here (Fig. 4C).Smith 12 hypothesized that the clade of Astriclypeidae originated in the Mediterranean Sea based on their common presence in the fossil record of the region.Based on our analyses (Figs. 3 and 4B), however, the clade of modern Astriclypeidae appears to have a reversed migration trend, originating from Asia and then spreading to the Red Sea (Fig. 4C).
The MRCA of Taiwanasteroidea likely started to expand its range from tropical origins to high latitudes along the Pacific coast with warming ocean currents during hothouse conditions.This superfamily gradually gave rise to a cold-water lineage following their range expansion to the northeastern Pacific and northwestern Atlantic during subsequent periods of cooling.A similar model of migration can also be observed in other invertebrates.The range expansions of gastropod superfamilies Turritelloidea and Buccinoidea in the northern Pacific have mainly been influenced by thermal deterioration 40 .The expansions of boreal Buccinidae, Beringiinae, and Turritellinae were caused by the progressive cooling beginning in the Late Eocene 40 .In the genus Littorina (Gastropoda: Littorinidae), speciation occurred in response to climatic cooling during the Cenozoic at higher latitudes on the Asian coast 41 .It should be noted that Ghiold and Hoffman 18 mentioned that the living Echinarachnius parma, currently inhabiting the northwestern Atlantic, likely migrated from the Northeast Pacific through the Arctic during the Late Pliocene climatic amelioration (Fig. 3).Extant Sinaechinocyamus, on the other hand, subsequently expanded to tropical area and is nowadays abundant in Taiwan.
The inferred historical biogeography of the Luminacea suggests that the common ancestor of this group originated from the EIWP in the early Cretaceous (Fig. 3).The Cassiduloida and Echinolampadoida, appear to have dominated the world's echinoid diversity, with 60% of all Eocene echinoids belonging to these groups 42 .Since then, the number of species decreased dramatically 7,14 .This implies that the evolution of cassiduloids and echinolampadoids following the Eocene was driven by a series o extinction events.The fossil record thus is fundamental to fully understand their evolutionary histories.Conversely, with dense taxon sampling of extant sand dollars sensu lato, their evolutionary histories are well illustrated.
Echinoids in general, particularly taxa inhabiting shallow subtidal to intertidal habitats, are sensitive to oscillating ocean temperatures 12,43 .The sand dollar sensu lato likely expanded from the origin tropical EIWP to adjacent areas between the late Cretaceous and Eocene, and before the Late Eocene-Oligocene Cooling.Shallow marine currents that regulate the heat flow and governs sea surface temperature is largely driven by wind directions in the troposphere and climate patterns on Earth 44,45 .The relatively higher sea temperatures and warm surface currents 46 during the late Cretaceous, Paleocene, and Eocene (Fig. 3) drove the dispersal of sand dollars from low-to high-latitudes.Two range expansions likely occurred during this period: the ancestral Laganidae expanded to the south (southern Australia and New Zealand), north (northwest Pacific), and west (tropical western Indian Ocean); and the MRCA of Taiwanasteroidea expanded to the northern Pacific realm (Fig. 3).Vicariance and secondary dispersals likely occurred at various regions as the climate shifted from hothouse to ice-age conditions after the Paleocene-Eocene Thermal Maximum (PETM) 47 (Fig. 3).For example, the Dendraster, Echinarachnius, and Scaphechinus seem to have adapted to cold-water and survived through modern northern Pacific and northwestern Atlantic; while the Mellitidae is currently the only group that inhabits tropical America 48,49 .The widely distributed genus, Clypeaster, likely originated in the western Tethys in the Eocene and has probably undergone several range expansions during the Miocene and spread worldwide.

Systematics
Infraorder Scutelliformes Haeckel, 1896.Diagnosis: One lantern support (auricle) on internal side of single interambulacral plate near peristome; interambulacral column ends with two plates at the apical disc; with lunules and/or noticeable marginal indentations; periproct is on oral side and near peristome.Superfamily Astriclypeoidea Lin in Lee et al. nov.. Diagnosis: Lunulate scutelliforms up to five lunules; periproct on oral surface; lunules with cross-linked lunule wall.
Remark: This superfamily is proposed here based on the strong molecular support (Fig. 2) and the predicted origination event occurs during the Eocene (Fig. 3).

Family Astriclypeidae Stefanini, 1912. Diagnosis: as for superfamily
Distribution: Indo-Pacific West (Fig. 1B) Remarks: Astriclypeus, Sculpsitechinus and Echinodiscus are included in this study Superfamily Taiwanasteroidea Lin in Lee et al. nov.. Diagnosis: Flattened scutelliforms without lunules.Remark: Wang 26 proposed the Superfamily Taiwanasteritida that includes the Family Taiwansteridae and the Family Fibulariidae.Both families are now considered not closely related and this superfamily as defined in Wang 26 is polyphyletic.Based on cladistic analyses, Wang 50 proposed the Superfamily Dendrasteracea, including Echinarachniidae, Dendrasteridae and Mellitidae, under the Suborder Scutellina, and Fibulariidae and Laganidae under the Suborder Laganina.Mooi 51 provided a detailed study on living species of Dendraster.During ontogeny, the position of periproct shifts from the margin toward the peristome in all three living Dendraster species.While the clade is well defined based on molecular evidence (Fig. 2), it is difficult to define morphologic synapomorphies for this superfamily because Dendraster has many autapomorphies.Although Sinaechinocyamus was small in size (< 2 cm), it originates in the Taiwan Strait based on the good fossil records in Taiwan since late Miocene 36 .The new clade Taiwanasteroidea proposed here include four extant families: Dendrasteridae, Taiwanasteridae, Scutellidae and Echinarachniidae.This is the only clypeasteroid clade that occurs on both West and East Pacific margins.Lee et al. nov.. Diagnosis: Lunulate scutelliforms with an anal lunule; with festooned lunule walls 25 and/or marginal indentations; periproct positioned between peristome and anal lunule.

Superfamily Mellitoidea Lin in
Remark: This superfamily is proposed here based on the strong molecular support (Fig. 2) and the predicted origination event occurs during the Miocene (Fig. 3).

Family Mellitidae Stefanini, 1912. Diagnosis: as for superfamily
Distribution: Their main habitats are along the western margins of North America and South America, and around the Caribbean Sea and Gulf of Mexico 48 .

Ethical approval.
No species of echinoderms collected in this study are listed in national laws as protected or endangered.Most of the collected organisms are not subjected to restriction by national or international laws and do not require special permission, except two specimens (see Table 2) collecting in 2017 in Kenting National Park, Taiwan (permit number 1071000145).2).The habitat of each collected specimen was shown in Table S2.Collected specimens were directly fixed in 95% ethanol prior to DNA extraction.Genomic DNA was extracted from spine muscle by scratching the oral side of the test using the QIAamp DNA Micro-Kit (Qiagen, Hilden).The specimens were kept as vouchers in National Taiwan University Museums (NTUM) (Table 2).
PCR amplification and sequencing.We performed polymerase chain reactions (PCR) to amplify two mitochondrial genes (cytochrome c oxidase I; cox1 and 16S ribosomal RNA; 16S), and two nuclear genes (28S ribosomal RNA; 28S and Histone H3; H3) from each collected specimen (Table 3).PCR reactions in a total volume of 25 µl, contained 20-80 ng of DNA template, 0.2 µM of each primer, and 12.5 µl 2 × EmeraldAmp GT PCR Master Mix (TaKaRa Bio.).PCR cycles included an initial denaturation step of 94 °C for 2 min, followed by 35 cycles of denaturation (94 °C, 30 s), annealing and elongation step (Table 3), and a final elongation step at 72 °C for 5 min.The purification and sequencing of PCR products were then carried out at Genomics Company (Taiwan).The obtained DNA sequence chromatograms were assembled using CodonCode Aligner v.6.0.2 (Codoncode Corporation, Dedham, MA, USA).The two overlapping cox1 sequences (cox1 part 1 and part 2) were assembled to form a longer cox1 sequence.

Phylogenetic analysis.
Following the current phylogenetic framework of the Echinoidea 7 , we included 25 species of Luminacea (Cassiduloida, Echinolampadoida, Clypeasteroida, and Scutelloida) and three non-Luminacea species from Spatangoida and Camarodonta in the present phylogenetic analyses.We included Colobocentrotus mertensii (Odontophora) as an outgroup to root the inferred phylogenetic trees.The newly obtained cox1, 16S, 28S, and H3 sequences together with their orthologous sequences from 17 echinoderm taxa retrieved from Genbank were included in the analyses (Table 2).We aligned the cox1 and H3 sequences by eye and 16S and 28S sequences by the MAFFT ver.7 multiple alignment program, applying default parameters 57 .The final trimmed alignments consisted of 1,269 bp, 578 bp, 1,145 bp and 309 bp DNA sequences for cox1, 16S, 28S, and H3, respectively.To gain the power in terms of phylogenetic resolution, we adopted an approach with simultaneous analysis for the multi-gene data.The separated gene datasets with a total of 28 common echinoderm taxa were herein combined for the subsequent analyses.Although some of the individual gene sequences of the included taxa in the concatenated dataset may come from different individuals or be missing at a particular gene locus, we do not expect this to cause a significant error in phylogenetic inference at the interfamilial level.
We performed a phylogenetic analysis using partitioned maximum likelihood (ML) method as implemented in RAxML v.8.0 58 with the GTR + G substitution model based on the compiled multi-gene dataset.The partitions were set by gene and by codon position for protein coding genes.The parameters including nucleotide substitutions and distribution of site rates in each partition were estimated independently.Thus, the estimated alpha value of the gamma distribution varies among partitions, with values of 0.070, 0.045, and 1.281 for three codon position of cox1, respectively, 0.212 for 16S, 0.063 for 28S, 0.888, 0.020, and 0.020 for H3.Nodal support was assessed by bootstrapping 59 with 1,000 pseudo-replicates.We also performed a Bayesian Inference (BI) analysis on the combined dataset with MrBayes v.3.2.6 60 on the CIPRES Science Gateway 61 .Substitution model for each partition was set following the best-fit partition scheme given by PartitionFinder 2 62 (Table S3).Four Markov chains were performed in each of two parallel runs for 20,000,000 generations, sampling every thousand generation and a heating temperature of 0.02.Convergence was evaluated using Tracer v.1.7 30 confirming all Effective Sample Size (ESS) values were over 200.The BI tree was summarized using a 50% majority rule tree.Support values from both ML and BI methods were used to evaluate the robustness of inferred phylogenetic relationships within Luminacea.

Divergence time estimation.
Divergence times of lineages were estimated using BEAST v.2.6.7 63 in a relaxed log-normal clock using the same dataset as for the phylogenetic analyses with eight well-documented echinoderm fossils (see below) for calibration.The Bayesian trees were estimated with a Yule model.We set a single partition for the dataset with GTR + G substitution model.Trees were linked across genes whereas clocks were set unlinked.We followed the previous study of sand dollar phylogeny 48 , by setting the distributions of priors for mean rates (ucld.mean)uniform, ranging from 0.001 to 1 and standard deviation (ucld.stdev)as exponential distribution with a mean equal to 0.3333.Four independent runs of 100 million MCMC generations each were performed and sampled every 1000th generations.Each run was initiated from a random starting time tree.We checked the parameter log files for convergence with Tracer v.1.7.The resulting trees files from the four independent runs were removed 10% of trees as burn-in in each run, and combined using LogCombiner v.2.6.7.The maximum clade credibility (MCC) tree with mean divergence times generated from BEAST was reconstructed using TreeAnnotator v.2.6.4.
Our phylogenetic tree was time-calibrated using echinoid fossils that provide hard minimum and soft maximum bounds through exponential distributions in which the 95% credibility interval was equal to the maximum age of the strata where the first known fossil excavated, or the maximum bound of the most recent common ancestor (MRCA) that has been estimated in previous studies.The fossils and references used in the analysis were shown in Table S1.

Ancestral area reconstruction.
The origin and patterns of geographical diversification of Luminacea taxa were assessed by ancestral area reconstruction using the Dispersal-Extinction-Cladogenesis (DEC) model implemented in RASP v. 4.2 31 with the time-calibrated consensus tree reconstructed using BEAST (see above).Terminal taxa were regarded as the representatives for their genera (Fig. 3).We defined eleven biogeographical units for the included genera according to previous studies 15,18 with minor modifications based on oceanic basin, continental shelf, surface currents, and regional endemicity 18

Biogeography network visualization. For further geographic visualization of biogeography we built
StrainHub networks from the phylogenetic data and metadata using the strainhub R package 39 .The network of connections between geographic provinces is based on the tree used in RASP analysis and the source/hub ratio (SHR).The size of the circles represents the SHR of the nodes in the network (Fig. 4B) with larger, circles indicating the relative importance of a geographic province as a source of lineages.The thickness of the lines indicates the frequency of movement of lineages from one body of water to another.

Figure 2 .
Figure 2. Phylogenetic relationships of the Luminacea inferred using partitioned Maximum Likelihood analysis based on 3301 bp long concatenated multi-gene sequences.Asterisk represents polyphyletic Scutelloida.Nodal supports are shown as bootstrap (BS) values in percentage (above) and posterior probabilities (PP) (below).Values below 60% in BS and 0.8 in PP are not shown.Colored rectangles highlight the resolved main clades.Figure was made with Microsoft PowerPoint (https:// www.micro soft.com/, version 2016).

Figure 3 .
Figure 3. Most-likely ancestral range reconstruction of the Luminacea using the Dispersal-Extinction-Cladogenesis (DEC) model on the simplified Bayesian phylogenetic tree inferred by BEAST v.2.6.7 30 based on cox1, 16S, 28S, and H3 data.Outgroups were omitted from this analysis.Nodes represent the median divergence times.Values near nodes represent posterior probabilities (PP).PP values below 0.95 are not shown.Shapes of lunule from the corresponding scutelloid clades are shown to the right.Colored, lettered boxes on the nodes represent the most likely areas of origin (lower right; a. Tropical eastern Indian Ocean and western Pacific Ocean (EIWP); b.Southern Australia and New Zealand (SANZ); c.Northwestern Pacific (NWP); d.Tropical western Indian Ocean (WIO); e. Northeastern Pacific (NEP); f.Northwestern Atlantic (NWA); g.Tropical eastern Pacific (EP); h.Tropical western Atlantic and Caribbean Sea (WA); i. Northeastern Atlantic and Mediterranean (NEA); j.Tropical East Atlantic (EA); k.South Africa (SAFR)) reconstructed by RASP v.4.2 31 .Filled squares represent the constrained and assigned age prior nodes.Black arrows indicate inferred events of range expansions.Graph of Phanerozoic paleotemperatures was modified from Scotese et al. 32 .Figure was made with Microsoft PowerPoint (https:// www.micro soft.com/, version 2016).Image credit: Jih-Pai Lin.
Equipment and settings.Images of sand dollars in Figs. 3 and S3 were taken by J.-P.L. using a Nikon D750 digital SLR camera fitted with a macro lens (AF-S VR Micro-Nikkor 105 mm f/2.8GIF-ED) and were edited in Adobe Photoshop CS6.Brightness and contrast were adjusted with the software ImageJ.

Table 2 .
Species and sequences used in the current phylogenetic analysis.
a, b, c : sequences being combined into single taxon in phylogenetic analyses.Voucher specimens are deposited at the marine invertebrate (inv-) collection of the National Taiwan University Museums (NTUM).New sequences generated in this study are indicated in bold.AB900169 Vol.:(0123456789) Scientific Reports | (2023) 13:10199 | https://doi.org/10.1038/s41598-023-36848-0www.nature.com/scientificreports/Sample collection and DNA extraction.A total of 12 specimens of echinoids were newly collected and exanimated in this study (Table

Table 3 .
Primers, annealing temperatures and number of cycles used for PCR and sequencing.