Evolution and Viability of Asian Horseshoe Crabs Appear Tightly Linked to Geo‐Climatic Dynamics in the Sunda Shelf

Abstract

ABSTRACT

Horseshoe crabs are unique living fossils that have remained almost unaltered through 400 million years of global change. They face rapid worldwide declines under increasing anthropogenic pressure. Using comprehensive geographic and genomic sampling combined with approaches that integrate DNA with environmental and climatic datasets, we assessed the population genetic structure, demographic histories, and vulnerability to future climate change in three out of four extant horseshoe crab species, all centered in Asia. Our study highlights that the Sunda Shelf, a complex and dynamic shallow‐marine landscape, has been the sole repository of most genetic diversity among all three Asian species, and therefore crucial to the long‐term survival of horseshoe crabs. Our study not only provides the first genomic baseline data for the evaluation of Asian horseshoe crabs’ conservation status but also identifies core habitats that potentially act as refugia and corridors for Asian horseshoe crab populations with impending anthropogenic global warming.

Introduction

Horseshoe crabs are one of the planet's foremost “living fossils”. They are known for their high morphological conservatism and a consistently slow evolutionary rate over hundreds of millions of years (Bicknell and Pates 2020). Although they are a fundamental building block of coastal marine ecosystems where they occur (Botton 2009), only four geographically restricted species remain today (Sekiguchi and Shuster 2009): the Atlantic horseshoe crab (Limulus polyphemus) along the Atlantic coast of the United States and the Gulf of Mexico, and three Asian horseshoe crabs, the mangrove horseshoe crab (Carcinoscorpius rotundicauda), coastal horseshoe crab (Tachypleus gigas), and tri-spine horseshoe crab (Tachypleus tridentatus), in coastal East, Southeast, and South Asia (Figure 1). Although anecdotal reports suggest declines of horseshoe crabs worldwide (John et al. 2018; Wang et al. 2020), important baseline population data remain scant and are mostly only available in economically advanced regions. Hence, the Atlantic horseshoe crab is more intensively studied compared with its Asian cousins (Luo et al. 2020), two of which (C. rotundicauda and T. gigas) are currently classified as data deficient in the IUCN Red List, whereas T. tridentatus is considered endangered.

[IMAGE OMITTED. SEE PDF]

The biology of all four extant horseshoe crab species has been tightly adapted to ephemeral coastline habitat and shallow-marine conditions since the late Paleozoic (Blażejowski et al. 2017; Bicknell and Pates 2020). The three Asian species’ ranges broadly overlap across the Southeast Asian Sunda Shelf, which is known as one of world's biodiversity hotspots for terrestrial (Myers et al. 2000), freshwater (He et al. 2018), and coastal marine (Hoeksema 2007; Schumm et al. 2019) flora and fauna. Owing to accelerated diversification rates fueled by habitat dynamics across glacial cycles, the Sunda Shelf has become a cradle of terrestrial and freshwater endemism following its first submergence at ∼400 ka (Cros et al. 2020; Husson et al. 2020; Salles et al. 2021; Sholihah et al. 2021; Garg et al. 2022). The region also exhibits high marine biodiversity but relatively low endemism (Costello et al. 2017), suggesting niche filling by dispersion (Ludt and Rocha 2015; Pinheiro et al. 2017), which may have turned the Sunda Shelf into a sanctuary for coastal marine biodiversity during interglacial periods.

To understand population structure, evolutionary histories, and climate change-driven vulnerability, we comprehensively sampled the three Asian horseshoe crab species across their natural range. Our results provide baseline data for conservation action geared toward the continued existence of a species group that has survived global change almost unaltered since before the age of the dinosaurs.

Materials and Methods

Covering the entire extant geographical distribution of the three Asian horseshoe crab species (Figure 1), we collected tissue samples from 251 individuals (115 C. rotundicauda, 74 T. gigas, and 62 T. tridentatus) across 52 sites in 11 countries (Table S1). Single nucleotide polymorphisms (SNPs), produced with double digest restriction-site-associated DNA sequencing, were identified using a modified pipeline from our previous horseshoe crab population genomic study (Tang et al. 2021). Based on various SNP datasets (Table S2), we first estimated present-day genetic differentiation among individuals and delimited populations within each species (see Appendix S1 for details). To investigate the evolutionary histories of these populations, we modeled habitat suitability based on 11 paleo-climatic variables (Table S3) available from the online database CHELSA (Karger et al. 2023) using MaxEnt v3.4 (Phillips et al. 2017). We also modeled the demographic dynamics over time using genomic coalescent simulations, especially as implemented in Stairway Plot v.2 (Liu and Fu 2020) and FastSimCoal v.2.7 (Excoffier et al. 2021). For FastSimCoal analyses, we tested divergences among population pairs with different scenarios of gene flow (Figure S1). For both programs, we generated folded site frequency spectra including invariant sites as input, and set the mutation rate to 1.918 × 10⁻⁸ per site per generation and the generation length to 14 years (Shingate et al. 2020) for all three horseshoe crab species (see Appendix S2 for details). Finally, we adopted a novel analytical framework (Chen et al. 2022) to assess the climate change-driven vulnerability of Asian horseshoe crab populations by the years ∼2050 and ∼2100 under two greenhouse gas concentration projections, namely the representative concentration pathways (RCPs) RCP2.6 and RCP8.5. The geo-climatic data (Table S3) were acquired from Bio-ORACLE v2.0 (Assis et al. 2018). In addition, resistance to dispersal under future climate scenarios was calculated based on the correlation between genetic distances and geo-climatic variables using the R package resGF (Vanhove and Launey 2023) to evaluate the feasibility for horseshoe crabs to mitigate the adverse effects of climate change through migration (see Appendix S3 for details).

Results

Population Structure and Differentiation

All three Asian horseshoe crab species exhibited substantial genetic differentiation across their geographical distributions according to principal component (PC) analyses (Figure 1 and Figure S2). These results were consistent across eight SNP datasets subjected to different quality filters (Figure S3). In C. rotundicauda, individuals from North Borneo are highly differentiated from all other populations along PC1, whereas the remaining individuals display a clinal east-west divergence along PC2. In T. gigas, individuals from the Indian Ocean and the Pacific Ocean are divided along PC1, with a higher dispersion of Indian Ocean individuals across both PC1 and PC2. In T. tridentatus, there is a clinal north-south divergence along PC1. Population assignments reveal that C. rotundicauda emerges with more deeply structured populations as compared to T. gigas and T. tridentatus, as indicated by a higher number of populations detected (Figure 1, Figure S4, Table S2).

Based on the population assignments of the three species and the boundaries of marine ecoregions (Spalding et al. 2007), we divided the natural range of the three Asian horseshoe crabs into seven geographical regions (Figure 1): (A) Northern Bay of Bengal + Andaman Sea Coral Coast, (B) Malacca Strait, (C) Gulf of Thailand + Southern Vietnam + Sunda Shelf, (D) Gulf of Tonkin + Southern China, (E) Java Sea + Makassar Strait, (F) Palawan/North Borneo, and (G) East China Sea + Central Kuroshio Current.

A suite of analyses converges on a consensus population assignment (Figure 1) in which individuals of C. rotundicauda are grouped into six populations (A–F); T. gigas emerges as a single population but likely diverging into two (Northern Bay of Bengal + Andaman Sea Coral Coast vs. elsewhere); T. tridentatus is grouped into two populations delimited by the East China Sea (D+E+F vs. G). Based on the consensus population assignments, we estimated parameters to characterize diversity within each population and differentiation among populations. Within C. rotundicauda, population D, which marks the northernmost boundary of the species distribution (Figure 1), exhibits a much lower nucleotide diversity compared to the other populations (Table S4). In T. gigas, individuals in the Bay of Bengal and the Andaman Sea are characterized by a higher nucleotide diversity than individuals elsewhere. In T. tridentatus, the northern population G has a lower nucleotide diversity than the southern one (D + E + F). Analysis of molecular variance (AMOVA; Table S5) suggests that divergence among populations is higher in C. rotundicauda (37.2%–53.9% across eight SNP datasets) compared to T. gigas (comparing Northern Bay of Bengal + Andaman Sea Coral Coast with elsewhere, 11.9%–17.4%) and T. tridentatus (12.5%–23.1%). Measures of relative (F_st) and absolute (D_xy) population differentiation in C. rotundicauda (average F_st = 0.62 and average D_xy = 0.14 across datasets) are generally larger compared to T. gigas (0.17, 0.04) and T. tridentatus (0.19, 0.06), which is consistent with the results of AMOVA (Tables S6 and S7).

Species Distribution Modeling and Demographic Reconstruction

Species distribution modeling suggests markedly different evolutionary trajectories of the three horseshoe crabs during the last glacial maximum (LGM; Figure 2). Comparing the present day with the LGM, habitat availability has drastically increased for the tropical C. rotundicauda and T. gigas, whereas habitat availability has slightly declined for the temperate T. tridentatus. During the LGM, the range of C. rotundicauda was divided into multiple refugia isolated from one another along the east coast of the Sunda Shelf (Figure 2a); the range of T. gigas was limited to a single refugium along the north coast of the Andaman Sea (Figure 2b). In contrast, the range of T. tridentatus covered most of its present range, with ample availability of habitats in the East and South China Sea (Figure 2c). For C. rotundicauda, habitat suitability is mainly determined by bathymetry (48.4%), distance to land (34.1%), and mean salinity at sea surface (9.7%); for T. gigas, it is mainly determined by bathymetry (42.3%), distance to land (37.1%), and mean temperature at sea surface (9.8%), whereas for T. tridentatus, habitat suitability is mainly determined by distance to land (62.9%), bathymetry (25.1%), and the range of air temperature (4.8%). None of these top-contributing variables exhibit a high correlation to any of the other variables.

[IMAGE OMITTED. SEE PDF]

We reconstructed effective population size over time and identified the most likely demographic scenarios for the consensus populations of the three Asian horseshoe crab species (Figure 2) by comparing simulated data and observed data (Table S8). In C. rotundicauda, the first divergence of extant populations occurred between population F and all others at ∼1.6 Ma (Figure 2a and Table S9). Subsequent divergences were staggered between 330 and 225 ka, following the peaks of two sequential glaciations (Figure 2a, Figure S5, Table S9). After these divergences, effective population sizes generally remained stable across populations with relatively independent fluctuations, which may have been subject to local environmental dynamics during glacial cycles (Figure 2a). Populations in the historic center of the distribution (A, B, and C) were subject to milder fluctuations than those at the margins (D, E, and F). T. gigas fluctuated and grew 10 times in effective population size over the last 1 million years without undergoing population differentiation (Figure 2b). The two populations of T. tridentatus diverged at ∼550 ka. Population G declined considerably after its divergence, while population D + E + F has maintained its effective population size. Gene flow (mainly from G to D + E + F) has continually occurred between the two populations after the first divergence (Figure 2c and Table S9).

Climate Change-driven Vulnerability

Genotypes of the three Asian horseshoe crab species have different associations with various environmental variables: the smaller C. rotundicauda is more sensitive to temperature, whereas the larger T. gigas and T. tridentatus are more sensitive to salinity and concentration of chlorophyll (Figures S6 and S7). Based on their genotype–environment associations, genomic offset values of the three horseshoe crab species, which reflect turnover of their allele frequencies in response to future climate change, are expected to increase by the end of this century especially under high emission scenarios (Figures S8–S10), highlighting horseshoe crabs’ vulnerability to global warming. In contrast, ecological niche modeling, incorporating additional topographical variables relevant to the ecological niche delimitation of horseshoe crabs, indicated that intensifying climate change would result in overall increases in niche suitability, with slight losses near the equator offset by gains at greater latitudes (Figures S11–S13). When combining genomic offset values with changes of niche suitability, the resulting genome-niche index turned out to be coherent with genomic offset patterns across all three horseshoe crab species (Figures S14–S16). Overall, all three species of Asian horseshoe crabs will face substantial pressure to adapt to the climate later during this century unless the current level of emissions can be reversed (under RCP2.6). In 2050 under RCP8.5, some populations are predicted to be more vulnerable to climate change than others: the northerly C. rotundicauda population in the Gulf of Tonkin and South China (D) and most populations of T. gigas (except the ones in the Bay of Bengal and Southern Vietnam) will be among the most highly impacted (Figure 3).

[IMAGE OMITTED. SEE PDF]

To understand how natural dispersal may mitigate horseshoe crabs’ vulnerability to climate change, we mapped the distribution of dispersal resistance. The patterns of resistance to dispersal across the three Asian horseshoe crab species are generally consistent between the two emissions scenarios (Figures S17–S19). Resistance to dispersal in C. rotundicauda is overall the highest among the three Asian horseshoe crab species, even considering that relatively low levels of resistance to dispersal are predicted along coastlines, except for the east coast of Borneo (Figure 3). In C. rotundicauda, 67.36% of SNPs exhibited a correlation with at least one of the 19 environmental variables that predicted 62.46% of the resistance pattern. High resistance to dispersal was identified to be associated with a large range of temperatures at both the minimum depth (10.64%, Figure S20) and sea surface (10.52%, Figure S21) as well as a small range of chlorophyll concentrations (9.45%, Figure S22). Resistance to dispersal in T. gigas was overall the lowest, with a relatively homogeneous distribution across its range. In T. gigas, 12.56% of SNPs exhibited a correlation with at least one of the 19 environmental variables that predict 69.63% of the resistance pattern. High resistance to dispersal was identified to be associated with a large range of current velocity at sea surface (10.32%, Figure S23), steep south-facing slopes (9.66%, Figure S24), and low mean temperatures at sea surface (9.11%, Figure S25). In T. tridentatus, resistance to dispersal was predicted to be low across coastal East Asia but high in Palawan and the north and east coast of Borneo. A total of 34.74% of SNPs in this species exhibited a correlation with at least one of the 19 environmental variables that predicted 75.17% of the resistance pattern. High resistance to dispersal was identified to be associated with low mean chlorophyll concentrations (7.59%, Figure S26), deep sea trenches (7.11%, Figure S27), and a small range of current velocity at the minimum depth (7.04%, Figure S28).

Discussion

Our study sheds light on the population structure of all three Asian horseshoe crab species and identifies genetically distinct populations that are suitable to serve as conservation units for future efforts to optimize the preservation of each species’ evolutionary potential for adaptation to local environmental dynamics. Our analyses consistently uncover Sundaic populations at the center of evolutionary dynamics, while peripheral populations—oftentimes those at cooler latitudes—appear to play a marginal role in the long-term persistence of each species (Figure 1). In Sundaic populations, the period of highest diversification (for the less dispersive C. rotundicauda) or expansion (for the dispersive T. gigas and T. tridentatus) commenced after ∼400 ka when the former subcontinent of Sundaland began to be periodically inundated (Figure 2) (Sarr et al. 2019; Salles et al. 2021). This pronounced geographical change, which marked the birth of the Indonesian Archipelago as we know it today, has resulted in a flourishing of horseshoe crabs and potentially other coastal marine species (Ludt and Rocha 2015), and is paralleled by findings of rampant diversification of terrestrial and freshwater taxa (Salles et al. 2021; Sholihah et al. 2021), highlighting the unique contributions of regional environmental forcing (e.g. tectonic, sea-level fluctuations) to global biodiversity. The complexity of the Sunda Shelf has allowed for the persistence of refugia that have facilitated the survival of horseshoe crabs through the vagaries of past climate fluctuations and will be indispensable for their survival during impending episodes of anthropogenic climate change (Figure 3). Unlike most other marine regions across the world, the shallow and stable shelf area in the Sundaic region has facilitated benthic long-distance dispersal of horseshoe crabs across long evolutionary periods of time, a potent mechanism to counteract the negative effects of rapid local changes in habitat suitability (Figure 3).

All three species have incurred gains in genetic diversity as a result of the Sunda Shelf's inundation (Figure 2). On the other hand, recent observed declines in horseshoe crab populations due to human activities (John et al. 2018) are probably not reflected in our dataset, given that horseshoe crabs have long generation times (Petit-Marty et al. 2022). Continuous monitoring of both census and effective population sizes remains essential for their management (John, Tang, and Eackles 2022). According to our genome-niche index estimates (Figure 3, Figures S14–S16), all three Asian horseshoe crabs are vulnerable to future climate change, especially in light of their long generation times (∼14 years), which severely slow their rate of genomic turnover and constitute a nearly insurmountable obstacle to the kind of rapid adaptation that will be required under future climate change scenarios. Possible range shifts and natural migration may mitigate the environmental stress on horseshoe crab populations in the future. In C. rotundicauda, populations along the coast of Southern China and the Gulf of Tonkin are subject to some of the highest evolutionary pressures to adapt to future climate change. Their likely path to survival may be to shift southwards (Figure 3). Appropriate protection and restoration of mangrove habitats will be imperative to maintain mangrove horseshoe crabs’ migratory potential to mitigate such evolutionary pressure. In T. gigas, which is most vulnerable to future climate change among the three Asian horseshoe crab species, the Sunda Shelf may provide sufficient refugial habitat around its periphery (e.g., around the Bay of Bengal, Malacca Strait, and Southern Vietnam) to maintain connectivity among populations (Figure 3). In T. tridentatus, impending climate change is predicted to produce only a moderate influence on species viability (Figure 3) as compared to human activities such as harvesting and coastline modifications (Laurie et al. 2019). Therefore, its conservation should focus on regulations to promote sustainable fishery practices and the restoration of coastal habitats, especially in areas with a relatively long history of intensive development (e.g., Japan, Taiwan, and China). Understanding the distribution of refugial habitats and natural migration routes will be helpful for coastal conservation planning specifically targeting hotspots for the viability of horseshoe crabs. Our study provides an important impetus and the necessary baseline data for the preservation of key habitats for horseshoe crabs’ future survival. As an important caveat, however, our work is only based on environmental factors and does not take into account future human activities that may directly alter habitats, such as coastal development. The survival of horseshoe crabs will therefore critically depend on interventions (e.g., legislation in fishery) based on local contexts.

Acknowledgments

This work was supported by the Marine Science Research and Development Program, National Research Foundation, Prime Minister's Office, Singapore (MSRDP-P19 and MSRDP-P41), and Ministry of Education, Singapore under a Tier 2 grant (MOE-T2EP30120-0018). Sample collection in Japan, the Philippines, and parts of Indonesia was supported by the Global COE Program (Center of Excellence for Asian Conservation Ecology as a Basis of Human-Nature Mutualism), MEXT, Japan, and a JSPS KAKENHI Grant (23770090). Sample collection in Vietnam was supported by the Vietnam Academy of Science and Technology (UQĐTCB.01/23-25). Sample collection in China was supported by the Hainan Natural Sciences Foundation Project (322MS154) and the Natural Science Foundation of Guangdong Province of China (2024A1515012696). Sample collection in India was supported by the Arcadia, Conservation Leadership Programme (CLP), UK. We thank Prashant Shingate and Vydianathan Ravi for providing suggestions on population genomic analyses. We thank Naila Khuril Aini, the Japanese Society for the Preservation of the Horseshoe Crab (Nihon Kabutogani wo Mamoru Kai), Hiroko Koike, Mitsuhiro Aizu, Ryo Nogami, Ame M. Garong, Irzal Azhar, Ali Mashar, Irfan Yulianto, Marivene Manuel-Santo, Dinh Dieu Thuy, Mijanur Rahman, and Omar Hasan for assisting in horseshoe crab collections.

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

All data are deposited under Bioproject PRJNA1127623 at NCBI. Codes used for all analyses are available at https://github.com/qt37t247/Horseshoe-crab-ddRAD. DNA and tissue samples used in this study are deposited at the Global Horseshoe Crab Biorepository (https://github.com/qt37t247/horseshoecrab_bio-repository) in the Lee Kong Chian Natural History Museum, Singapore.

References

Assis, J., L. Tyberghein, S. Bosch, H. Verbruggen, E. A. Serrão, and O. De Clerck. 2018. “Bio‐ORACLE v2.0: Extending Marine Data Layers for Bioclimatic Modelling.” Global Ecology and Biogeography 27: 277–284.

Bicknell, R., and S. Pates. 2020. “Pictorial Atlas of Fossil and Extant Horseshoe Crabs, With Focus on Xiphosurida.” Frontiers in Earth Science 8: 98.

Blażejowski, B., G. Niedźwiedzki, K. Boukhalfa, and M. Soussi. 2017. “Limulitella Tejraensis, a New Species of Limulid (Chelicerata, Xiphosura) From the Middle Triassic of Southern Tunisia (Saharan Platform).” Journal of Paleontology 91: 960–967.

Botton, M. L. 2009. “The Ecological Importance of Horseshoe Crabs in Estuarine and Coastal Communities: A Review and Speculative Summary.” In Biology and Conservation of Horseshoe Crabs, edited by J. Tanacredi, M. Botton, D. Smith, 45–63. Boston, MA: Springer.

Chen, Y., Z. Jiang, P. Fan, et al. 2022. “The Combination of Genomic Offset and Niche Modelling Provides Insights Into Climate Change‐Driven Vulnerability.” Nature Communications 13: 4821.

Costello, M. J., P. Tsai, P. S. Wong, A. K. L. Cheung, Z. Basher, and C. Chaudhary. 2017. “Marine Biogeographic Realms and Species Endemicity.” Nature Communications 8: 1057.

Cros, E., B. Chattopadhyay, K. M. Garg, et al. 2020. “Quaternary Land Bridges Have Not Been Universal Conduits of Gene Flow.” Molecular Ecology 29: 2692–2706.

Excoffier, L., N. Marchi, D. A. Marques, R. Matthey‐Doret, A. Gouy, and V. C. Sousa. 2021. “fastsimcoal2: Demographic Inference Under Complex Evolutionary Scenarios.” Bioinformatics 37: 4882–4885.

Garg, K. M., B. Chattopadhyay, E. Cros, et al. 2022. “Island Biogeography Revisited: Museomics Reveals Affinities of Shelf Island Birds Determined by Bathymetry and Paleo‐Rivers, Not by Distance to Mainland.” Molecular Biology and Evolution 39: [eLocator: msab340].

He, F., V. Bremerich, C. Zarfl, et al. 2018. “Freshwater Megafauna Diversity: Patterns, Status and Threats.” Diversity and Distributions 24: 1395–1404.

Hoeksema, B. W. 2007. “Delineation of the Indo‐Malayan Centre of Maximum Marine Biodiversity: The Coral Triangle.” In Biogeography, Time, and Place: Distributions, Barriers, and Islands, edited by W. Renema, 117–178. Dordrecht, The Netherlands: Springer.

Husson, L., F. C. Boucher, A.‐C. Sarr, P. Sepulchre, and S. Y. Cahyarini. 2020. “Evidence of Sundaland's Subsidence Requires Revisiting Its Biogeography.” Journal of Biogeography 47: 843–853.

John, B. A., B. R. Nelson, H. I. Sheikh, et al. 2018. “A Review on Fisheries and Conservation Status of Asian Horseshoe Crabs.” Biodiversity and Conservation 27: 3573–3598.

John, B. A., Q. Tang, and M. Eackles. 2022. “The Need to Establish Horseshoe Crab Global Biorepository.” In International Horseshoe Crab Conservation and Research Efforts: 2007‐2020, edited by J.T. Tanacredi, M.L. Botton, P.K.S. Shin, Y. Iwasaki, S.G. Cheung, K.Y. Kwan, and J.H. Mattei, 345–350. Cham: Springer International Publishing.

Karger, D. N., M. P. Nobis, S. Normand, C. H. Graham, and N. E. Zimmermann. 2023. “CHELSA‐TraCE21k–High‐Resolution (1 km) Downscaled Transient Temperature and Precipitation Data Since the Last Glacial Maximum.” Climate of the Past 19: 439–456.

Laurie, K., C. P. Chen, S. G. Cheung, et al. 2019. “Tachypleus tridentatus” (errata version published in 2019). The IUCN Red List of Threatened Species. T21309A149768986. [DOI: https://dx.doi.org/10.2305/IUCN.UK.2019-1.RLTS.T21309A149768986.en].

Wang, C. C., K. Y. Kwan, P. K. S. Shin, et al. 2020. “Future of Asian Horseshoe Crab Conservation Under Explicit Baseline Gaps: A Global Perspective.” Global Ecology and Conservation 24: [eLocator: e01373].