1 Introduction

The horseshoe crabs, as well-known classic 'living fossils', have survived for nearly 500 million years (Van et al., 2010). Their survival is inseparable from their feeding habits, morphological structure, and physiological and metabolic characteristics, including omnivory, strong immunity, and tolerance to environmental changes. Tachypleus tridentatus, Carcinoscorpius rotundicauda, Tachypleus gigas, and Limulus polyphemus are the four species of horseshoe crabs in the world (Wang et al., 2020). L. polyphemus and T. tridentatus have been captured since the 1980s to provide pharmaceutical manufacturers with Limulus/Tachypleus amebocyte lysate (LAL/TAL) by using their blood as the raw material (Novitsky, 2015). In the modern medical system, every injection, every intravenous infusion, and every medical device that needs to be implanted into the human body need the LAL/TAL reagent for endotoxin testing. L. polyphemus and T. tridentatus are listed as Vulnerable (VU) (Smith et al., 2016) and Endangered (EN) species (Laurie et al., 2019), respectively, on the International Union for Conservation of Nature (IUCN) Red List. High tide spawning sites and intertidal juvenile nursery habitats are critical to the horseshoe crab's survival (Laurie et al., 2019). The overharvesting and habitat loss with the development of coastal economic zones have severely affected the survival of horseshoe crabs. All four horseshoe crab species worldwide are now endangered (Wang et al., 2020). Both T. tridentatus and C. rotundicauda distributed in China have been strictly protected by law as National Second-class Protected Wildlife Species since February 2021. Recently, the significant population declines have gained in notoriety due to the integration of science, management, education, and policy for strengthening conservation research. Habitat protection has been considered one of the best strategies to conserve horseshoe crabs (Chen et al., 2015; Laurie et al., 2019).

All living animals depend on food for their survival. Feeding habits, digestion, and absorption are directly related to important life activities such as growth and reproduction. The horseshoe crabs can change their nutritional sources and upgrade their position in the food chain via growth (Gaines et al., 2002). The digestive tract of horseshoe crab is partially formed at hatching and becomes functional before feeding begins. The digestive tract plays a vital role during the digestion in the animals. The structure and function of the digestive tract of the horseshoe crab are closely related to its omnivorous feeding habit (Hong, 2011). The external morphology and anatomy of the digestive tracts of T. tridentatus and C. rotundicauda are similar (Chen et al., 2016). The basic internal gross anatomy of horseshoe crabs has been described. The digestive tract is divided into foregut, midgut, and hindgut three parts (Chatterji et al., 1988; Hong, 2011). The foregut includes the cavum buccale, esophagus, stomach/proventriculus, and pylorus. The hindgut consists of the rectum and anus (Chen, et al., 2016). There is no remarkable difference in structure between midgut and hindgut (Xie, 2018).

At present, many research reports have been published concerning the feeding habits and digestive metabolism of horseshoe crabs, generally focusing on their food sources and the anatomy of the digestive tract (Gaines, et al., 2002; Carmichael et al., 2009; Hong, 2011; Chen, et al., 2016; Fan et al., 2017). However, the mechanisms of feeding, digestion, and absorption in horseshoe crabs remain unclear, and information on the transcriptome of digestive organs for T. tridentatus and C. rotundicauda is limited. The transcriptome evolution of different species can be comprehended by collecting and comparing transcriptome data (Wang et al., 2021). In this study, high-throughput Illumina Solexa sequencing and gene annotation were used to characterize the transcriptome of the digestive tract and mass mucosa related to the digestive metabolism of T. tridentatus and C. rotundicauda. The first comparative analysis of large-scale gene expression profiles between T. tridentatus and C. rotundicauda juveniles is reported herein. Thus, the transcriptome analysis reported in this article is helpful for understanding the molecular mechanisms of the adaptive evolution related to the digestive metabolism of T. tridentatus and C. rotundicauda.

2 Methods 2.1 Horseshoe Crabs and Sample Preparation

Experimental animal samples were collected in accordance with the ethical principles of animal experimentation. The materials were stomach (S), midgut (G), and fecal mucosa (F) isolated from three individuals of T. tridentatus and C. rotundicauda, respectively. The animals were all healthy juveniles with 3.3 cm ± 0.5 cm of prosomal width. The animals were obtained from the intertidal zone in the west bay of northern Beibu Gulf, Guangxi, China at the end of March 2019. The juveniles were transported to the laboratory and acclimated to laboratory conditions for three weeks. All horseshoe crab juveniles were cultured in aquarium tanks (dimensions: 120 cm × 40 cm × 25 cm) equipped with a water filtration system, thermostatic heaters, and ultraviolet sterilizers. A 4-cm sediment layer was provided underneath. Seawater was maintained under the following conditions: temperature 26 – 30℃, salinity 32 – 33, pH 7.6 – 7.9, dissolved oxygen 6 – 7 mg L⁻¹, and a light: dark photoperiod of 12 h: 12 h. During the acclimation period, the horseshoe crab juveniles were fed a daily ration of 2% of their biomass with frozen brine shrimp and clam meat. The juveniles were randomly transferred from the acclimation tank to the aquaria used for the experiments and kept for two days without food until the beginning of the experiment. Three tissue samples from three horseshoe crab juveniles were collected, and each tissue sample was cut into small pieces and put into 2-mL RNase-free centrifuge tube. An appropriate volume of RNA Later (RNA Stabilization Reagent, QIAGEN) was added, and the sample was stored at −80℃. For easy identification of the samples, three stomach samples of T. tridentatus were named TTS1, TTS2, and TTS3; three midgut samples of T. tridentatus were named TTG1, TTG2, and TTG3, and three fecal mucosa samples of T. tridentatus were named TTF1, TTF2, and TTF3. The three stomach samples of C. rotundicauda were named CRS1, CRS2, and CRS3; the three midgut samples of C. rotundicauda were named CRG1, CRG2, and CRG3; and the three fecal mucosa samples of C. rotundicauda were named CRF1, CRF2, and CRF3.

2.2 RNA Extraction and Quality Control

Total RNA was extracted from the 18 tissue samples with the RNeasy MiNi Kit (QIAGEN, USA) in accordance with the manufacturer's instructions. The quality of RNA samples was controled with several aspects. First, RNA degradation and possible contamination were monitored by 1% agarose gel electrophoresis. Then, the purity, concentration, and integrity of total RNA were respectively measured using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), the Qubit ® RNA Assay Kit in Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA), and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

2.3 Library Preparation for Transcriptome Sequencing

Library preparation for transcriptome sequencing was performed by conventional methods of molecular biology. RNA sequencing samples were prepared based on 1.5 µg RNA per sample. Sequencing libraries were constructed using a NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, USA) according to the manufacturer's instructions. Briefly, oligo (dT)-attached magnetic beads were used to purify mRNA from total RNA. Library quality was assessed with the Agilent Bioanalyzer 2100 system. Libraries were sequenced by the Illumina Hiseq platform following the manufacturer's protocols, and paired-end reads were obtained.

2.4 Transcriptome Assembly

High-quality clean reads obtained by removing reads containing poly-N and reads containing Illumina adapters, and low quality reads from raw data were subjected to subsequent analysis based on the methods of quality control (Zhang et al., 2017). In the present study, at least 6 Gb clean data were obtained for each sample. Trinity's utility (Grabherr et al. 2011) was used to perform de novo transcriptome assemblies.

2.5 Data Analysis 2.5.1 Comprehensive gene annotation

We performed gene functional annotation using the following databases. The NCBI nonredundant protein sequences (NR) database was used for identifying putative mRNA functions. The Protein family (Pfam) database was used for classifying protein sequences into families and domains. The possible functional classifications and molecular pathways were predicted with databases including NCBI nonredundant nucleotide sequences (NT), Clusters of Orthologous Groups of proteins (KOG/COG), a highquality manually annotated and reviewed protein sequence database (SWISS-PROT), the KEGG Ortholog database (KO), and Gene Ontology (GO).

2.5.2 Differential expression analysis

For the samples in the present experiments, differential gene expression analysis of two experimental groups was conducted using the DESeq R package (1.10.1). The Benjamini and Hochberg approach for controlling the false discovery rate was applied to adjust the resulting P values. Genes with adjusted P-values below 0.05 were determined to be differentially expressed. GO enrichment analysis of the differentially expressed genes (DEGs) was conducted by the GOseq R packages (Young et al., 2010).

3 Results 3.1 Summary of Sequencing Data

High-quality sequencing and mapping results were obtained. The summary statistics of the sequencing data are shown in Table 1. After removing the low-quality reads, we obtained clean reads ranging from 40.37 million to 52.95 million, and clean bases ranging from 6.06 Gb to 7.94 Gb from the stomach, midgut, and fecal mucosa samples of C. rotundicauda. The number of clean reads from three experimental tissues of T. tridentatus ranged from 43.24 million to 56.13 million. Correspondingly, the number of clean bases ranged from 6.49 Gb to 8.42 Gb (Table 1). The average GC content of the transcriptome data of T. tridentatus stomach (TTS) was 37.36%, a value slightly lower than that of C. rotundicauda which was 40.74%. The average GC contents of the transcriptome data of T. tridentatus midgut (TTG) and fecal mucosa (TTF) were 44.29% and 40.71%, respectively, slightly higher than those of C. rotundicauda, which were 40.08% and 38.81% (Table 1). Clean reads (73.62% – 79.92%) from 18 samples could be mapped to the reference genome, suggesting that the quality of our sequencing data was high. The genome assembly resulted in 373069 unigenes with an N50 of 1314, while the transcriptome assembly resulted in 736378 unigenes and an N50 of 2121. These 373069 unigenes could be assigned to the protein database for annotation information. Next, we analyzed the length distributions of the unigenes and transcripts in the assembly samples (Fig.1). The N50 value of transcript length was longer than 1100 bp, and the N50 value of unigene length was longer than 900 bp.

3.2 Functional Annotation and Classification of the Assembled Unigenes

For the purpose of predicting the functions of the unigenes, all 373069 sequences were analyzed by the database, and the success rates of the annotated data are listed in Table 2. Among the 373069 unigenes, more than 30% showed significant matches in the PFAM database, GO database, NT database, NR database, or Swiss-Prot database. A total of 132594 unigenes (35.54%) showed significant matches in both the PFAM and GO databases.

There were 10482 genes in the digestive tracts of the horseshoe crabs that exhibited significant matches with those of Oncorhynchus tshawytscha. The top 10 species among these matched unigenes were Oncorhynchus tshawytscha, Pygocentrus nattereri, Xenopus laevis, Latimeria chalumnae, Rhincodon typus, Xenopus tropicalis, Callorhinchus milii, Nanorana parkeri, Astyanax mexicanus, and Scleropages formosus. The ratio of the species classification with the highest similarity to the annotated genes is shown in a pie chart (Fig.2). These species were distributed among fishes, amphibians, and reptiles. This result may provide clues for the biological evolution of horseshoe crabs.

The GO annotations of the assembled unigenes from the transcriptomes of T. tridentatus and C. rotundicauda digestive tract are shown in Fig.3. There were 132594 contigs assigned to three GO terms classified as cellular component (CC), biological process (BP), and molecular function (MF). These were assigned to 26, 20, and 10 subcategories, respectively. For biological process, cellular process (76617, 50%), including cellular macromolecule metabolic process (11%) and protein metabolic process (5.5%) were the most represented, followed by metabolic process (70900, 45.3%), single-organism process (62048, 39.7%), biological regulation (28216, 18%), regulation of biological process (26597, 17%), and response to stimulus (18746, 5.2%). Under the molecular function category, binding (64517, 41.3%) including nucleic acid binding (13.1%) was the most abundant, followed by catalytic activity (59773, 38.2%) and transporter activity (10947, 7%). Among the cellular component terms, cell (37831, 22.1%) and cell part (37831, 22.1%) occupied the same dominant subcategory (Fig.3). These annotations represent a profile for gene expression of T. tridentatus and C. rotundicauda, suggesting that there were diverse protein coding genes in the digestive tract closely related to digestion, transportation, and nutrient absorption. These results offer many research resources for further investigation of cellular, metabolic, regulatory, binding, catalysis, and transport mechanisms. In addition to a large number of the annotated genes being assigned to digestion, absorption, and transportation processes, some important genes involved in innate immune systems and nervous systems were also found.

Fig.4 illustrates the KOG annotation of the assembled unigenes from the transcriptome of T. tridentatus and C. rotundicauda. Totally 51910 unigenes were assigned into 26 KOG families. Among the KOG categories, 'translation, ribosomal structure and biogenesis' was the largest category with 7926 unigenes (15.3%), followed by 'posttranslational modification, protein turnover, chaperones' (6539 unigenes, 12.6%), 'general function prediction only' (5579, 10.7%), 'signal transduction mechanisms' (4669 unigenes, 9%), 'energy production and conversion' (4489 unigenes, 8.6%), and 'amino acid transport and metabolism' (4256 unigenes, 8.2%), indicating the specific molecular responses and mechanisms of transportation and transformation of food nutrients.

Among 373069 assembled unigenes, a total of 4768 unigenes were mapped to KEGG pathways (Fig.5) and categorized into 32 subclasses. The most enriched pathways were translation (599 members), environmental information signal processing / transduction (554 members), amino acid metabolism (503 members), carbohydrate metabolism (488 members), endocrine system (393 members), folding, sorting and degradation (285 members), immune system (237 members), nervous system (197 members), nucleotide metabolism (189 members), digestive system (180 members), metabolism of cofactors and vitamins (179 members), and energy metabolism (167 members). This demonstrated that our transcriptome database had various unigenes related to metabolism and environmental information signal processing / transduction. A well-categorized and annotated resource of the transcriptome in digestive tract of T. tridentatus and C. rotundicauda can become a valuable database for investigating their specific bioprocesses and identifying important functional genes related to the physiological processes of digestion and absorption.

3.3 Differentially Expressed Genes Between T. tridentatus and C. rotundicauda

The stomach (TTS, CRS), midgut (TTG, CRG), and fecal mucosa (TTF, CRF) transcriptomes of T. tridentatus and C. rotundicauda were analyzed for the presence of DEGs in comparisons of the two horseshoe crab species. The total numbers of upand down-regulated DEGs in each tissue are listed in Table 3. Venn diagrams depict the overall distribution of DEGs in each of the comparisons. Differential expression analysis revealed that there were 65889 DEGs between the T. tridentatus stomach and C. rotundicauda stomach, while there were 2001 unigenes differentially expressed between T. tridentatus midgut and C. rotundicauda midgut; 1109 unigenes were up-regulated, and 892 were down-regulated. In a total of 7162 DEGs in the comparison between T. tridentatus fecal mucosa and C. rotundicauda fecal mucosa, 7160 unigenes were up-regulated, and only two genes were down-regulated. Venn diagrams (Fig.6) depict the overlap between different sets of DEGs in the three comparisons of three tissues in the two species. Compared with C. rotundicauda, 99.95% of the DEGs in the T. tridentatus stomach were down-regulated, while 99.97% of the DEGs of T. tridentatus fecal mucosa and 55.42% of the DEGs of the T. tridentatus midgut were up-regulated.

3.4 GO Functional Classification of DEGs

After identifying the DEGs, the top 20 terms of the GO enrichment analyses in each basic functional category were selected, and their basic functions combined with GO annotation were described. A good understanding of the biological function of the digestive tract groups between T. tridentatus stomach (TTS) and C. rotundicauda stomach (CRS) could be obtained by the standard GO classification. The DEGs of GO functional classification were assigned into three functional groups: biological process, molecular function, and cellular component. Three categories of GO functional classification comparing the tissue groups (TTS – CRS, TTG – CRG, and TTF – CRF) are shown in Figs.7A, B, C.

Among these GO terms, most DEGs in the TTS-CRS group involved in the biological process category were metabolic process (24283 genes), organic substance metabolic process (19471 genes), cellular metabolic process (18088 genes), primary metabolic process (18505 genes), and nitrogen compound metabolism (14787 genes); most of the DEGs were classified into the cellular component and macromolecular complex (7911 genes), while in the molecular function category the corresponding annotation was in catalytic activity (21013 genes). For the DEGs in the TTG – CRG group, most were assigned to the biological processes category: organic substance metabolic processes (589 genes), protein metabolism (241 genes), and organonitrogen compound biosynthesis (184 genes). In the molecular function category, most of the DEGs were involved in hydrolase activity (245 genes), and most that were assigned to the cellular component were in cell cytoplasm (181 genes). For the DEGs in the TTF – CRF group, most DEGs that were assigned to the biological process category were in organic substance metabolic process (2311 genes), primary metabolic process (2205 genes), and cellular metabolic process (2126 genes). In the molecular function category, most of the DEGs were involved in catalytic activity (2319 genes); the following groups were assigned to the cellular component (1377 genes, GO: 0005623) and cell part (1377 genes, GO: 0044464). The DEGs of the stomach involved in the GO terms were much more frequent than those of the midgut and fecal mucosa, and DEGs had different distribution characteristics in the stomach than in the other two tissues. Most DEGs of the midgut were involved in hydrolase activity, protein metabolism and cell cytoplasm, while most of the DEGs of the stomach and fecal mucosa were involved in catalytic activity and primary metabolic processes.

4 Discussion

Horseshoe crabs are ancient marine arthropods with a long evolutionary history. It is very possible that they have benefited from their special anatomical structure and function, especially the omnivorous feeding habits and effective innate immunity. However, their genetic mechanisms underlying their ability to adapt to environmental changes during their long evolutionary history are still unclear. There is no doubt that as a typical consumer of the coastal and marine ecosystems, the key factors for horseshoe crabs' survival include feeding and the related digestion and absorption.

Here, we report the first transcriptome assembly for the stomach, midgut, and fecal mucosa of T. tridentatus and C. rotundicauda juveniles using the Illumina RNA-seq platform. We obtained from 40.37 million to 52.95 million clean reads and from 6.06 Gb to 7.94 Gb of clean bases from the samples. The genome assembly resulted in 373069 unigenes with an N50 of 1314, while the transcriptome assembly resulted in 736378 unigenes and an N50 of 2121. These 373069 unigenes could be assigned to the protein databases for annotation. Most unigenes were involved in important biological processes related to digestion, absorption, and transportation.

Some annotated genes involved in immune systems and nervous systems were also found in the digestive tract. For example, multiple genes were involved in the immunityrelated JAK-STAT signaling pathway. The functions of the JAK-STAT signaling pathway include defense response to viruses, positive regulation of the immune response, signal transduction, response to drugs, oxidation-reduction reactions, and other important processes. Thus, the digestive tract of the horseshoe crab is a complex system with multiple functions that play important roles in survival and development.

There were a large number of genes in the digestive tract of the horseshoe crab that were similar to those of fishes, amphibians, and reptiles, as identified by annotation in the NR, NT, KO, SwissProt, PFAM, GO, and KOG databases. This may be due to horseshoe crabs having similar living habits, anatomical structure, or function of digestive organs. The structure of the digestive tract of the horseshoe crab is related to its omnivorous habit, and the function of secreting mucus is similar to that in other animals. The structure and function of the digestive tract of the horseshoe crab are closely related to its omnivorous feeding habit.

The previous research on comparative genomics shows that T. tridentatus and the Atlantic horseshoe crab L. polyphemus have the most orthologues shared among horseshoe crab species (Liao et al., 2019). In this comparison between T. tridentatus and C. rotundicauda juveniles, the differences in gene expression in the same tissues between the two species also provide potential molecular mechanisms with which to explore the niche differentiation of two co-occurring horseshoe crab species in southwestern China, and the mechanism of differential expression warrants further investigation.

The DEG assay showed that most DEGs in the stomach transcriptome data of T. tridentatus vs. C. rotundicauda were down-regulated, while the DEGs in the fecal mucosa transcriptome were mostly up-regulated. The different transcriptome patterns between organs might be related to the digestion process during feeding. The digestive tract is divided into foregut, midgut, and hindgut. The stomach being enlarged to store and grind food is the main characteristic of the foregut. The anatomical structure and function are similar between the midgut and hindgut (Chen et al., 2016). They are related to the digestion of organic matter in food, absorption of the hydrophilic and hydrophobic nutrients, the secretion of intestinal fluids, and the transportation of food through the digestive tract, which can help horseshoe crabs adapt to an omnivorous diet. The cell processes and activities in the midgut, stomach, and fecal mucosa can be different during digestion. Further investigation showed that most DEGs of the midgut were involved in hydrolase activity, protein metabolism, and cell cytoplasm, while most DEGs of the stomach and fecal mucosa were involved in catalytic activity and primary metabolic processes, and most DEGs of the stomach were assigned to the cellular component of macromolecular complexes. Most DEGs of the fecal mucosa were assigned to the cellular component and cell parts categories.

The transcriptome data of fecal mucosa of horseshoe crabs may serve as a reference for environmental DNA (eDNA) research on horseshoe crab conservation. The analysis of eDNA has significant potential for monitoring horseshoe crab populations, thus influencing conservation efforts. The eDNA technology is a good biological utility for conservation monitoring without requiring the collection of the living horseshoe crabs. When individuals interact with the environment, DNA is expelled with the excretions and accumulates in the surrounding environment. Sources of eDNA include the horseshoe crab's feces, mucus, and shells. High-throughput DNA sequencing methods including metagenomics, metabarcoding, and single-species detection can be used to analyze these samples for rapid monitoring and application to horseshoe crab conservation. The transcriptome data of fecal mucosa herein shows that it can be a good source of experimental materials for eDNA technology indicating environmental traces of horseshoe crabs.

The present study provides a better understanding of the adaptation of two horseshoe crab species to an omnivorous diet and presents an available transcriptome resource for future research. The horseshoe crabs require substantial food resources, a suitable hydroclimate, and deeper water regions of the bay and the adjacent continental shelf (Tanacredi et al., 2009). While the horseshoe crab's life history with the distinct stages in complex habitats make it vulnerable to various threats at each stage of its development, horseshoe crabs have survived on Earth for more than 450 million years. Therefore, most genes annotated from the digestive tract were assigned to key life processes including metabolism, biogenesis, signal transduction, transportation, immune response, regulatory, binding, catalysis, secretion, and environmental information signal processing. The digestive tract may also play a key role in host adaptation to complex and diverse environments and resistance to the invasion of pathogenic microorganisms in the multiple inshore and estuarine environments. Highly effective adaptation significantly decreases the need for differentiated phenotypic variants associated with environmental impacts, thereby providing the conditions for long-term evolutionary success. The horseshoe crab is an ancient creature that is both simple and sophisticated regarding its ability to adapt to the changes in the global environment.

5 Conclusions

The transcriptome assemblies for the stomach, midgut, and fecal mucosa of T. tridentatus and C. rotundicauda juveniles from the northern Beibu Gulf in March 2019 were reported here for the first time. The genome assembly resulted in 373069 unigenes with an N50 of 1314, while the transcriptome assembly resulted in 736378 unigenes and an N50 of 2121. The unigenes annotated from the digestive tract of horseshoe crabs showed the highest similarity to those of fishes, amphibians, and reptiles, and most unigenes were closely related to the functions of metabolism, binding, translation, biogenesis, catalysis, signal transduction, energy production, transportation, immune response, secretion, and environmental information signal processing. Compared with C. rotundicauda, nearly all of the (99.95%) DEGs in the T. tridentatus stomach were downregulated, while almost all (99.97%) of the DEGs of T. tridentatus fecal mucosa and about half (55.42%) of the DEGs of the T. tridentatus midgut were up-regulated. Most DEGs of the midgut were involved in hydrolase activity, protein metabolism, and cytoplasm, while most DEGs of the stomach and fecal mucosa were involved in catalytic activity and primary metabolic processes, and most DEGs of the stomach were assigned to the cellular component of macromolecular complexes. Most DEGs of the fecal mucosa were assigned to the cellular component of cells and cell parts. The transcriptome data of the fecal mucosa can provide a reference for environmental DNA research for endangered species. The results will also facilitate the investigation of molecular mechanisms of the digestive tract related to feeding habits and environmental traces of horseshoe crabs.

Acknowledgements

This work was supported by the Scientific Research Project of Huaqiao University (No. 605-50X18005), the Guangxi BaGui Youth Scholars Program, and the Guangxi Recruitment Program of 100 Global Expert.