The healthy juvenile Chinese mitten crabs were taken from Tianjin Xieyuan Aquaculture Co., Ltd. The average body length, body width, and weight of crabs were 24.28 ± 1.54 mm, 26.97 ± 1.65 mm and 8.41 ± 1.47 g, respectively. The crabs were cultured in a plastic incubator (70 × 40 × 50 cm) at 20 ± 1°C with natural light for more than 7 days. The crabs were fed compound feed twice a day. The daily feeding amount was about 5% of the body weight. The siphon method was used to clean up the residual bait and inject new water. The daily water change is more than 50%. Three crabs were selected to collect the muscle samples from pre- and post-molt stages, respectively. The muscle tissues in the walking legs of crabs in pre- and post-molt stages were collected. The post-molt stage was determined by observing the crab’s exoskeleton. In the post-molt stage, the crab’s body is still soft, and the exoskeleton has not been completely calcified, which are obvious features that can be determined by direct observation. The pre-molt stage was determined with microscope observation. The third maxilliped was cut carefully with anatomical scissors and placed downward on a glass slide with clean water and then observed with an Olympus (OCKX41, Japan) optical microscope (magnification of 200–400×). The most obvious feature of the pre-molt stage is the separation of epithelium and exoskeleton, forming a transparent blank area, which marks the start of animal molting (Kang et al., 2012). Therefore, when the blank area between the epithelium and the exoskeleton is observed under the microscope, it is considered that the crab is in the pre-molt stage ( Supplementary Figure S1 ). The crabs in pre- and post-molt stages were put on the ice plate for low-temperature anesthesia. The muscle tissue in the walking legs of the crabs was quickly collected and put in liquid nitrogen immediately and then preserved in a −80°C refrigerator. The muscle in one crab was separated into two samples for ATAC-seq and RNA-seq, respectively. Cell activity of the muscle sample was detected with Trypan blue assay and counted.

Six muscle samples were used for ATAC-seq, numbered as Mu_A_1, Mu_A_2, and Mu_A_3 (samples for post-molt stages) and Mu_D_1, Mu_D_2, and Mu_D_3 (samples for pre-molt stages). ATAC-seq was performed as previously reported (Corces et al., 2017). Briefly, nuclei were extracted from muscle samples, and the nuclei pellet was resuspended in the Tn5 transposase reaction mix. The transposition reaction was incubated at 37°C for 30 min. Equimolar adapter 1 and 2 were added after transposition. PCR was then performed to amplify the library. After the PCR reaction, libraries were purified with the AMPure beads, and library quality was assessed with Qubit. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina NovaSeq platform, and 150 bp paired-end reads were generated. Nextera adaptor sequences were first trimmed from the reads using a skewer (0.2.2). Q20 and Q30 of the raw data and clean data were calculated. These reads were aligned to the reference genome downloaded from the NCBI BioProject database (accession number PRJNA555707) using BWA, with standard parameters. These reads were then filtered for high-quality (MAPQ ≥ 13), non-mitochondrial chromosome and properly paired reads (longer than 18 nt). All peak calling was performed with MACS2 using “MACS2 call-peak – nomodel – keepdup all – call-summits.” For simulations of peaks called per input read, aligned and de-duplicated BAM files were used without any additional filtering.

Differential analysis of the two molting stages was performed using the DESeq2 R package (1.20.0). The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Peaks with p-value ≤ 0.05 and |log2.FoldEnrich| ≥1 were assigned as differential peaks.

ChIPseeker software was used to stat the distribution of peak in each functional area. The correspondence between peak and each functional area is in the priority order of promoter, untranslated region (UTR), UTR, exon, intron, downstream TTS, and distal intergenic regions. That is, if the position of a peak may be the exon of one gene and the intron of another gene, according to the above priority, it will be considered as the exon rather than the intron. Goseq and KOBAS software were used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation for the genes associated with differential peaks.

Total RNA from E. sinensis muscle samples was sequenced using the Illumina high-throughput sequencing technology by Novogene Inc. The total RNA was extracted using the TRIzol method (Invitrogen) for transcriptome analysis. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of a Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total amount of 1 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using a NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, USA) following the manufacturer’s recommendations. Library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s protocols. After cluster generation, the library preparations were sequenced on an Illumina NovaSeq platform, and 150 bp paired-end reads were generated.

The original image data were transferred into sequencing data by base calling, which were defined as raw reads. Clean data (clean reads) were obtained by removing reads containing an adapter, reads containing poly-N, and low-quality reads from the raw data. At the same time, Q20, Q30, and GC contents of the clean data were calculated. All the downstream analysis was based on the clean data with high quality. Then, E. sinensis genome was downloaded from the NCBI BioProject database (accession number PRJNA555707) and used as the reference genome. Index of the reference genome was built using Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5. FeatureCounts v1.5.0-p3 was used to count the number of reads mapped to each gene. Then, FPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene. FPKM is the expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced.

Differential expression analysis of two molting stages was performed using the DESeq2 R package (1.20.0). The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with p_adj ≤ 0.05 and |log2.Fold_change| ≥ 1 were assigned as differentially expressed. Goseq and KOBAS software were used for GO and KEGG annotation for the DEGs. The ATAC-seq and RNA-seq data have been deposited to Gene Expression Omnibus (GEO) database with the accession number GSE206233.

The expression profile of the differential peak-related genes in ATAC-seq results was compared with the DEGs in RNA-seq results. The upregulated DEGs in RNA-seq were compared with the associated genes with the upregulated peaks in ATAC-seq. The downregulated DEGs in RNA-seq were compared with the associated genes with the downregulated peaks in ATAC-seq. The genes with the same expression profile in the two methods were taken for GO enrichment analysis.

TMHMM (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0) (Krogh et al., 2001) was used for the prediction of the TMHs for the candidate GPCR genes. PreidictProtein (https://predictprotein.org/) (Bernhofer et al., 2021) and SWISS-MODEL (https://swissmodel.expasy.org/) (Waterhouse et al., 2018) were used to predict the secondary and three-dimensional structure of the proteins encoded by the candidate genes.

The expressions of candidate GPCRs were verified with real-time quantitative PCR (qRT-PCR). Three muscle samples were used in the qPR-PCR validation for each gene as the biological duplication, with each sample containing the mixed muscle tissues from three crabs to avoid individual differences. The primer sequences used in qRT-PCR are listed in Supplementary File S1 . The RNA was extracted with TRIzol reagent (Invitrogen, San Diego, CA, USA). The reverse transcription of RNA was conducted with PrimeScriptTM RT Master Mix (Perfect Real Time) test kit (Takara, Japan). Fluorescence quantitative detection TB Green^® Premix Ex TaqTM II(Tli RNaseH Plus)*2 test kit (Takara, Japan) was used for fluorescence quantitative detection. ABI Applied Biosystems was used in this experiment. The housekeeping gene is β-actin. The 20-μl reaction system contains 1 μl cDNA, 0.4 μl each primer (10 μM), 10 μl TB Green^® Premix Ex TaqTM II (2×), 0.4 μl ROX Reference Dye (50×), and ddH₂O up to 20 μl. The reaction program was set as 95°C for 30 s, followed by 40 cycles of 5 s at 95°C, 30 s at 60°C, and 60°C for 1 min. The data were quantitatively analyzed with 2^–ΔΔCt method.

ATAC-seq was used to detect the landscape of genomic chromatin accessibility in the muscle during molting. After the raw data have been obtained from ATAC-seq, the data were trimmed to get the clean data for the following analysis. The result of quality control for raw and clean data is shown in Table 1 . The clean data were then mapped to the genome, and the results are shown in Table 2 . The mappability is above 80%, which is a quite high ratio. The unique mapped reads is 63% on average for all the samples. The results of mapped read distributions across the gene bodies and transcriptional start site (TSS) showed the good quality of the ATAC-seq ( Figures 1A, B ).

To further confirm the quality of ATAC-seq, Pearson and Spearman correlation analysis was performed based on the signals of merged peaks from all samples. The closer the correlation coefficient is to 1, the higher the similarity of expression patterns between samples. As shown in Figures 1C, D , the similarity of the three samples in each group (D and A group) was quite high. After the mapped reads were obtained, the duplicate reads were deleted, and the peaks were called. The results of peak calling are shown in Table 2.

The “FoldEnrich” values in the two groups (D and A groups) were compared, and 7,163 differential peaks were found, in which 3,404 were upregulated and 3,759 were downregulated. The functional regions of each peak on the genome were annotated. The genome-wide functional regions were divided into promoter, UTR, exon, intron, downstream gene start site (TTS), and distal intergenic regions. The TSS are usually enriched in the promoter region. The functional distribution of differential peaks is shown in Figure 2 . The distal intergenic and intron regions took the largest part in the differential peaks. The percentages of peaks enriched near the TSS region were 14% and 15% in all the up- and downregulated peaks.

The results of GO enrichment analysis for the differential peak-related genes showed that in the biological process, these genes were mainly concentrated in the regulation of metabolic process and gene expression; in the cellular component, the genes were mainly located in the membrane and extracellular; and in the molecular function, the genes were enriched in the catalytic activity, phosphotransferase activity, and protein kinase activity. In the KEGG enrichment analysis, a large number of genes were distributed in the endocytosis, protein processing in endoplasmic reticulum, purine metabolism, and many signaling pathways ( Figure 3 ). The enriched functions of these differential peak-related genes indicate that plenty of specific signaling regulations occur in the muscle during molting, and the regulatory process mainly takes place in the membrane and extracellular. In addition, the large amount of genes with catalytic activity reflects that the metabolic process also undergoes great changes.

Total RNA from muscles in pre- and post-molt stages was extracted to obtain the E. sinensis transcriptome data. The RNA samples had no genomic contamination, protein and impurity contamination, and color abnormality. The final quality testing result was level A, which was the highest level for eukaryotic transcriptome RNA. Total RNA was used as the input material for the RNA sample preparations and then sequenced. The results of RNA-seq are shown in Table 3 . The average percentages at Q20 and Q30 of all the samples are 95.97% and 90.68%, respectively. The average GC content is 55.08%. More than 80% of the clean reads of each sample were mapped to E. sinensis genome. The distribution of the mapped reads on the genome is shown in Figure 4 . The proportion of reads aligned to the exon region was the highest in all the samples. The reads aligned to the intron region may come from the precursor mRNA or the intron retained by variable splicing events. The reads aligned to the intergenic region may be influenced by ncRNA or a few DNA fragments, or the imperfect genome annotation.

In the differential analysis, a total of 1,684 genes were found to be differentially expressed in the A vs. D group, in which 1,388 were upregulated and 296 were downregulated. Through the GO enrichment analysis of these DEGs, it was found that in the biological process, the DEGs were mainly concentrated in several metabolic processes, biological regulation, and gene expression; in the cellular components, the DEGs were mainly located in the membrane; and in the molecular function, the DEGs were mostly enriched in binding, catalytic activity, protein binding, and structural molecule activity. In the KEGG pathway enrichment analysis, lots of genes were enriched in ribosome, amino sugar and nucleotide sugar metabolism, carbon metabolism, biosynthesis of amino acids, purine metabolic, and mitogen-activated protein kinase (MAPK) signaling pathway ( Figure 5). The enriched functions of the DEGs reflect that the metabolic process changes a lot during molting. In addition, some specific signaling regulations such as MAPK signaling pathway were quite different in pre- and post-molt stages.

The associate genes related to differential peaks from ATAC-seq were compared with the DEGs from RNA-seq. Finally, 784 common genes were found, in which 616 had the same expression profile ( Table 4 ). The other 168 genes were common but inconsistent in expression profile, in which 22 genes were upregulated in ATAC-seq but downregulated in RNA-seq, and 78 genes were downregulated in ATAC-seq but upregulated in RNA-seq. Compared to D stage, most of the DEGs were upregulated in A stage. The GO enrichment analysis of the consistent DEGs in biological process, molecular function, and cellular components is shown in Figure 6 . The metabolic process contains the most DEGs, most of which are concentrated in organic substance, nitrogen compound, and macromolecule metabolic processes. The regulation of cellular process contains 44 differential genes, including regulation of metabolic process, gene expression, and transcription. In terms of the cellular component, the most DEGs were distributed in the cell membrane. For the molecular function, the number of DEGs with protein binding, ion binding, and catalytic activities is the largest. All these indicate that the expression levels of many genes involved in metabolism, information transmission and transcriptional regulation have changed during molting.

The molting process of E. sinensis is a physiological process completed by the coordination of multiple organs. Different organs exchange information through signal transmission at each molting stage, so as to accurately control the activities and launch time of multiple organs participating in the molting process. In order to explore the molecular regulation mechanism on the change of muscles during molting, it is important to investigate its signaling transduction process. Therefore, among the DEGs between D and A stages, we focused on the genes related to the signaling transduction. According to GO enrichment analysis, in the consistent DEGs, 19 genes were annotated to the “signal transduction (GO:0007165),” of which 13 are located in the G-protein-coupled receptor signaling pathway. Further analysis on the molecular function of these genes showed that they all had transmembrane signaling receptor activity, of which 10 genes were located on the membrane. These genes were all upregulated, which suggests that some specific GPCRs may play an important role in the process of receiving and transmitting signals in the post-molt stage in muscle.

GPCRs mediate the majority of cellular responses to external stimuli, including light, odors, hormones, and growth factors. GPCRs are integral membrane proteins that contain seven transmembrane α-helices (TMHs) (Weis and Kobilka, 2014). In order to further analyze whether the above 10 membrane-located genes encode GPCRs, TMHMM, PredictProtein, and SWISS-MODEL were used to predict their transmembrane structure, secondary structure, and three-dimensional structure, respectively. The prediction results and FPKM of these genes are shown in Table 5 ( Supplementary Figures S2 – 4 ). According to the analysis of transmembrane structure predicted by TMHMM, five of these genes encode proteins with seven or more TMHs. Through the analysis of secondary structure and three-dimensional structure, U165 and U196 have seven TMHs, which is consistent with the prediction results of TMHMM. U116 is predicted to have seven TMHs according to TMHMM and PredictProtein, while it is predicted by SWISS-MODEL to have two adjacent seven-TMH structures, which is consistent with the characteristics of family C GPCRs. Family C GPCRs are unusual in that, besides the GPCR-defining seven-THM domain, they possess relatively large amino-terminal extracellular domains that form obligate dimers (Koehl et al., 2019). The protein encoded by U120 is predicted to have eight TMHs by TMHMM, but in the prediction of the PredictProtein and SWISS-MODEL, only seven of the α-helices are THMs; the other one does not have the transmembrane characteristics. The protein encoded by U258 is predicted to have eight THMs according to TMHMM. While similar as U120, in the predictions of PredictProtein and SWISS-MODEL, only seven α-helices are TMHs. In addition, although U33 is predicted to have six THMs by TMHMM and SWISS-MODEL, it should have seven THMs according to PredictProtein. It is worth noting that the similarity between the reference template sequence and U33 is <25% in SWISS-MODEL prediction. Therefore, the three-dimensional structure prediction results of U33 is of little reference value. There are a total of six genes predicted to have seven THMs by at least one sort of software. Therefore, we considered these genes as the candidate GPCRs for further experimental validation. The other four proteins (encoded by U519, U241, U9, and U387) are excluded as GPCRs, as they are impossible to have seven THMs in the prediction results of all the three sort of software.

The above six candidate GPCRs were verified by qRT-PCR. As Figure 7 shows, all the seven genes were expressed in both A and D stages. The relative expression levels of all these genes in the A stage were higher than those in the D stage, showing an upregulated trend in A vs. D group, which was consistent with the results from ATAC-seq and RNA-seq. The ratio of standard deviation to mean value is in the range of 12%–26%, which is in a reasonable scale. qRT-PCR results verified the results of high-throughput sequencing. Therefore, these proteins can further verify their role as GPCRs in muscle discontinuous growth during molting by gene knockout in subsequent experiments.

Molting is a critical developmental process for crustacean, which is a complex process regulated by the combination of multiorgans. The molting frequency directly influences the growth rate of E. sinensis. Currently, many organs have been found to participate in the molting process, such as hepatopancreas (Huang et al., 2015), eyestalk (Li et al., 2016), muscle (Wang et al., 2003; Tian et al., 2017), Y organ (Stewart et al., 2013), etc. Some genes were also reported to be involved in the molting process (Zhang et al., 2011; Li et al., 2015). However, the regulation mechanism of the molting process is still unclear. Muscle is one of the organs with obvious changes in the molting process. The study of the mechanism on the initiation and termination of molting in the muscle is important to understand and control the molting cycle of E. sinensis. Presently, the research on muscle changes mostly focuses on the observation and analysis at the tissue level. There is little research at the molecular level, especially for the research on the signal transmission of muscle in the molting process.

In this work, the chromatin accessibility and gene expression of muscle during pre- and post-molt stages were investigated with ATAC-seq and RNA-seq, respectively. ATAC-seq is used for the detection of chromatin accessibility, which reflects the transcriptional activity of chromatin. It is an important method to study the regulation of gene expression and has been widely used in human, plant, animals, and yeast (Hendrickson et al., 2018; Liu et al., 2019; Shashikant and Ettensohn, 2019; Yang et al., 2020). The precise identification of specific open DNA regions in the genome is not only capable of exploring the regulatory elements in the genome but also an effective way to analyze the genes expression by identifying the peak-related genes. The integration of ATAC-seq and RNA-seq is capable of providing more accurate expression profile of genes and proposing more instructive suggestions based on the high-throughput data analysis. It also supplies a better solution for the in-depth analysis of the mechanism on cellular activities in the process of biological growth and development (Lowe et al., 2019). There have been some related studies on the integration of ATAC-seq and RNA-seq to analyze the gene expression of various organisms. For example, Yang et al. established the differential expression profile of light-induced primordia formation in Sparassis latifolia by integrating the assay for ATAC-seq and RNA-seq (Yang et al., 2020). Ackermann et al. determined the genetic expression of human alpha- and beta-cells based on chromatin accessibility and transcript levels, which allowed for the detection of novel alpha- and beta-cell signature genes not previously known to be expressed in islets (Ackermann et al., 2016). In this work, ATAC-seq and RNA-seq techniques were used to analyze the muscle samples in the pre- and post-molt stages of molting circle and determine the specific genes functioning in the muscle before and after molting.

The results of enrichment analysis based on ATAC-seq and RNA-seq reflect that the metabolism and regulation in muscle both undergo great changes during molting. Some enriched GO annotations and KEGG pathways obtained by these two methods are consistent. Lots of the DEGs were enriched in biological regulation process. Most DEGs were located in the membrane, and many of them have catalytic and binding activities. For the pathway distribution, the DEGs were commonly enriched in endocytosis, protein processing in endoplasmic reticulum, carbon metabolism, purine metabolism, and MAPK signaling pathway. The enrichment of DEGs in the biological regulation process and signaling pathway indicates that the regulation mechanisms in pre- and post-molt stages are quite different. The initiation and termination of molting are regulated by different regulatory elements.

The integration analysis of ATAC-seq and RNA-seq showed 538 upregulated genes and 78 downregulated genes, which indicates that more specific genes were active in post-molt stage. The result of GO analysis showed that these DEGs were distributed in a variety of biological process, especially in regulation of metabolic process, gene expression, and transcription. Most DEGs were located in the membrane with protein binding, ion binding, or catalytic activity. In the results of ATAC-seq and RNA-seq, the common DEGs with the same expression profile account for 78.57% (616/784) of all the common DEGs identified in the two methods. This indicates that ATAC-seq and RNA-seq may have a common problem of high-throughput sequencing, that is, false positive results. It further demonstrates that the integration analysis of the two high-throughput sequencing techniques can effectively reduce false positives and obtain more accurate results.

In order to explore the signal transduction mechanism of muscles in the molting process, the DEGs in the signal transduction process were further analyzed. We found that in the 19 DEGs annotated to the signal transduction process, 13 are located in the G-protein-coupled receptor signaling pathway and all upregulated. This result shows that GPCRs play a leading role in muscle signal transmission during post-molt stage. The other six genes are zonadhesin, ankyrin-3, probable nuclear hormone receptor, Toll-like receptor Tollo, hormone receptor, and mitochondrial uncoupling protein 4. Ankyrin-3 and Toll-like receptor Tollo are upregulated in D stage, which indicates that they are specifically active in pre-molt stage, while the other four genes are upregulated in A stage like GPCRs, indicating that they are specifically functioning in post-molt stage.

Currently, there are few studies on the existence of GPCR in E. sinensis. Only Xu et al. cloned a gene belonging to the GPCR family named EsGPCR89 in E. sinensis and found that it played roles in cerebral antimicrobial function via hemocytes (Qin et al., 2019). In terms of the signal transduction, Tian et al. suggested that a transforming growth factor-β type I receptor, named EsTGFBRI, probably play roles in the molting-related muscle growth in E. sinensis (Wei et al., 2015). Wang et al. cloned an activin receptor IIB (ActRIIB) gene in E. sinensis and showed that ActRIIB signaling pathway might be involved in the regulation of growth mass and protein and lipid metabolism in muscle during molting (Wang et al., 2020b). Some studies have found that ecdysone signal transduction in Drosophila is mediated by GPCR (Srivastava et al., 2005), but there is no such research in crustaceans. In this work, it is found for the first time that GPCRs are involved in the molting process of E. sinensis and play an important role in muscle signal transduction.

As a result, six candidate GPCR genes (U116, U165, U196, U120, U258, and U33) were identified by structural analysis and qRT-PCR verification after the integration analysis of ATAC-seq and RNA-seq. U116 encodes metabotropic glutamate receptor 7, which is a GPCR activated by glutamate. Its ligand binding causes a conformation change that triggers signaling via guanine nucleotide-binding proteins (G proteins) and modulates the activity of downstream effectors in human (Wu et al., 1998; Song et al., 2021). U258 is a probable GPCR Mth-like 4 gene, which were found negatively regulating oxidative stress in Tribolium castaneum (Li et al., 2014). Besides, methuselah (Mth) is associated with lifespan, stress resistance, and reproduction in Drosophila melanogaster (Araujo et al., 2013). U165 and U196 are both GPCR Mth2 genes, which are involved in maintaining genome stability in human cells against oxidative stress (Hashiguchi et al., 2018). U33 and U120 were predicted to be GPCRs by GO annotation, while their specific functions have not been determined. The effects of these candidate GPCR genes can be further determined by gene knockout to analyze their specific functions during molting.

Although the effect of GCPRs on ecdysis is unclear in E. sinensis, it has been explored in other animals. Bai et al. identified 111 GPCRs through genome analysis of T. castaneum and screened 8 GPCR encoding genes through large-scale RNA interference experiments. The knockdown of these GPCRs resulted in serious development stagnation and molting failure, which shows that GPCRs are important for molting and development of T. castaneum. These eight GPCRs include D2R, cirl, Stan, mthl, FZ, SMO, bursicon, and ccapr (Bai et al., 2011), where mthl is encoded by U258 in this study. In addition, the 111 GPCRs identified by them also include several metabotropic glutamate receptors, which are encoded by U116 in this study. Zhao et al. found that GPCRs are related to the regulatory function of 20-hydroxyecdysone (20E) through the study on the Lepidoptera insect Helicoverpa armigera (Cai et al., 2014). GPCRs control the transduction of 20E in cell membrane, transmit 20E signal, and control 20E to enter cells (Wang et al., 2015). In another study, they found that 20E activates phospholipase Cγ1 via the GPCRs and cytosolic tyrosine kinases to regulate Ca²⁺ influx and protein kinase C phosphorylation of transcription factor ultraspiracle and subsequently influence gene transcription for molting and metamorphosis (Liu et al., 2014). These studies provide clues for the further study of the GCPRs in E. sinensis.

In addition, through RNA-seq analysis, the researchers also found the existence of GPCRs in the Y-organ of Gecarcinus lateralis (Tran et al., 2019) and the central nervous system tissue and Y-organ and epidermal tissue of Carcinus maenas (Oliphant et al., 2018), indicating that GPCR may be a widely existing regulatory element in crabs, but their regulatory mechanism in crabs is still unclear.