C. japonica and P. trituberculatus were collected from a marine benthic survey at the coast of Yantai City (Yellow Sea, China) in June 2018. C. sapidus was purchased from local seafood market of St. Leonard, Maryland, US. The muscles of crab pincers were sampled and preserved in 95% ethanol for DNA extraction. The genomic DNA were extracted using traditional phenol-chloroform method. The DNA was broken into fragments of 350 bp using an ultrasonicator (Covaris Inc. Woburn, US), then the library was constructed through the steps of end repair, addition of A at the 3’ end, addition of connector, target fragment recovers using electrophoresis, PCR amplification of the target fragments. The libraries were sequenced on the Illumina HiSeq X platform at Gene Denovo Co. (Guangzhou, China). The raw data were generated after base calling from the sequencing with each base a quality score. The ratios of bases with quality score above 20 and 30 called Q20 and Q30 respectively were calculated. The raw data also include the adapter bases and unknown bases (N). The clean reads were obtained by removing those reads containing adaptor, the reads of containing >10% N, and low-quality reads (containing >40% base with Q <= 10).

The K-mer analysis was performed to predict the genome size, heterozygosity and repeat sequences. K-mers are substrings of length k contained within a biological sequence in bioinformatics. In sequence assembling, the term K-mer refers to all of a sequence’s subsequences of length k. The genome size can be estimated using the equation below:

Where the K-mer number is the total number of K-mer in the read set, the K-mer depth is the expected value. The K-mer depth frequency distribution is assumed to follow a Poisson distribution, therefore the expectation of K-mer depth can be estimated from the depth frequency distribution curve. The heterozygosity ratio and repeat sequences were estimated based on the methods described by Liu et al. (2013). All the K-mer analysis were performed using Jellyfish (Marçais and Kingsford, 2011) and GenomeScope (Vurture et al., 2017; http://qb.cshl.edu/genomescope/). In this study, we used k = 21 for K-mer analysis which is most frequently used. The clean reads were assembled into contigs in software SOAPdenovo2 by applying the de Bruijn graph structure (Luo et al., 2012), and then the contigs were spliced into scaffolds. Genome-wide Microsatellite Analyzing Tool Package (GMATA, Wang and Wang, 2016) was used to identify all the microsatellite loci, statistics the location, length, repetitions of the repeat motif, and design the primers for ideal microsatellite makers. The parameters were set as below: min-length of motif is 2, the max length is 6, and the min-repeat times is 5. The maker design parameters are min amplicon size = 120 bp, max amplicon size = 400 bp, optimal annealing temperature = 60°C.

We employed two approaches two measure the genome similarity among the three crab species, 1) the Average Nucleotide Identity (ANI, Ciufo et al., 2018; Jain et al., 2018) and 2) the number and cumulative length of conserved regions. The algorithm of ANI is described by the formula:

The ID%, Length of the Alignment, and Length of the shorter fragment terms refer to the individual blast hits from genome–genome comparisons (Gosselin et al., 2022). The ANI was calculated using FastANI v1.33 (Jain et al., 2018) with the parameter set-up of K-mer size = 16 and Fragment Length = 3,000. The genome conserved regions were identified using MUMmer 4.0 (Delcher et al., 2003; Marçais et al., 2018). The minimum length of a cluster of matches is 65, the minimum length of a single exact match is 20, and the maximum gap between two adjacent matches is 90. The alignment results were further filtered using identify > 85% and minimum aligned length > 1,000 bp and >5,000 bp respectively.

The nucleotide sequence of the identified conserved regions (>5,000 bp) that appear in all the three species were extracted using BEDTools (Quinlan and Hall, 2010). Because a chromosome-level assembled genome of P. trituberculatus has been available (Tang et al., 2020; NCBI Reference No. ASM1759143v1), the gene names of those conserved regions can be acquired through BLAST (Altschul et al., 1997), and then the gene ontology (GO) terms were obtained using gene name from NCBI database (https://www.ncbi.nlm.nih.gov). The GO enrichment analysis was performed using the OmicShare tools, a free online platform for data analysis (https://www.omicshare.com/tools).

The gene that has the longest conserved region in all three crab species was selected for further analysis. The gene characteristics and the location of the conserved region on the gene were identified by BLAST with P. trituberculatus genome (Tang et al., 2020) shown in NCBI Sequence Viewer (Rangwala et al., 2021). The protein domains were predicted using SMART (a Simple Modular Architecture Research Tool; Letunic et al., 2021; http://smart.embl-heidelberg.de). The homeotic gene fragments from other known-genome Decapoda species were obtained by BLAST the coding sequence of the similar region. A phylogenetic tree was constructed using MEGA 11 (Tamura et al., 2021) based on Neighbor-joining method (Saitou and Nei, 1987).

The summary for the sequencing data is shown in Table 1. The minimum Q20 and Q30 are 95.09% and 89.02% which meet the sequencing accuracy requirement of illumine platform that the Q20 and Q30 should be at least 90% and 85%. The estimated genome sizes based on the 21-mer analysis (Figure S1, Table 2) are 737.6 Mb, 994.4 Mb and 993.4 Mb for P. trituberculatus, C. japonica, and C. sapidus respectively. The observed heterozygous rate is larger than 1% for all three species. The repeat sequence content ranges from 58.4% to 62.6% (Table 2). The preliminary genome assembly results are shown in Table 3. Total of 485,821 scaffolds with the total length of 762 Mb were obtained for P. trituberculatus, 655,217 scaffolds with the total length of 1Gb were obtained for C. japonica, and 706,990 scaffolds with the total length of 858 Mb were obtained for C. sapidus. The N50 length is short for all three species, ranging from 1,397 (C. sapidus) to 2,131 (P. trituberculatus). The GC content for the assembled genome sequences ranges from 39% - 42% for the three species, which fall into the GC preference range of Illumina sequencing (25% - 65%). Total of 1.3 million, 1.0 million and 1.7 million microsatellite loci were identified from the assembled genome sequences of P. trituberculatus, C. japonica and C. sapidus (Table S1). Of the five motif repeat types, the dinucleotide repeat occurs most frequently in all the three species. The top grouped dinucleotide motifs are GT/AC in both P. trituberculatus and C. japonica, while it is different in C. sapidus (TC/GA, Figure S2). The microsatellite markers and their primers developed from these repeat motifs are listed in Table S2.

The ANI between P. trituberculatus and C. japonica (82.3%) is close with that between P. trituberculatus and C. sapidus (82.6%), while the ANI between C. japonica and C. sapidus is relative smaller (80.3%, Table 4). The MUMmer results revealed 16,414 conserved fragments (alignment length > 1,000 bp, identity > 85%) with total length of 26.5 Mb between P. trituberculatus and C. japonica, 12,745 conserved fragments with total length of 20 Mb between P. trituberculatus and C. sapidus, and 9,366 conserved fragments with total length of 14.7 Mb between C. japonica and C. sapidus (Table S3). The results were further filtered using alignment length > 5,000 bp, the numbers of conserved fragments are 126, 67, and 58 respectively (Table 4). Among these similar fragments, total of 40 were found in all three crab species. There are 34 fragments from 28 genes that have been annotated. The gene names are shown in Table 5. These genes are mostly involved in 22 biological processes, 6 molecular functions and 13 cellular components (Figure 2A). The top 20 significantly enriched ontologies are shown in Figure 2B, including regulation of cardiac conduction, circulatory system process, heart contraction, blood circulation, muscle structure development, synaptic membrane, etc.

The longest conserved region among the three species was found in the gene of ankyrin-2 (ANK2, Table 5). The whole length of P. trituberculatus ANK2 (PtANK2) is 391,211 bp, which has 34 isoforms with total of 25 - 43 exons (Tang et al., 2020). The identified conserved region in ANK2 consists of two separate fragments. The longer fragment (CF1) locates from 29,981,076 to 29,994,164 on Chromosome 46, which mostly locates in the intron area and contain 1 – 2 short exons. The shorter fragment (CF2) is a sequence of 10,602 bp that contains the longest exon of PtANK2. That exon appears in all isoforms (Figure 3). The PtANK2 protein contains a string of ankyrin repeats, a ZU5 domain and a DEATH domain in the first 1/3 region, the rest region is characterized by low complexity proteins where both CF1 and CF2 locate (Figure 3).

The protein sequences translated from the long exon of CR2 were found highly conserved in other known-genome Decapoda species. Protein sequences from other 9 species including the Chinese mitten crab Eriocheir sinensis (Accession No.: PRJNA636904), the snow crab Chionoecetes opilio (PRJNA602365), the blue king crab Paralithodes platypus (PRJNA555178), the coconut crab Birgus latro (PRJNA704570), the American lobster Homarus americanus (PRJNA655509), the red swamp crayfish Procambarus clarkii (PRJNA727411), the Kuruma shrimp Penaeus japonicus (PRJDB11151), the whiteleg shrimp Penaeus vannamei (PRJNA438564), and Chinese white shrimp Penaeus chinensis (PRJNA691453) were extracted from their genome database in NCBI. The phylogenetic tree of the 12 species based on neighbor-joining method coincides with the current taxonomy and evolutionary relationship of Decapods (Figure 4). The P. trituberculatus and C. sapidus clustered first then connected to C. japonica. The three swimming crabs then grouped with the other two Brachyura crabs C. opilio and E. sinensis successively. B.latro and P. platypus forms a clade that belongs to Anomura. H. americanus and P. clarkia make up a cluster of Astacidea that consists of lobsters and crayfish. All Brachyura, Astacidea are infraorders of Pleocyemata suborder. The three prawn species belong to the other suborder Dendrobranchiata of Decapoda. The genetic distance reflected by tree branch length between Astacidea and prawns is smaller than that between either crab infraorder and prawns, indicating the lobsters and crayfish first diverged from prawn, this is consistent with “Carcinization” process (McLaughlin and Lemaitre, 1997).

The heterozygosity is one of important characteristics of diploid or polyploid organism’s genome, which is caused by the genetic variations between different chromosomes. High heterozygosity has been a challenge for the estimation of genome size and genome assembly (Kajitani et al., 2014; Gao et al., 2018). For diploid or polyploid organisms, only the sequences from one set of chromosomes should be counted and assembled. However, the sequences from high heterozygous regions will be also included, therefore may double the region. On the K-mer depth frequency distribution curve, it often manifests as multiple peaks or tailing, which will cause deviation when estimating the expectation of the K-mer depth. The heterozygosity is usually used to divide genomes into low heterozygosity (0.5% ≤ heterozygous ratio <0.8%) and high heterozygosity (heterozygous ratio ≥0.8%) (Li et al., 2019). In this study, high heterozygosity was detected in all three crab species (>1.2%). The genome heterozygosity is generally high (ranging from 1% to 4%) for most crabs and bivalves according to the existing records (Cai et al., 2019; Bachvaroff et al., 2021; Penaloza et al., 2021; Tang et al., 2021), compared to fish genome with usually lower than 0.5% (Chen et al., 2020; Xu et al., 2020; Choi et al., 2021).

Organisms’ genome is also characterized of a certain fraction of repetitive sequences that may originated from DNA duplication and recombination during the evolution. The content of the repeat sequences increases from virus (<1%) to maize (77%; Gao et al., 2018). The human genome of about 66%–69% is repetitive (de Koning et al., 2011). When the repeat content is larger than 50%, it is recognized as high repetitive genome (Tang et al., 2015). In this study, the repeat sequences of all the three crab species are above 50%, the C. japonica and C. sapidus are larger than 60%. The assembly becomes very difficult for high repetitive genome when only using the next-generation short-read sequencing. The high heterozygosity and high repeat content of all the three crab species resulted in that the assembled scaffolds are very short. For high-quality assembly of a high heterozygosity and high repeat content genome, the combination method of second-generation short-read sequencing (Illumina) and third-generation long-read sequencing (PacBio) supplemented with Hi-C technology are usually used (Li et al., 2019; Tang et al., 2020; Bachvaroff et al., 2021).

Microsatellites is one kind of tandem repeat sequences typically with motif of 1–6 nucleotides (Richard et al., 2008). It is widely distributed on the whole genome, characterized by high polymorphism, co-dominant inheritance and easy to capture (Bhattarai et al., 2021). The microsatellites are very dense in the genomes of all the three species (~1,320/Mb). The microsatellite loci developed from this study provide useful tools for population genetics and selective breeding, such as fingerprinting, linkage mapping, parentage identification, or marker assisted selection etc.

In taxonomy, both P. trituberculatus and C. sapidus are from Portuninae subfamily and they have more similar morphological appearance, compared to C. japonica that is from a different subfamily - Thalamitinae. The previous phylogenetic evolutionary trees constructed using mitochondrial DNA coding sequences also showed that P. trituberculatus and C. sapidus clustered together with small genetic distance and were in the Portuninae branch, which had an explicit separation with the Thalamitinae species that C. japonica is included (Mantelatto et al., 2018; Evans, 2018). Moreover, our phylogenetic tree constructed using a highly conserved ANK2 protein fragment is consistent with the phylogenetic and taxonomy relationship in which P. trituberculatus is closer with C. sapidus rather than C. japonica. However, for the comparison of genome similarity among three crab species based on the genome draft of this study, a high ANI was observed between P. trituberculatus and C. japonica that are from different subfamilies. It is very close to the ANI between P. trituberculatus and C. sapidus, the species within subfamily. In contrast, the C. sapidus genome is less similar with C. japonica that may reflect a normal subfamily-level differentiation. In addition, both the number and the cumulative length of the identified conserved regions (>1,000 bp and >5,000 bp) between P. trituberculatus and C. japonica are larger than those between P. trituberculatus and C. sapidus, which may indicate the genome of P. trituberculatus is more similar with C. japonica compared to C. sapidus.

We speculate that the high genome similarity between P. trituberculatus and C. japonica might be attributed to their similar living environments. Both P. trituberculatus and C. japonica occupy the same areas across the western coast of Pacific, while C. sapidus lives in the west coasts of Atlantic which is totally separate in geography from P. trituberculatus’ inhabit range, where the hydrological conditions, climate and biological communities are different. The environment cannot change the genome directly but may play as a selection pressure to impact the genome evolution. Based on the theories of natural selection (Darwin, 1859) and modern synthesis (Fisher, 1958; Wright, 1932; Dobzhansky, 1937), mutations randomly occurred in the genome of organisms, the genome evolved by retaining the mutations that adapt to selection pressures and wiping out the ones that do not adapt.

However, the inconsistency between entire genome similarities and phylogenetic evolutionary relationship inferred from previous studies might be caused by the difference of rate of evolution between coding DNA (cDNA) and non-coding DNA (ncDNA). The organism’s genome consists of coding cDNA and ncDNA sequences. The cDNA encodes proteins that are the basic component for the cell structure, function, and regulation of the body’s tissues and organs, but only comprises 1% of the entire genome (Rands et al., 2014), whereas the other 99% does not encode any proteins which had been ever thought as useless “junk” (Ohno, 1972). The cDNA has much slower mutation rate than ncDNA (Kimura, 1983) because they are so related to the life activities of the organisms and therefore the mutations are hardly retained. Compared to cDNA, the super high occupation ratio of ncDNA on the genome means much higher mutation occurring frequency and more mutations. Because they do not encode protein, their mutations would have higher possibility to be retained. Despite the mutations are considered occurring neutrally based on the neutral theory (Kimura, 1983), they are subject to selection as long as the genome region where they are located are functional. It is becoming clear that the ncDNA plays essential roles in the regulation of gene expression (Ludwig, 2016). A comparative genome study between D. melanogaster and its related sibling species D. simulans demonstrated that a large fraction of ncDNA regions is functionally important and subject to both purifying selection and adaptive evolution. Adaptive changes to ncDNA might have been considerably more common in the evolution of D. melanogaster (Andolfatto, 2005). High mutation occurrence and subjecting to selection for ncDNA make it possible that the similarity in the entire-genome is higher than some conserved cDNA sequences between related species that live in the same environment.

The identified longest conserved region in all three crab species is contained in the arkyrin-2 gene (ANK2). The ANK2 encodes a member of the ankyrin protein family that are adapters facilitating the local organization of integral membrane proteins with cytoskeletal elements (Cunha and Mohler, 2008; National Center for Biotechnology Information (NCBI), 2022). The ankyrin protein also play important roles in activities such as cell motility, activation, proliferation, contact and the maintenance of specialized membrane domains (National Center for Biotechnology Information (NCBI), 2022). This protein is active in many cell types, particularly in the brain and in cardiac muscle. It interacts with ion channels and ion transporters that generates the electrical signals controlling the heartbeat and maintain a normal heart rhythm (Mohler and Bennett, 2005; Cunha and Mohler, 2008). Such function of ANK2 coincides with the GO enrichment result those the most significant gene ontologies are cardiac-related. Most ankyrins are typically composed of three structural domains: a consecutive ANK repeats (EMBL-EBI ID: IPR002110) that associate with ion channels, transporters, and cell adhesion molecules, such as the Na+/Ca2+ exchanger and Na/K-ATPase; a central region with a highly conserved spectrin binding domain (ZU5, EMBL-EBI ID: IPR000906) that participates in induction of apoptosis and binding of melanoma-associated antigen; and an Alpha-helical domain (Death, EMBL-EBI ID: IPR000488) that is related in sequence and structure to the death effector domain and the caspase recruitment domain. None of the similar regions that we identified spans the domain area, this is because the proteins of this area are encoded by very short exons, the genetic variations in introns are more than those in exons, the longer exons mainly distribute in the second half area, thus the similar regions (>1,000 bp) are mostly detected in this area. The mRNA of ANK2 contain multiple isoforms. The number of isoforms is very high in some known Decapoda genomes (P. trituberculatus: 38; H. americanus:32; P. japonicus: 34). These isoforms are generated by alternative splicing. Alternative splicing is the primary mechanism through which the genome generates mRNA and protein diversity from a given coding repertoire. Alternative splicing plays important roles in regulating development, physiology and homeostatis, providing an alternative program in response to environmental stress (Mastrangelo et al., 2012; Kelemen et al., 2013; Liu and Guo, 2017). The P. trituberculatus ANK2 contains 25 ~ 43 exons, almost all of which are missing in some isoforms except the longest exon close to the 3’ end. That exon is included in one of the two longest conserved regions (CR2, Figure 3).

The protein sequence translated by this exon is found very conserved in other Decapoda species. The phylogenetic analysis using this protein fragment from 12 Decapoda species revealed a clear taxonomy and evolution relationship. The specification of crabs referring to a “Carcinization” or “Brachyurization” (McLaughlin and Lemaitre, 1997). Under that hypothesized scenario, the prawns (Dendrobranchiata) first diverged from Decapod ancestors, then a group evolved from swimming to crawling to walking (Reptantia, the ancestors of lobster/crayfish and crabs), during this process, their abdomen and tail gradually become short and flat, and start folding toward below the carapace, finally a clade of this group formed the nowadays’ crab shape. In the phylogenetic tree, the genetic relationship between lobster/crayfish and crabs is closer than that between lobster/crayfish and prawns or between crabs and prawns, indicating later divergence of these two clades, however, the genetic distance between lobster/crayfish and prawns is smaller than that between crabs and prawns indicating the lobster/crayfish first diverged from prawns and their appearance is also more like the spawns than crabs. The evolution of this ANK2 gene may occur accompanying the evolution process of the species.