Present and Potential Future Distributions of Asian Horseshoe Crabs Determine Areas for Conservation

Vestbo, Stine; Obst, Matthias; Quevedo Fernandez, Francisco J.; Intanai, Itsara; Funch, Peter

doi:10.3389/fmars.2018.00164

ORIGINAL RESEARCH article

Front. Mar. Sci., 14 May 2018

Sec. Marine Conservation and Sustainability

Volume 5 - 2018 | https://doi.org/10.3389/fmars.2018.00164

This article is part of the Research Topic 5th International Marine Conservation Congress View all 23 articles

Present and Potential Future Distributions of Asian Horseshoe Crabs Determine Areas for Conservation

$\r\nStine Vestbo*$ Stine Vestbo¹^*

Matthias Obst²

Francisco J. Quevedo Fernandez³

Itsara Intanai⁴

Peter Funch¹

¹Genetics, Ecology, and Evolution, Department of Bioscience, Aarhus University, Aarhus, Denmark
²Department of Biological and Environmental Sciences, University of Gothenburg, Gothenburg, Sweden
³School of Computer Science and Informatics, Cardiff University, Cardiff, United Kingdom
⁴Department of Biology, Faculty of Science and Technology, Prince of Songkla University, Pattani, Thailand

Conservation of horseshoe crabs has recently received increasing attention as several populations are in decline. However, scarce information on their distributions in Southeast Asia is impairing conservation efforts. In this study, we sought to improve our understanding of the geographical range and distinct populations of the three Asian horseshoe crabs species in order to identify optimal conservation areas. We mapped the geographic range of Carcinoscorpius rotundicauda, Tachypleus gigas, and T. tridentatus using recent data from field work, literature, Global Biodiversity Information Facility (GBIF), and unpublished data from our scientific network. The data were correlated with 23 different environmental variables of potential ecological importance for horseshoe crabs using the openModeller webservices, including new tidal variables. Ecological niche models were generated using two algorithms, Maximum Entropy and support vector machine, for the three species under present conditions, and projected into a climate change scenario of 2050. The niches of the Asian horseshoe crabs were mostly determined by tidal regime, chlorophyll A concentrations, depth, distance to land, and sea surface temperature. According to our predictions, horseshoe crabs in Southeast Asia are not expected to experience any severe change in extent and distribution of suitable habitat in the future. In order to conserve Asian horseshoe crabs, we suggest establishing Marine Protected Areas at locations where distinct populations and several species occur, such as northern Vietnam, China, Borneo, and southern Japan.

Introduction

Coastal and intertidal areas are currently under threat globally due to a range of anthropogenic activities, including infrastructure development and coastal protection, as well as effects of climate change, such as rising sea levels pushing coastal areas closer to anthropogenic structures, resulting in “coastal squeeze” (Defeo et al., 2009). Coastal areas are characterized by high productivity and species richness (Ketchum, 1972; Ray, 1991), and humans depend strongly on a variety of species that are living and reproducing in the coastal zone, such as fishes (Allison et al., 2009; Bell et al., 2009), sea cucumbers (Purcell et al., 2013), and horseshoe crabs (Chatterji, 1994), thus the conservation of such species is of great importance.

During the past decade, horseshoe crab conservation has received increasing attention (Botton et al., 2015). Horseshoe crabs are marine chelicerates of the order Xiphosura, containing only four extant species. Three species occur in the coastal waters of India, Southeast Asia, China, and Japan [Carcinoscorpius rotundicauda (Latreille, 1802), Tachypleus gigas (Müller, 1785), and T. tridentatus (Leach, 1819)]. One species [Limulus polyphemus (Linnaeus, 1758)] occurs in the coastal waters of the eastern continental shelf of North America and in the Gulf of Mexico (Sekiguchi and Shuster, 2009). Populations of horseshoe crabs world-wide are currently in decline (e.g., Cartwright-Taylor et al., 2011; Kwan et al., 2016; Smith et al., 2017), and the IUCN Red List of Threatened Species now describes L. polyphemus as vulnerable (Smith et al., 2016). Although the Asian horseshoe crab species are thought to be in greater decline than L. polyphemus, these species are currently listed as data deficient (World Conservation Monitoring Centre, 1996a,b,c). While the geographic range of L. polyphemus is known (Sekiguchi, 1988; Anderson and Shuster, 2003; Sekiguchi and Shuster, 2009; Faurby et al., 2011), the precise distribution of the three Asian species remains less clear (World Conservation Monitoring Centre, 1996a,b,c), thus impairing conservation efforts.

The global decline in horseshoe crab populations is mainly due to anthropogenic activities. For L. polyphemus the decline has primarily been due to overharvesting and loss of breeding habitat (Smith et al., 2017), as this species breeds on beaches with slight slopes; a habitat commonly used for real estate construction and development (Nordstrom, 2004). Coastal armoring as a response to erosion is an additional factor reducing available breeding habitat (Jackson et al., 2015). Wild caught horseshoe crabs are used commercially in medicine, where substances from their hemolymph (Carcinoscorpius, Limulus, and Tachypleus Amoebocyte Lysate, CAL, LAL, and TAL) are used to test if drugs, blood products, and pharmaceutical devices are free from bacterial contamination (Levin and Bang, 1968; Rao and Bhagirathi, 1989; Levin et al., 2003). L. polyphemus is released after collection of LAL, but the two Tachypleus species are commonly used for chitin production after hemolymph collection, thus harvesting of these species results in 100% mortality (Gauvry, 2015). However, a synthetic, commercially available version of the CAL test has been developed (Ding et al., 1995; Lonza, 2016). The three Asian horseshoe crab species are fished for human consumption (Botton, 2001), and were previously imported to North America to be used as bait (Smith et al., 2016), although this import is now prohibited (IUCN, 2013). As for L. polyphemus, the expansion of urban infrastructure into coastal areas, as well as coastal armoring present additional anthropogenic habitat disturbances for the Asian horseshoe crabs, thus in some areas, such as Japan and Peninsular Malaysia, horseshoe crab breeding habitats are now almost completely lost due to coastal infrastructure development (Botton, 2001; Nelson et al., 2016). In response to the declining horseshoe crabs populations, Marine Protected Areas (MPAs) have been established, especially in the USA (Carl N. Shuster Jr. Horseshoe Crab Reserve, Delaware Bay) and Japan (Saikai National Park in Nagasaki). However, in most of Southeast Asia MPAs protecting critical horseshoe crab habitats are still scarce.

Climate change might have additional detrimental effects on the distribution and population sizes of horseshoe crabs. Rising sea levels could reduce suitable horseshoe crab breeding habitats, since the retreating shorelines would bring the horseshoe crabs closer to the human infrastructure near the coasts. This is the case for Delaware Bay, New Jersey, USA, where sea levels have risen at a high rate during the twentieth century and hence decreased horseshoe crab breeding habitat (Loveland and Botton, 2015). Climate change would also affect ecological factors that influence the developmental success of horseshoe crab eggs and larvae, such as salinity and temperature (Jegla and Costlow, 1982; Laughlin, 1983; Ehlinger and Tankersley, 2004; Zaleha et al., 2011), oxygen levels (Palumbi and Johnson, 1982; Funch et al., 2016), beach geochemistry, wave energy, and erosion (Botton et al., 1988; Penn and Brockmann, 1994; Jackson et al., 2008). One of the species, C. rotundicauda, is found in habitats characterized by mangroves, which are also vulnerable to the effects of climate change, primarily rising sea levels (Gilman et al., 2008).

Ecological niche modeling (ENM, also known as species distribution modeling) allows for predictions of suitable habitats for a given species by identifying the environmental factors driving the species' distribution, such as temperature, salinity, or depth (Peterson et al., 2011). Furthermore, ENM can be used to project the distribution of suitable habitats of a species into future scenarios (e.g., future climate change scenarios), allowing forecasts of possible species range shifts. For example, the method has been used to predict future distributions and range shifts for Mediterranean fishes (Albouy et al., 2013) as well as changes in suitable habitat for commercially important fish and invertebrates along the Atlantic coast of North America (Kleisner et al., 2017). ENM can be regarded as an important tool for conserving threatened species, and has proven useful for the assessment and planning of protected areas (Kremen et al., 2008; Stirling et al., 2016). Thus, the scattered information on Asian horseshoe crab distributions and their unknown conservation status would benefit from gathering recent occurrences as well as a modeling of suitable habitats in the region. The objectives of this study were to (1) identify important environmental factors that influence the distributions of the three Asian species of horseshoe crabs, (2) update the information on current species distributions, and (3) identify possible high priority regions for their conservation, where climate change has minimal effect and where species ranges overlap.

Materials and Methods

Species Data and Study Area

The geographic study area included all coastal waters between 10.0° S to 35.0° N and 77.0° W to 135.0° E, thus including the continental shelves of the Bay of Bengal, India and Indonesia and between Indonesia and Japan. This area, the geographical mask, was chosen based on the existing knowledge of horseshoe crab distribution limits in Asia (Sekiguchi, 1988), and used for geographical filtering of occurrence records, collecting background points, and building and projecting the models. We assembled occurrence records (presence-only) for all three species in this region from our own observations, collaborators, scientific networks as well as through publishing a scratchpad site at http://horseshoecrabs.myspecies.info/ (Table 1, Supplementary Table 1). For many species, numerous distribution records exist in the literature, and we manually geo-referenced additional occurrence data from these sources (Supplementary Table 2). Subsequently, we added data from the Global Biodiversity Information Facility (http://gbif.org) in September, 2015 (Supplementary Table 3). A taxonomic data refinement workflow (Mathew et al., 2014) was used to check the synonyms, download, visualize, filter, and integrate occurrence records for all species. We inspected all records and excluded those of dubious identity (e.g., on land, in deep water, or without clear reference of origin), as well as those older than 1995 and those outside the geographical mask. All new data (i.e., excl. GBIF records) were submitted to the Ocean Biogeographic Information System (http://www.iobis.org/) and are available under the link https://doi.org/10.14284/293.

TABLE 1

Table 1. Occurrence data used in the analysis. Environmentally unique points (EUPs) were calculated from the compiled observations (records) using the respective geographical mask with a resolution of 300 arc seconds and 23 environmental variables.

Environmental Data

Based on available information on general ecology of horseshoe crabs (Bonaventura et al., 1982; Sekiguchi, 1988; Shuster et al., 2003), we included 23 environmental variables in the analysis, relating to sea surface temperature, chlorophyll A concentration, salinity, oxygen and nutrient concentrations, distance to land, water depth, pH, and tidal regime (Table 2). Chlorophyll A variables were included since chlorophyll levels acts as a proxy for phytoplankton primary production, which in turn reflects food availability for filter-feeding bivalves, which is one of the most important food items for horseshoe crabs (Chatterji, 1994). We used 13 marine layers from Bio-Oracle (http://www.bio-oracle.ugent.be/) with a resolution of five arc-minutes (Table 2; Tyberghein et al., 2012). These data layers are generated from monthly satellite data (Aqua-MODIS and SeaWiFS; https://oceancolor.gsfc.nasa.gov) as well as in situ measured oceanographic data from the World Ocean Database, 2009 (Boyer et al., 2009), and overlap with the species occurrence data in their temporal origin. In addition, we used five marine layers for both present and 2050 climate scenarios from AquaMaps (http://www.aquamaps.org/download/main.php) with a resolution of 30 arc-minutes (Table 2; Kaschner et al., 2008). Present-day (i.e., year 2010) datasets from AquaMaps were built from long-term averages of temporally varying environmental variables (Ready et al., 2010), while future data sets for 2050 were derived from the IPSL-CM4 A2 climate change scenario (IPCC, 2007).

TABLE 2

Table 2. Environmental variables investigated for ecological niche modeling of the three Asian horseshoe crabs.

We additionally generated five present-day (i.e., year 2010) global tidal variables with the Finite Element Solution oceanographic model (FES2012), provided by Noveltis, Legos and CLS Space Oceanography Division and distributed by AVISO+ (http://www.aviso.altimetry.fr/) (Table 2, more detailed information on how the global tidal variables were created can be found in Supplementary Data Sheet 1). FES2012 is a fully revised version of the original global hydrodynamic tide solution model (Lyard et al., 2006) with improvements regarding longer altimeter time series, more accurate ocean bathymetry, major non-linear tides, and more accurate tidal currents. It includes overall 32 tidal constituents distributed on 1/16° grids (amplitude and phase), corresponding to 3.75 arc-minutes (Carrère et al., 2012). The tidal variable layers can be accessed through a webservice (https://www.biodiversitycatalogue.org/services/37), and is hence available for future research on distribution modeling of marine intertidal organisms. In addition, we deposited the raster files together with the algorithm for calling the FES program and creating the ASCII grids at the Swedish Environmental and Climate Data Repository (www.ecds.se) under the identifier 87196c86-8e36-4908-8f16-424e61aa313c.

Analysis of Environmental Variables

Environmental values were obtained using version 2 of the BioClim workflow (http://purl.ox.ac.uk/workflow/myexp-3725.2) available at the BioVeL portal, https://portal.biovel.eu (Hardisty et al., 2016). The workflow was used to retrieve environmentally unique points (EUPs) from the species occurrence files for a set of 23 environmental data sets (Table 2). Since all points within a raster cell share the same values for the environmental variables, EUPs are obtained by filtrating the occurrence points, so that each raster cell contains no more than a single point (Nix, 1986). The workflow returned environmental values from the specified layers that matched with the species occurrence records, and it was executed in batch mode (called data sweep function) to repeat all calculations with the same parameter settings for all species.

We analyzed differences in environmental tolerances between all three horseshoe crab species by applying a Welch's test for unequal variances to the retrieved environmental values (Supplementary Table 4). This test was chosen because the data did not conform to parametric test assumptions. In addition, pairwise comparisons of environmental tolerances between species were performed using the Games-Howell Post-Hoc test (Supplementary Table 5). Both statistical tests were performed using SPSS version 22 (IBM Corporation, 2013). We subsequently performed a correlation analysis (data not shown) as well as a principal component analysis (PCA) on the environmental values using the R statistical environment 3.0.2 (R Core Team, 2013). The analyses were performed in order to identify ecologically important variables that explain the variation in the data set, and which can be used as predictor variables (Metzger et al., 2005; Porfirio et al., 2014). The PCA was estimated as a 23-dimensional hypervolume and used to compare the environmental space occupied by the three species. We used all non-correlated variables that were represented (with >5%) in the first two components of the PCA to build the niche models, except for Mean SST, which was included in the future projection (Supplementary Tables 6–8).

Ecological Niche Modeling

We used version 20 of the ecological niche modeling (ENM) workflow (http://purl.ox.ac.uk/workflow/myexp-3355.20) in batch mode (called data sweep) to describe and compare the geographical space occupied by the three species and to estimate the distribution of suitable habitat with favorable biotic, environmental, and geographical conditions. Demographic or dispersal properties that may also be used in species distribution modeling were not considered (Peterson et al., 2011; Reiss et al., 2014). For more information on the ENM workflows see Holl et al. (2013), Leidenberger et al. (2015a,b), and De Giovanni et al. (2015). We executed parallel analyses with two ENM algorithms by means of the openModeller webservice suite (de Souza Muñoz et al., 2011). These are i) Maximum Entropy v. 1.0 (Phillips et al., 2006; Phillips and Dudík, 2008) and ii) support vector machine v. 0.5 (Schölkopf et al., 2001). These two methods were chosen because Maximum Entropy (Maxent) is a presence-background algorithm that has been shown to perform with high predictive accuracy even if there are moderate errors in the geo-referencing (Elith et al., 2006; Graham et al., 2008), and support vector machine (SVM) is a machine-learning algorithm that is not analyzing characteristics of statistical distributions and hence does not require independent observation data, and thereby overcomes potential problems with spatial autocorrelation in the species occurrence data (Drake et al., 2006).

Models were created using each species' maximum distribution range within the mask and a set of predictor variables identified in the PCA. Variables selected for present day projections (high resolution variables) and future trend analysis (low resolution variables) are listed in Table 2. Models were created based on EUPs (Table 1) with the following specifications. Maxent models were set to run with 10.000 background points (including input points) drawn from the mask. Feature selection was automated, allowing the algorithm to combine feature types when fitting a model, and perform 500 iterations. Tolerance for detecting model convergence was set to 0.00001, while sample threshold was set to 80 (product), 10 (quadratic), and 15 (hinge). SVM models were set to execute the C-SVC algorithm with radial basis kernels, gamma values 1/k (where k is the number of layers), and a cost value of 1. All models were set to produce a probabilistic, instead of binary output. Predictive models may encounter some problems with extrapolations (Zurell et al., 2012), i.e., when the environmental hyperspace of the scenario is not represented by true observations. These were minimized by using only few essential variables in the future predictions, thereby avoiding complexity and interaction between too many factors, as well as through using a consensus approach (Table 2; Pearson et al., 2006).

For each species, we ran both algorithms across two sets of environmental layers (specified above and in Table 2); one with high resolution (using Bio-Oracle and FES2012 layers) and one with low resolution (using AquaMaps layers). The high-resolution models were used for present day projections of suitable habitat, while the low-resolution models were used to deduct climate driven changes in suitable habitat. Models were tested using 10-fold cross-validation based on the area under the curve (AUC) value and omission error rate (false negative rate), and subsequently projected using present and (in case of AquaMaps layers) 2050 climate scenarios (Table 3). The results of the ENMs were visualized as maps showing the distribution of suitable habitat of a species in a particular geographic region at a particular time (de Souza Muñoz et al., 2011). The sweep function allowed automated batch processing of all species for a given algorithm and parameter set. Overall, we executed 12 niche models, i.e., two algorithms, three species, and two environmental datasets (Table 3).

TABLE 3

Table 3. Results of model tests for all models of current distribution (based on Bio-Oracle and tidal variables) and future trend analysis (based on AquaMaps variables).

GIS Analysis

Raster maps created by the niche modeling algorithms were processed using the qGIS software package v. 2.6 Brighton (Quantum GIS Development Team, 2014). Maps of individual models are available in the Supplementary Material (Supplementary Images 1–6). For the consensus analysis presented in the maps, we transformed each individual model into a binary output using the lowest presence threshold (LPT) in order to include any signal of suitable habitat. Thereafter the maps of both models were merged into a single consensus model highlighting only areas where both algorithms agreed on the prediction of suitable habitat. Likewise, heat maps showing the predicted changes between 2050 and present-day distributions were created as consensus maps identifying only areas where both algorithms agreed on predicted loss or gain of suitable habitat. Finally, we compared high-resolution models (used for present day projections) with low-resolution models (for future change projections). We found that both model settings predict suitable habitat for similar biogeographic regions with few deviations mentioned in the discussion below.

Results

Environmental Factors Defining Horseshoe Crab Distributions

The PCA showed that C. rotundicauda and T. gigas largely share the same ecological niche, and that the niche of T. tridentatus to some extent overlaps with the niches of the other two species (Figure 1). Further, the PCA showed that the distributions of the three species are mainly driven by chlorophyll concentration, tidal regime, temperature, depth and distance to land (Figure 1, Table 2). The combination of the environmental tolerances represents the specific fundamental ecological niche for each of the three species. We found no difference between the three Asian species in any of the tidal variables (Figures 2A–C). All species experienced varied annual cycle deviation as well as cycle durations of varied lengths. The mean minimum sea surface temperature found for T. tridentatus was 17.0°C and hence 10 degrees lower than those for T. gigas and C. rotundicauda (Figure 2D). T. tridentatus experienced a larger mean temperature range of 12.7°C compared to T. gigas and C. rotundicauda, which experienced a mean range of 3.9 and 4.0°C, respectively (Figure 2E). The water oxygen levels followed the same pattern as sea surface temperatures, with a higher mean for T. tridentatus (5.0 mL/L) compared to T. gigas and C. rotundicauda (4.4 and 4.5 mL/L, respectively, Figure 2F). Neither the maximum chlorophyll A concentration nor the range of chlorophyll A concentration differed between the species (Figures 2G,H). We found a mean salinity of 33.5 PSU for T. tridentatus, which was slightly higher compared to T. gigas (31.8 PSU) and C. rotundicauda. (32.3 PSU, Figure 2I), although the narrow salinity range for C. rotundicauda was likely an artifact caused by inaccessible data from inland waters (see Discussion). The mean nitrate level was higher for T. tridentatus (1.9 μmol/L) than the other two species (0.6 μmol/L for both) (Figure 2J). Mean phosphate differed between all three species, with T. gigas experiencing the highest (0.28 μmol/L) and T. tridentatus the lowest (0.20 μmol/L) mean level (Figure 2K). T. tridentatus differed from the other two species in occupying a niche with higher mean pH value (Figure 2L).

FIGURE 1

Figure 1. Niche overlap analysis of the three Asian horseshoe crab species showing the principal component analysis plots of the environmental “hypervolume” for 23 predictor variables [CR, C. rotundicauda ( yes ); TG, T. gigas ( yes ); TT, T. tridentatus ( yes )]. Ellipses represent 68% of the hypervolume for each species; points represent presence of each species at environmentally unique locations. The positions of variable names have been adjusted to prevent overlapping.

FIGURE 2

Figure 2. Comparison of the three horseshoe crab species (CR, C. rotundicauda; TG, T. gigas; TT, T. tridentatus) showing the variation of 12 key environmental variables. Variable and unit are shown on the y-axis. Gray boxes indicate 50% of the sample points and are limited by the 1st (Bottom) and 3rd quartiles (Top). Black, horizontal lines within the gray boxes display the medians. Comparisons between C. rotundicauda and the two Tachypleus species for salinity were excluded because of absence of environmental data for inland waters. Asterisks show significant (<0.05) differences between species.

Species Distributions

The sampled distribution of C. rotundicauda includes the coasts of eastern India and Bangladesh, the coasts of Malaysia, Thailand, Cambodia, and southern Vietnam (Figure 3), but excludes the eastern coast of Vietnam, however, the distribution continues along the coast of northern Vietnam and southern China. C. rotundicauda is also present along the coasts of the Indonesian islands Sumatra and Java, as well as Borneo. The distribution of T. gigas is similar to that of C. rotundicauda, except that T. gigas seems to be absent from the coasts of Bangladesh, Cambodia, northern Vietnam, and China (Figure 4). T. tridentatus is distributed along the coasts of southern Japan, China, Taiwan, and northern Vietnam, in addition to the Philippines, and the islands of Borneo and Java (Figure 5).

FIGURE 3

Figure 3. Distribution maps for Carcinoscorpius rotundicauda with environmentally unique points, indicating the sampled distribution range for the species. Upper map (A) shows present-day projections of suitable habitat in red based on consensus of Maxent and SVM algorithms. Lower map (B) shows changes between the present-day and 2050 projections of suitable habitat based of consensus of Maxent and SVM algorithms, with green cells indicating loss of suitable habitat and red cells indicating gain of suitable habitat. Support values for the individual models are given in Table 3.

FIGURE 4

Figure 4. Distribution maps for Tachypleus gigas with environmentally unique points, indicating the sampled distribution range for the species. Upper map (A) shows present-day projections of suitable habitat in red based on consensus of Maxent and SVM algorithms. Lower map (B) show changes between the present-day and 2050 projections of suitable habitat based of consensus of Maxent and SVM algorithms, with green cells indicating loss of suitable habitat and red cells indicating gain of suitable habitat. Support values for the individual models are given in Table 3.

FIGURE 5

Figure 5. Distribution maps for Tachypleus tridentatus with environmentally unique points, indicating the sampled distribution range for the species. Upper map (A) shows present-day projections of suitable habitat in red based on consensus of Maxent and SVM algorithms. Lower map (B) show changes between the present-day and 2050 projections of suitable habitat based of consensus of Maxent and SVM algorithms, with green cells indicating loss of suitable habitat and red cells indicating gain of suitable habitat. Support values for the individual models are given in Table 3.

Distribution of Suitable Habitat

Both models, SVM and Maxent, demonstrated high predictive power, as seen by the high AUC values (all values >0.9) (Table 3). The Maxent model performed marginally better than SVM for all species except T. tridentatus. Omission errors were very low for both models (<3.5%), except for Maxent's future prediction for C. rotundicauda (6.2%).

Suitable habitats for C. rotundicauda are distributed along the coastline of mainland Southeast Asia, from the Bay of Bengal, including Sri Lanka, to southern China, except from the eastern coasts of Peninsular Malaysia and Vietnam (Figure 3A). Suitable habitats are also present in the Malacca Strait between Sumatra and Peninsular Malaysia, Southeast Sumatra, southern coast of Java, the coast of Kalimantan (Indonesian Borneo), and the coast of Sarawak (Malaysian Borneo). Additionally, small areas of suitable habitat are found in Sabah (Malaysian Borneo), West Papua, and in Manila Bay, Philippines. The distribution of suitable habitat differs from the sampled distribution, as it includes the coasts of Myanmar, southern Borneo, eastern Sumatra, and southern Java.

The distribution of suitable habitat for T. gigas is similar to that of C. rotundicauda, however, for T. gigas it does not go further east along the Southeast Asian mainland than to southeast Vietnam (Figure 4A), nor does it include the south coast of Java. Compared to the current known distribution, the distribution of suitable habitat includes the coast of Myanmar, Bangladesh, east coast of Sumatra, and southern Borneo.

The distribution of suitable habitat of T. tridentatus spans from northern Vietnam, along the coast of China, and into the Sea of Japan (Figure 5A). It also includes the coast of northwest Taiwan and Japan, including the Seto Inland Sea. Small areas of suitable habitat are present along the coasts of Southwest Myanmar, East Sumatra, South Java, Borneo, South Vietnam, Bay of Bangkok, and Manila Bay.

Changes in the Distribution of Suitable Habitat Under Climate Change Scenarios

We included three environmental variables, mean sea surface temperature, mean chlorophyll A concentration, and distance to land, for the modeling of the future distributions of suitable habitat (Table 2). Our models predicted that the habitat suitability for C. rotundicauda will increase in several coastal areas throughout Southeast Asia. More suitable habitats will arise in more eastern regions; the Philippines and along the coasts of several Indonesian islands: Sulawesi, Lombok, Sumbawa, and Flores (Figure 3B). Similarly, for T. gigas several areas along the coasts of Southeast Asia will become more suitable in the future (Figure 4B). The Bay of Bengal and North Sumatra will be less suitable. T. tridentatus will lose suitable habitat in the future, however, none of the lost areas are overlapping with the species' current occurrences (Figure 5B). In summary, none of the Asian horseshoe crabs seem to undergo a major range shift driven by prospective changes of the environmental factors tested here.

Discussion

Carcinoscorpius rotundicauda

We found that C. rotundicauda occurs along the coasts of Bangladesh, Cambodia, Vietnam, East Thai-Malay Peninsula, and China, in addition to the countries and coasts already surveyed by Sekiguchi (1988) (India, West Thai-Malay Peninsula, Singapore, and Indonesia). However, contrary to Sekiguchi (1988) we have not found evidence for its presence at Palawan, Philippines (Figure 3). According to the occurrence record, C. rotundicauda is absent along the coast of Myanmar, even though our models indicate suitable habitat in this area. Indeed, there is evidence that horseshoe crabs are present in this region, however, the available information points to L. polyphemus, and is likely a misidentification (Thapanand-Chaidee et al., 2010). We thus conclude that C. rotundicauda probably is present along the coast of Myanmar, and that missing records of horseshoe crabs are due to lack of studies in the area. The sampled distribution of C. rotundicauda shows a gap along the east coast of Vietnam, extending approximately 1,500 km in the high-resolution projections (based on Bio-Oracle variables, Figure 3) and approximately 400 km in the low-resolution projections (based on AquaMaps variables). The absence of the species in this region might be explained by the narrow continental shelf present along the coastline, which provides less horseshoe crab habitat, as well as water characterized by a high degree of mixing with deep water (Rojana-anawat et al., 2001). The distribution gap could also be related to food availability, as the east coast of central Vietnam has a lower macrozoobenthos biomass compared to the southern and northern coasts (Nguyen and Dao, 1995). However, the lack of mangroves along the east coast of central Vietnam may be the most plausible explanation for the distribution gap (Hong and San, 1993). C. rotundicauda is commonly known as the mangrove horseshoe crab, as it is often found in habitats characterized by mangroves and mudflats in contrast to the two Tachypleus species, which can be found at open sandy beaches (e.g., Cartwright-Taylor et al., 2011; Robert et al., 2014; Jawahir et al., 2017), thus a scarcity of mangroves would have a negative effect on the presence of C. rotundicauda.

When interpreting our model results, it is important to mention that species distribution models do have caveats, which may lead to false or uncertain projections of actual species distribution ranges (Jarnevich et al., 2015). We have tried to mitigate such caveats and uncertainties through an unbiased sampling of the biological data from independent sources, a careful selection of predictor variables, and through modeling with several environmental data sets and analytical algorithms. Nevertheless, the indication of a gap in the distribution of C. rotundicauda along the Vietnamese east coast should still be treated as a hypothesis that needs further validation with additional sampling and modeling.

The future distribution of suitable habitat of C. rotundicauda includes several coastal areas of Indonesian islands east of the sampled distribution. Since horseshoe crabs reside on the continental shelves (Sekiguchi and Shuster, 2009), it has been argued that the distributions of the Asian horseshoe crabs are delimited by the eastern edge of the continental Sunda Shelf (roughly corresponding to Wallace Line). Thus, even though areas of future suitable habitat are present east of the Sunda Shelf, we would not expect C. rotundicauda to be able to reside here, unless humans have mediated their dispersal. The models used in this study find an approximation of the fundamental niche of the horseshoe crabs, not the realized niche, thus geographical barriers are not necessarily detected.

We found that C. rotundicauda, as well as the other two Asian horseshoe crab species, inhabit coastal areas with varied tidal regimes. This is in contrast to L. polyphemus, as this species can be found in areas with mixed semi-diurnal tides (i.e., tidal cycles of approximately 12 h) (Rudloe, 1985; Barlow et al., 1986). The narrow range of salinity found for C. rotundicauda is most likely an artifact, as we know that the three Asian species occur in river estuaries, and C. rotundicauda has been found at salinities as low as 10 PSU (Chatterji, 1999; Mishra, 2009; Chen et al., 2015). Unfortunately, we could not include occurrences from inland rivers in Southeast Asia in our models, due to lack of environmental data from these areas. The narrow temperature range of C. rotundicauda reflects its tropical distribution, and can be related to its preference for temperatures between 20 and 30°C (Srijaya et al., 2014), which seems to be necessary for regular ecdysis (Lee and Morton, 2005). We find that C. rotundicauda is present in waters with lower pH levels compared to seawaters of normal pH level (around 8.2). This corresponds well with previous findings suggesting that pH levels found in horseshoe crab habitat are low due to the presence of estuaries (mixing of freshwater with pH levels of 7–7.5), photosynthesis by algae, as well as anoxic decomposition of organic material (Jawahir et al., 2017).

Tachypleus gigas

The distribution of T. gigas found in this study is consistent with that found by Sekiguchi (1988), with the exception that we also found occurrences of T. gigas along the east coast of the Thai-Malay Peninsula (Figure 4). Our high-resolution projections show limited suitable habitat at the southern extension range of the species (i.e., around Java), but the low-resolution projections clearly indicate suitable habitat along the northern coastline of Java.

As with C. rotundicauda, we suggest that T. gigas is present along the coast of Myanmar even though occurrence data are missing, due to the lack of research in the area. The presumed presence of T. gigas in this area needs to be confirmed with real observations in the future. Alternatively, there is the possibility that other factors not included in the models, such as biological interactions or substrate conditions, may actually refrain T. gigas form inhabiting this region. Our models predict additional future suitable habitats in more eastern regions of Indonesia, but again we argue that T. gigas is unable to disperse into these areas because of the deep waters east of the continental Sunda Shelf.

T. gigas is known to occur in areas with lower salinities than our results show (15–26 PSU) (Cartwright-Taylor et al., 2011). However, T. gigas is not as tolerant to low salinities as C. rotundicauda (Chatterji, 1999), and it has been shown that salinities lower than 20 PSU prevent its eggs from hatching (Zaleha et al., 2011). Apart from salinity, the ecological tolerances of T. gigas found in this study are similar to those of C. rotundicauda, which reflect their overlapping distributions.

Although we found overlapping niches and distributions for C. rotundicauda and T. gigas, the local habitat preferences are different between the species on a smaller spatial scale, i.e., C. rotundicauda is associated with mangroves and mudflats and T. gigas with sandy sediments at more exposed beaches, and adding an environmental factor such as sediment characteristics to our models might separate the niches of the two species (Jawahir et al., 2017).

Tachypleus tridentatus

In contrast to the T. tridentatus distribution reported by Sekiguchi (1988), we found no recent occurrence data on this species at the coasts of Sumatra, Sulawesi, South Vietnam, or Luzon, Philippines (Figure 5). However, we found several occurrences along the coast of South China. Although the coast of Myanmar is an area of suitable habitat for T. tridentatus, we would not expect it to occur here, as the area lies far west from any sampled occurrence sites. But according to our high-resolution model, the species might occur at the coast of Sumatra, as suggested by Sekiguchi (1988), even though we did not find any recent records of its presence here. There are no occurrences of T. tridentatus within the Yellow Sea, and we were unable to find historical evidence suggesting previous presence here, although the species is distributed at the coasts of the Korean Peninsula and Japan as well as in the East and South China Sea. Our high-resolution models indicated suitable habitat in the Yellow Sea, but this was not confirmed by the low-resolution models, leaving some degree of uncertainty in our projections for this area. The reason why T. tridentatus may be absent in the Yellow Sea could be related to a variety of factors, e.g., unsuitable sediment type or breeding habitat. If this discontinuity reflects a true distribution gap, it would suggest the presence of allopatric populations of T. tridentatus along each side of the Yellow Sea and future studies should aim to find more evidence for a distribution gap in this area. The future prediction of T. tridentatus shows less suitable habitat in several areas, however, none of these areas lie within the current known distribution of the species.

All Asian horseshoe crab species are known to occur at lower salinities than reported by this study, including T. tridentatus, as juveniles have been found at salinities of 8–24 PSU at a nursery beach in Hong Kong (Chiu and Morton, 2004). The temperature range of T. tridentatus is wider than of the two other Asian species, which reflects the wider latitudinal range of T. tridentatus, however, ecdysis in juveniles is halted at temperatures below 22°C (Lee and Morton, 2005). Thus, temperature is likely an important factor limiting the distribution of T. tridentatus to the north. Our results show that T. tridentatus experiences a high level and narrow range of pH, which might be explained by the interaction between temperature and pH. Lower temperatures result in higher pH (Zumdahl and Zumdahl, 2000), and thus the colder waters around Japan might have a higher pH compared to warmer waters surrounding more southern Asian countries. The chlorophyll A levels found for T. tridentatus does not differ from the two other Asian horseshoe crabs. Our results indicate that chlorophyll A concentration, in addition to temperature, is an important ecological factor determining the distributions of the Asian horseshoe crabs. This is supported by Hsieh and Chen (2009), who found that high densities of juvenile T. tridentatus were correlated with high chlorophyll A levels and high densities of polychaetes. They argued that chlorophyll A concentration reflected the amount of microalgae and hence the abundance of food available for higher trophic levels including polychaetes and that it was likely that polychaetes were an important food source for the juvenile horseshoe crabs. Thus, food resources for juvenile T. tridentatus increase with chlorophyll A concentration.

Potential Areas for Protection of Horseshoe Crabs

We suggest four criteria for identifying areas apt for conservation of horseshoe crabs in general, inspired from previous conservation studies (Petit et al., 1998; Hannah et al., 2007; Botero-Delgadillo et al., 2012). Marine Protected Areas (MPAs) with focus on horseshoe crab conservation should be prioritized in (1) geographic regions with disconnected distribution ranges, indicating physically distinct (i.e., allopatric) populations; (2) regions separated by large population breaks known from previous population genetic or phylogeographic studies; (3) areas with overlapping distributions of two or more species of horseshoe crabs; and (4) areas where the species experience little loss of suitable habitat under future climate scenarios. As previously noted, the potential distribution of suitable habitat, as well as the future potential distribution, found in this study should be considered as hypotheses of where Asian horseshoe crabs can be found and not as definitive results (Jarnevich et al., 2015). Thus, further investigations and models of the relationship between the species and the areas suggested here should be conducted prior to making decisions regarding conservation.

Our results show a wide gap in suitable habitat between C. rotundicauda populations in the southern and northern parts of the Vietnamese coast (Figure 3A). This, as well as the lack of mangrove habitat in the region (Hong and San, 1993), suggests the existence of allopatric populations of C. rotundicauda on the southern and northern coasts of Vietnam. Additionally, there is a gap in the distribution of suitable habitat at the northeast coast of Peninsular Malaysia, and there is evidence that genetically distinct populations of C. rotundicauda reside on each side of the Thai-Malay Peninsula (Obst et al., 2012; Adibah et al., 2015). The suitable habitat distribution also shows two separate distribution areas of C. rotundicauda on Borneo and Java, which could be allopatric to the three distinct distribution sites at mainland Southeast Asia. Consequently, based on our models, we hypothesize that there are probably five separated distribution ranges of C. rotundicauda in Southeast Asia, i.e., India to South Peninsular Malaysia, East Thailand to South Vietnam, North Vietnam to South China, Java, and Borneo, and none of these areas will become less suitable in the future according to our study. Based on our criteria 1, 2, and 4, we suggest the establishment of MPAs in coastal areas of West and East Thailand, North Vietnam and South China, and Borneo, in order to conserve C. rotundicauda in Southeast Asia (Figure 6).

FIGURE 6

Figure 6. Map of Southeast Asia and Japan showing regions which are optimal for the establishment of MPAs in order to conserve the three Asian species of horseshoe crabs. Rectangles with diagonal lines and green color indicate areas where MPAs for horseshoe crab conservation are absent, whereas rectangles with vertical lines and purple color indicate areas where MPAs for horseshoe crab conservation are already established.

Contrary to C. rotundicauda, we do not expect allopatric populations of T. gigas to exist across the Thai-Malay Peninsula as there is a continuity of suitable habitat for T. gigas in this region (Figure 4A). This assumption is supported by the short genetic distance between T. gigas collected from Vietnam and the Andaman Sea (Obst et al., 2012). The contrast in habitat continuity between T. gigas and C. rotundicauda can be explained by different habitat preferences for the two species. C. rotundicauda is more confined to mangroves and river deltas, while T. gigas is more frequent in truly marine environments, which is likely to result in higher dispersal along the coast. However, there is evidence of limited dispersal between the T. gigas populations in Northwest and Southwest Peninsula Malaysia (Rozihan and Ismail, 2011). Thus, based on our results, we hypothesize three distinct distribution areas for T. gigas: India to South Vietnam, Java, and Borneo. Additionally, C. rotundicauda and T. gigas are co-occurring at several locations along the coasts of the Southeast Asian mainland, hence according to criteria 1, 3, and 4 we suggest establishment of MPAs at the locations along the Southeast Asian mainland as wells as Borneo (Figure 6).

For T. tridentatus our results show four distinct distribution ranges; Japan to South Korea, China to North Vietnam, Borneo to Philippines, and Java (Figure 5A). T. tridentatus has overlapping distribution ranges with the two other Asian species; T. tridentatus and C. rotundicauda co-occur in South Vietnam as well as in the Hainan and Guangxi region of the Chinese coast, and all three Asian species co-occur along the coasts of Borneo and Java. According to criteria 1, 3, and 4, establishment of MPAs should hence be recommended in Japan, North Vietnam and South China, and Borneo (Figure 6). In Japan, MPAs are already established (http://mpatlas.org/explore/), and some of these include important breeding sanctuaries for horseshoe crabs, such as the Saikai National Park in Nagasaki. For areas where C. rotundicauda occurs sympatrically with T. gigas or T. tridentatus, it is important to consider the different habitat preferences for the species, i.e., low-salinity mangroves and river estuaries for C. rotundicauda, and high salinity sandy beaches for T. gigas and T. tridentatus.

Based on our analysis we propose that at least five regions should have high priority for horseshoe crab conservation: four in Southeast Asia, based on geographically separated areas with suitable habitat as well as the co-occurrence of two or more species, and one in Japan for the protection of distinct populations of T. tridentatus (Figure 6). We recommend that future MPAs are examined on site to ensure their suitability and that specific habitats, i.e., mangroves, mudflats, or gentle slope beaches, are present in order to provide the appropriate habitat for the horseshoe crab species in question (e.g., Kwan et al., 2016), as some factors of potential importance for horseshoe crab habitat suitability, such as sediment characteristics and presence of anthropogenic structures, could not be included in our analysis. The work presented here provides an initial step in the conservation of the Asian horseshoe crabs species, a research area in need of increasing attention considering the importance of these coastal organisms to human subsistence and health.

Future Studies

Destruction of breeding habitats pose a significant threat to horseshoe crabs (e.g., Mishra, 2009; Zaldívar-Rae et al., 2009; Nelson et al., 2016) and adult horseshoe crabs are harvested on the beaches during spawning for commercial exploitation. Since our study focused on the distribution range of adult and juvenile horseshoe crabs, we can only draw few conclusions about the distribution and preferred environmental conditions of horseshoe crab breeding sites. However, the approach implemented in this study could also be used to inventory and predict potential breeding sites using scientific networks and fieldwork, but should then include additional environmental variables. Horseshoe crabs nest on beaches with slight slopes and in mangroves and estuaries. Hence data on habitat types, beach inclination, sediment type and granularity, wave action, and organic contents of the sediment are essential variables that need to be measured and included for such analysis. Additionally, our study includes tidal regime environmental variables in the ecological niche modeling that can be used in future research on the distributions of other intertidal organisms, and thus could be of importance for the conservation of intertidal communities at a global scale.

Conclusions

The distributions of the three Asian horseshoe crab species have been assembled with this work, and marine areas of suitable habitat have been identified. Additionally, our models predict that C. rotundicauda and T. gigas in the future will experience increased habitat suitability. The information on ecological variables important for the three species, as well as the current and future distributions of suitable habitat provided here, can be used as an initial step in determining where MPAs should be established for conservation purposes and to determine the IUCN conservation status of the three Asian horseshoe crab species.

Author Contributions

SV, PF, and MO: designed the study; SV, MO, II, and PF: collected the data; MO and FQ: performed the modeling and analysis; SV: led the writing of the first draft in close collaboration with PF and MO. All authors contributed substantially to the revisions.

Funding

This study is supported by the EU's Seventh Framework Program project BioVeL (www.biovel.eu, grant no. 283359), the Horizon 2020 program ENVRI+ (grant agreement no. 654182), the Swedish LifeWatch program funded by the Swedish Research Council (grant no. 829-2009-6278), and by the Danish National Research Foundation Niels Bohr professorship project Aarhus University Research on the Anthropocene (AURA).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Dr. Shin Nishida, Dr. Hwey-Lian Hsieh, and Joseph Klein for valuable feedback and help to improve this paper. SV and PF consider this work a contribution to the Danish National Research Foundation Niels Bohr professorship project Aarhus University Research on the Anthropocene (AURA).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2018.00164/full#supplementary-material

Supplementary Data Sheet 1. Explanation of the construction of global tidal variables.

Supplementary Image 1. Heat map for C. rotundicauda using Maxent.

Supplementary Image 2. Heat map for C. rotundicauda using SVM.

Supplementary Image 3. Heat map for T. gigas using Maxent.

Supplementary Image 4. Heat map for T. gigas using SVM.

Supplementary Image 5. Heat map for T. tridentatus using Maxent.

Supplementary Image 6. Heat map for T. tridentatus using SVM.

Supplementary Table 1. Data provider names for occurrence records obtained by own observations, collaborators, and scientific network.

Supplementary Table 2. Data references for occurrence records obtained by literature.

Supplementary Table 3. Data provider names for obtained GBIF occurrence records. The name (or acronym) in use by the institution having custody of the object(s) or information referred to in the record.

Supplementary Table 4. Test statistics for analysis of differences in environmental conditions between four horseshoe crab species using Welch's test for unequal variances.

Supplementary Table 5. Multiple comparisons of difference between means of environmental variables for four horseshoe crab species (CR – C. rotundicauda, TG – T. gigas, TT – T. tridentatus) using Games-Howell Post-Hoc test. Significant mean difference at the 0.05 level in bold.

Supplementary Table 6. Loadings for the different environmental variables shown for all PCA components for C. rotundicauda.

Supplementary Table 7. Loadings for the different environmental variables shown for all PCA components for T. gigas.

Supplementary Table 8. Loadings for the different environmental variables shown for all PCA components for T. tridentatus.

References

Adibah, A. B., Ng, W. L., and Tan, S. G. (2015). The Malay Peninsula as a barrier to gene flow in an Asian horseshoe crab species, Carcinoscorpius rotundicauda Latreille. Biochem. Syst. Ecol. 60, 204–210. doi: 10.1016/j.bse.2015.04.026