The two Asian horseshoe crab species, C. rotundicauda and T. tridentatus, inhabit different microhabitat and sites in Hong Kong, where C. rotundicauda is usually found in mudflats near mangal habitats whilst T. tridentatus is found in sandier substratum, suggesting potential differential prey preference and/or environmental requirements. Samples of C. rotundicauda were collected from Wu Shek Kok, Sha Tau Kok, Hong Kong (22°32′N, 114°12′E) on September 10 and October 23, 2019, and samples of T. tridentatus were collected from Ha Pak Nai, Yuen Long, Hong Kong (22°25′N, 113°56′E) on September 11 and 24, 2019 (Figure 1). These sampling sites were chosen based on previous study of horseshoe crab populations in Hong Kong of which showed that the largest population of juvenile C. rotundicauda and T. tridentatus could be found in Wu Shek Kok, Sha Tau Kok, and Ha Pak Nai, Yuen Long, respectively. Wu Shek Kok locates on the north-eastern coastline and contains a large patch of dwarf mangrove (Figure 1). This mangrove supported various marine organisms, especially gastropods and bivalves. Ha Pak Nai is a sandy shore in Deep Bay at Northwest Hong Kong. This site is rich in seagrasses that provide nutrients for many marine fauna and oysters were commonly found (Figure 1; Kwan et al., 2016). Horseshoe crabs were collected by hand during low tide and were individually placed in plastic boxes with filtered seawater. They were allowed to stand for 1–3 h until egestion. The feces were collected using a dropper, and preserved in 95% ethanol immediately upon collection and stored in 4°C before DNA extraction. Scaled photos of each horseshoe crab were taken after collection of fecal samples for measurement of six morphometric parameters (i.e., the widths and lengths of prosoma and opisthosoma, paired dorsal compound eye distance and telson length; Supplementary Figure 1) in order to estimate their instar stages based on previous laboratory culture experiment (Chen et al., 2010). The horseshoe crabs were released back into their natural habitats after the feces collection. In each trip, five in situ sediment samples (5 ml in volume each) were collected at the area that horseshoe crabs were collected, to examine the prey availability in the habitat by eDNA analyses. The sediment samples were kept on ice during transportation back to the laboratory and were pooled and mixed well with a spatula to obtain one 25 ml pooled sediment sample in each sampling trip. As a result, a total of four pooled sediment samples were included in the present study (two trips for each of the two sites, Wu Shek Kok and Ha Pak Nai). All pooled sediment samples were stored at −80°C before further treatment.

DNA was extracted from the feces and sediment using QIAamp PowerFecal kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol. The same extraction kit was applied on two different types of samples to minimize discrepancies during DNA extraction and allow more reliable comparisons on metabarcoding data obtained. Five DNA extraction replications were carried out for each pooled sediment sample. The nuclear 18S rRNA gene was examined because this marker can be readily amplified across a broad spectrum of food items from different phyla, including algae which were previously suggested to be a major carbon sources for horseshoe crab species (Kwan et al., 2015; Fan et al., 2017). The V4 region of the 18S (∼450 bp) was amplified using the eukaryotic universal primers 3NDF (5′-GGCAAGTCTGGTGCCAG-3′) and 18S ER1 (5′-GACTACGACG-GTATCTRATCRTCTT-CG-3′) (Bråte et al., 2010). To prevent the amplification of host (horseshoe crab) DNA, a horseshoe crab specific peptide nucleic acid (PNA) probe HSC_18S (5′-NH2-GGCGGTTACTTCCTGGC-Lys-3′) was developed based on the sequence region unique to horseshoe crabs identified from the alignment of 18S sequences of the four horseshoe crab species and other related animals (e.g., spider and scorpion) downloaded from National Center Biotechnology Information (NCBI) GenBank. The PNA probe acted as a blocking primer that differentially bound to the sequence of horseshoe crab and inhibited the amplification.

The PCR amplification of each fecal sample was carried out with a reaction mix containing 7.5 μL of template DNA, 1× PCR buffer, 1 mM dNTPs, 10 nmol of each primer, 20 μg bovine serum albumin (BSA) (New England Biolabs, Ipswich, United States), 12.5 nmol of PNA probe, 0.25 U TaKaRa Taq^TM DNA Polymerase (Takara Bio Inc., Shiga, Japan) and sterile distilled water to a total volume of 25 μL. The PCR amplification of each sediment sample followed the same protocol, except no PNA probe was added and template DNA amount was reduced to 2.5 μL. Positive (with horseshoe crab genomic DNA) and negative (with sterile distilled water) controls were included in each batch of PCR reaction to validate the efficiency of PNA probe and reveal any occurrence of false-positive (contamination) and false-negative. The thermal cycle consisted of a 2-min initial denaturation at 94°C, 40 cycles of 15 s at 94°C and 2-min at 63°C, and a finial extension of 3 min at 72°C. Three PCR replicates were done for each sample to obtain a 75 μL PCR product to minimize PCR bias or randomness, and were purified by Gel/PCR DNA fragments kit (Geneaid, Taiwan, China) following the manufacturer’s protocol. All the purified PCR amplicons were sent to commercial company GENEWIZ (Suzhou, China) for high-throughput sequencing using Illumina MiSeq. The purified PCR amplicons were ligated with adaptors incorporated with a 6-bp barcoding region for multiplexing the samples from different locations and species. Preparation of libraries, quality assessment, and paired-end 2 × 300 bp sequencing on Illumina MiSeq platform were conducted by the commercial company.

Demultiplex of raw reads were done by GENEWIZ. Paired-end raw reads of each sample were assembled using the standard bioinformatics analysis pipeline of QIIME2 (version: qiime2-2019.10; Bolyen et al., 2019), DADA2 plugin was used for quality and sequence-length filtering, and removal of chimera and singletons. Similarity threshold of 97% was adopted for taxonomic assignment of operational taxonomic units (OTUs). Taxonomic assignment of all the OTUs representative sequences were done by both QIIME-compatible 18S SILVA SSU r119 database (Yilmaz et al., 2014) and NCBI GenBank. Sequences with ≥97% similarity were assigned to class-level, for sequences exhibiting similarity results lower than 97% but higher than 80% were assigned to phylum-level. Sequences were classified as unidentified taxa when exhibiting similarity results lower than 80% or no match was found in both databases. OTUs with less than 1% of the total read in each sample and sequences of horseshoe crabs were removed before statistical analyses. All the common OTUs were assigned into different taxonomic groups based on the classification results.

The relative read abundance (RRA) were used as a proxy of (1) proportion of intake by the horseshoe crabs for fecal samples and (2) proportion of prey availability in the habitats for sediment samples, while frequencies of occurrence (%FOO) of different comestibles were calculated to reflect the prevalence of ingested taxa among individuals. Non-metric multi-dimensional scaling (NMDS) and permutational multivariate analysis of variance (PerMANOVA) were performed in RStudio (RStudio Team, 2020) based on the Bray–Curtis dissimilarity matric, using the “ordinate”, “adonis”, and “ggplot2” functions from R packages, to visualize and to examine the differences in (1) dietary compositions between two species, (2) dietary preferences among instar stages, (3) prey availability between two study sites, and (4) dietary compositions and the prey availability in their habitat at the taxonomic group level. The unidentified taxa were excluded in both NMDS and PerMANOVA analyses, as the sequences could be originated from a wide range of organisms that was inappropriate to categorize as a single group. NMDS and PerMANOVA were performed on both RRA and %FOO metrics. Two matrixes were retained for the NMDS: (1) square-root transformed data of RRA or %FOO and (2) sampling sites, instar stages, and sample types.

PCR amplifications of the 18S were succeeded in a total 107 DNA samples extracted from 87 fecal samples (46 from C. rotundicauda and 41 from T. tridentatus) and 20 sediment samples (5 replicates from each of the 4 pooled sediment samples). The 87 horseshoe crabs of both species ranged from 6th to 9th instars according to Chen et al. (2010; Table 1). A total number of 2,164,120 filtered reads (mean = 20,225 ± 7,675) were obtained from 107 PCR products, including 942,756 reads from 46 C. rotundicauda fecal samples (mean = 20,495 ± 3,880); 1,011,472 reads from 41 T. tridentatus fecal samples (mean = 24,670 ± 7,413); 114,834 reads from 10 Wu Shek Kok sediment samples (mean = 11,483 ± 5,129); and 94,929 reads from 10 Ha Pak Nai sediment samples (mean = 9,493 ± 6,710). All paired FASTQ files are available in NCBI BioProject (PRJNA758222).

The final OTUs table consisted of 401 OTUs from 21 phyla, assigned from the 107 samples. A total of 160 OTUs (16 phyla) were detected in C. rotundicauda fecal samples; 114 OTUs (14 phyla) were detected in T. tridentatus fecal samples; 62 OTUs (6 phyla) were detected in Wu Shek Kok sediment samples; and 129 OTUs (16 phyla) were detected in Ha Pak Nai sediment samples. A total of 23 OTUs were common in the fecal samples of both species; 16 OTUs were unique in Wu Shek Kok sediment and C. rotundicauda fecal samples; 15 OTUs were unique in Ha Pak Nai sediment and T. tridentatus fecal samples; no OTU was common only in both sediment samples from Wu Shek Kok and Ha Pak Nai; no OTU was shared among all the samples (Figure 2). Furthermore, three OTUs were unique in fecal samples of both species and Wu Shek Kok sediment samples; two OTUs were unique in C. rotundicauda fecal samples and both sediment samples from Wu Shek Kok and Ha Pak Nai, and no OTU was unique in other relations. Of the total number of OTUs detected, 360 could be assigned down to the class or subclass level, and all of the OTUs were assigned to 14 taxonomic groups based on the classification results (Supplementary Table 1). The sequencing depth was sufficient enough to detect the OTUs occurring in each sample, as reflected by the plateau in rarefaction curves (Supplementary Figure 2).

The results from the two sampling days for each of the two species were largely consistent so the data from the two trips were analyzed together since the sample size in the present study was not sufficient for revealing any temporal pattern. The analyses of RRA showed that oligochaetes were the major prey group for both of the two horseshoe crab species in Hong Kong, which contributed to the largest proportion of reads in both C. rotundicauda (41.6%) and T. tridentatus (32.4%) (Figures 3A, 4). These results also matched with the percentage of frequencies of occurrence (%FOO; Figure 5, showing that oligochaetes could be recovered in 91.3% of C. rotundicauda fecal samples and 82.9% of T. tridentatus fecal samples. Bivalves and crustaceans were other prey groups that contributed to a relatively high proportion of reads (8.6 and 8.4%, respectively) in the diet of C. rotundicauda, where bivalves were recovered in 21.7% of C. rotundicauda individuals and crustaceans were found in 52.2%. Surprisingly, a relatively high proportion of anthozoans (22.8% of RRA) was detected in the diet of T. tridentatus, and it was found in 65.9% of individuals analyzed. This resulted in a significant difference between the dietary compositions of C. rotundicauda and T. tridentatus, both in term of RRA and %FOO (p < 0.001, Table 2; NMDS, Figure 6). One noteworthy observation was that the difference of diet compositions between different instar stages is significant in C. rotundicauda (RRA and %FOO: p < 0.01) but not in T. tridentatus (RRA and %FOO: p < 0.5, Table 2 and Figure 3B). Additionally, SAR (stramenopiles, alveolates, and rhizarians) contributed 24.4% of reads in C. rotundicauda and 13.3% in T. tridentatus (Figures 3A, 4), where SAR was recovered in 91.3% of C. rotundicauda diet and 73.2% of T. tridentatus (Figure 5). However, these were probably a result of incidental consumption or secondary predation, rather than both species were actively feed on SAR. Apart from SAR, unidentified taxa contributed 1.2% of reads in C. rotundicauda and 13.3% in T. tridentatus (Figures 3A, 4), where unidentified taxa were recovered in 13% of C. rotundicauda diet and 78% of T. tridentatus (Figure 5).

Bivalves and gastropods were the most abundance prey groups in Wu Shek Kok (51.3 and 38.3%, respectively) as shown by eDNA analyses on the sediment samples, while other metazoan (including taxa from platyhelminthes, gastrotrichs, cnidarians, kinorhynchs, and rotifers) and SAR contributed most reads in Ha Pak Nai sediment samples (29.4 and 26.0%, respectively; Figures 3C, 4). There was a significant difference in prey availability between the two sampling sites (RRA and %FOO: p < 0.001, Table 2), this result was also reflected by NMDS (Figure 6) and the analysis of RRA (Figures 3C, 4). More importantly, the differences between the dietary composition of C. rotundicauda and the prey availability in Wu Shek Kok were significant (RRA and %FOO: p < 0.001, Table 2). It was the same case between diet composition of T. tridentatus and the prey availability in Ha Pak Nai (RRA and %FOO: p < 0.001, Table 2). Less than 10% of reads from the major prey items of both species (i.e., oligochaetes, crustaceans, and anthozoans) were detected in sediment samples, indicating that both horseshoe crab species are selective feeders.

Our DNA metabarcoding results show that juveniles of C. rotundicauda and T. tridentatus are carnivores which mainly consume benthic animals such as oligochaetes, bivalves, and crustaceans. This finding is largely consistent with previous knowledge of prey consumption detected using morphological examinations and stable isotope analyses (Zhou and Morton, 2004; Pahang et al., 2012; Kwan et al., 2015). Oligochaete worms often constitute an important proportion of the meiofauna from estuarine mangrove sediments in subtropical and tropical waters (Erséus, 2002), so it is not surprising that they can serve as major food source for juvenile horseshoe crabs. However, our data did not find any evidence of seagrass as food sources for juveniles of T. tridentatus and C. rotundicauda in Hong Kong and supported previous studies using stable isotope analyses (Kwan et al., 2015; Fan et al., 2017). Although algal reads were detected in the fecal samples of some individuals (5 out of 87 individuals: 3 C. rotundicauda and 2 T. tridentatus), those reads only occupied a small proportion (<0.5%) of the total reads (Figures 3A, 4). Thus, it is probably a result of incidental consumption or secondary consumption, rather than evidence for omnivorous nor herbivorous behaviors for the crabs. Although the possibility of technical artifacts (e.g., differential digestibility and PCR amplification efficiency amongst food items) and/or seasonal variations in dietary composition in resulting failure of recovery of seagrass cannot be ignored, we justified it was unlikely in our study. The animal tissues are generally easier to be digested compared to algal tissues and we could detect a considerable amount of algal DNA from sediment samples collected from Ha Pak Nai, suggesting the technical artifact is not significant in this circumstance. On the other hand, our study was conducted in late summer, which overlapped with the studied period by Kwan et al. (2015) and several species of seagrass were observed in the Ha Pak Nai during our sampling trip, in consistent with the results of eDNA analyses. The carbon values from the seagrasses reported in previous stable isotope studies might be originated from the preys of horseshoe crabs as many marine meiofauna feed on seagrass detritus. Therefore, our data confirmed with previous studies that juveniles of C. rotundicauda and T. tridentatus were carnivores (Zhou and Morton, 2004; Pahang et al., 2012), which also matched with the morphology of their mouth that gnathobases (chitinous bristles) are adapted for holding and shredding the tissues of benthic animal (Wyse, 1971; Hayes and Barber, 1982; Razak and Kassim, 2018).

Remarkably, our molecular approach enabled to detect unexpected interactions between horseshoe crabs and other marine taxa. We showed that the sea anemone was one of the major prey groups of juvenile T. tridentatus, which has never been reported in any other previous studies. Several species of burrowing anemone, including Paracondylactis hertwigi and Edwardsia japonica, are commonly found in soft shores in western water of Hong Kong, including Ha Pak Nai where the juveniles of T. tridentatus were collected in the present study (Morton, 1983). The closest match for the Anthozoa 18S sequences detected in the horseshoe crab metabarcoding data was Edwardsia (>98%), suggesting that the detected anthozoan likely belongs to the family Edwardsiidae which consisted predominantly sediment burrowing anemone. This provides further evidence that the horseshoe crabs preyed on the burrowing anemone in the sandy habitat where the anemone was present. The anemone was only detected in a single individual of C. rotundicauda probably because they are rare in mangal habitat. Similarly, previous dietary studies on other marine taxa also reported that DNA metabarcoding enable us to reveal cryptic predation relationships on jellyfish and other cnidarians [e.g., jellyfish is a major prey group of Adélie penguin (Jarman et al., 2013); black-browed and Campbell albatross selectively feed on jellyfish (McInnes et al., 2017)], which could be easily overlooked by gut microscopic examinations and stable isotope analysis. Thus, our results recapitulate the omitted role of cnidarians in many food webs and provide new insights into the energy flow in the estuarine mudflat ecosystems.

Our results provide further evidence that the two horseshoe crab species are selective feeders as suggested by gut content analyses (Chiu and Morton, 2004; Zhou and Morton, 2004; Pahang et al., 2012). If horseshoe crabs are non-selective opportunistic feeders that consume whatever available in the habitat, we would expect their diet composition would closely mirror the benthic composition in the sediment. Yet comparison of metabarcoding results from the fecal samples and eDNA analyses on sediment showed that oligochaetaes were the major food items for both two species, but the relative abundance of oligochaetes was low in both habitat (<1% of RRA), indicating significant discrepancies between the ingested food and the prey availability in their corresponding habitats. Therefore, our findings are consistent with previous studies that suggest juvenile of both C. rotundicauda and T. tridentatus are selective feeders (Chiu and Morton, 2004; Zhou and Morton, 2004; Pahang et al., 2012), and provide no support to the hypothesis of generalist as proposed by stable isotope analyses (Kwan et al., 2015; Fan et al., 2017). Interestingly, although mollusks are commonly found as major prey for various horseshoe crab species, especially for the adult (e.g., Botton, 1984; Botton and Ropes, 1989; Chatterji et al., 1992; Pahang et al., 2012), we found the consumption of mollusks was uncommon in the juveniles of the two studied species, with only approximately 10% contribution to the diet in C. rotundicauda on average and only observed detected in 20% of the individuals in spite of mollusks were highly abundant in its habitat (approaching 90% in the eDNA analyses), and <10% in T. tridentatus (in term of both RRA and FOO). We anticipate the mouth part of juvenile horseshoe crab is not strong enough to crush the mollusk shells, so they prefer softer preys that are easier to manipulate such as annelid worms and sea anemone as shown here. This may imply a possible shift to more abundant and nutritious mollusks as major prey items once the horseshoe crabs grow up. This hypothesis on dietary change at different life stages remains to be verified in future studies.

Carcinoscorpius rotundicauda and T. tridentatus live in different habitats: C. rotundicauda is more common in muddy habitat in close proximity to mangrove with finer grain size, whilst T. tridentatus is usually found in more exposed bottom with larger grain size. This implies that dietary preference is anticipated to be one of the factors that lead to the separation in habitat between the two species. We found that the two species exhibit overlapping but significantly different diet compositions. Both species consume oligochaetes as the major food source that may reflect conservatism in functional morphology and ecological niche in closely related species. The major difference reported here emphasizes that C. rotundicauda consumes more crustaceans whilst T. tridentatus ingests a larger amount of anthozoans. However, these differences in dietary composition did not correlate with the prey availability as reflected by eDNA analyses, suggesting that establishment of habitat preference may not be based on the prey preference. This implies that the habitat requirement and ecological partitioning of horseshoe crab species may not be directly affected by prey availability, but other physical parameters, for instance sediment grain size is proposed to be an important factor that affect the ecological partitioning of C. rotundicauda and T. tridentatus (Li, 2008). Thus, habitat protection for both mangrove mudflats and sandy shores should be enforced in order to conserve both populations of C. rotundicauda and T. tridentatus.

This study demonstrated how DNA metabarcoding approach could be adopted to provide new insights into the feeding habit of horseshoe crabs, and hence to improve our knowledge on estuarine trophic ecology and species conservation. Admittedly, the sample size of the current study was relatively small but the non-invasive protocol by using fecal sample collection enables a standardized and accurate investigation that could be easily scaled up to verify current finding and answer further biological questions without the need of sacrificing the threatened animals. The PNA blocking primer developed in this study could also be applied for the other horseshoe crab species that are not included here. Therefore, we would anticipate the protocol can be more widely applied to conduct larger spatial and temporal studies. Furthermore, some evidences of dietary difference among instars were observed in C. rotundicauda, but not T. tridentatus. However, this may also be an artifact of small sample size analyzed for each instar stage in the present study, and more individuals from different instar stages are required for the evaluation of the intra-specific and inter-specific feeding habit to enhance our understanding on these iconic animals of high conservation value.

Future optimization on the protocol will focus on minimizing amplification bias, which is widely criticized as one of the major sources of error in translating the metabarcoding data into the dietary compositions (Alberdi et al., 2018). The inclusion of mitochondrial marker, in particular the cytochrome c oxidase subunit I (COI) that is considered as a barcoding gene for metazoans, can help to reduce the bias arisen from PCR amplification and provide better taxonomic resolution on the identities of prey below family level that cannot be achieved by the more conserved 18S. As aforementioned, we adopted the more conservative 18S marker in the present pilot study because we attempted to maximize the prey items that can be recovered and to test the hypothesis whether seagrass represents one of the main food items for horseshoe crabs. Based on the current data, we show that horseshoe crabs are likely to be selective carnivores, and thus future research will focus on more detailed prey identities down to species level to have a better resolution on the prey composition data. Accordingly, a local inventory sequence database for different fauna will also be necessary, which is not only applicable to horseshoe crabs, but also facilitate the metabarcoding dietary studies in other endangered and ecologically/economically important species. We believe that this technique could provide new directions for us to understand of the trophic relationships between organisms, and therefore facilitate future design of conservation strategies.