Transcriptome analysis elucidates key changes of pleon in the process of carcinization*

Ya’nan YANG, Zhaoxia CUI, Tianyi FENG, Chenchang BAO, Yuanfeng XU

School of Marine Sciences, Ningbo University, Ningbo 315832, China

Abstract When megalopa molting to the first juvenile crab stage, the crabs undergo carcinization morphogenesis. To study the key physiological and morphological processes in carcinization, we performed a comparative transcriptomic analysis between the cephalothoraxes and the pleons of megalopa and the first juvenile crab stage in Chinese mitten crab. The results reveal that the major physiological and morphological changes in the pleon were related to energy metabolism (oxidative phosphorylation and AMPK pathways),ventral nerve cord fusion (apoptosis-related pathways), and metamorphosis (transcription factors, Hedgehog and Hippo pathways). We also discovered that the key Hox genes abdominal-B and abdominal-A might regulate morphological changes, especially in the degeneration of the fifth pair of pleopods, and ganglion fusion, respectively. Studying the regulatory mechanisms of carcinization may help us better understand the developmental biology of the juvenile crabs.

Keyword: Eriocheir sinensis; carcinization; pleon; physiology; morphology

1 INTRODUCTION

Metamorphosis is an important developmental transition in a large number of animal taxa.Metamorphosis refers to the process in which an immature individual undergoes drastic anatomical and physiological changes to develop into an adult in the life cycle (Medina, 2009). In the brachyuran decapod, crabs go through a carcinization morphogenesis process, involving remarkable morphological and anatomical transformations. In freshwater crabs, these occur embyonically, since all true freshwater crabs lack free-swimming larval stages and exhibit direct development (Davie et al.,2015). Due to carcinization, crabs have the shortened and compressed pleon, folded beneath the cephalothorax, and inserted between the pereiopods or in a special cavity (McLaughlin and Lemaitre,1997; Guinot and Bouchard, 1998; Scholtz, 2014),the degeneration of pleopods (Lee et al., 1994), and the fusion of ventral nerve cords (Davie et al., 2015).In some crabs, the secondary sexual characteristics appear in this process, including the formation of gonopores and gonopods (Payen, 1975;Lee et al., 1994; Davie et al., 2015).

The Chinese mitten crab,Eriocheirsinensisis one of the most commercially important aquacultural species in China (FAO, 2020). Thanks to the development of hatchery and larviculture techniques,larvae are available from captive cultures, and aquaculture is not reliant on wild sources; as a result,the Chinese mitten crab aquaculture spreads nearly all over China (Sui et al., 2011). The larval development ofE.sinensisincludes five zoeal stages and one megalopa stage before the megalopa molts to the first juvenile crab stage (Anger, 1991). When the megalopa molting to the first juvenile crab stage, the crabs undergo carcinization morphogenesis, such as the reduction and folding of the pleon beneath the cephalothorax (Guinot and Bouchard, 1998) and the fusion of ventral nerve cord and thoracic ganglion(Song et al., 2017). To date, carcinization ofE.sinensishas been studied at a transcriptomic level (Song et al.,2015), and some genes such as proliferating cell nuclear antigen (PCNA), cuticle and ferritin had been identified in this process (Li et al., 2011). It remains unclear whether the described genes are related to the key physiological and morphological processes in carcinization morphogenesis.

To study the molecular mechanism of pleon metamorphogenesis during carcinization, we performed a comparative transcriptomic analysis for the cephalothoraxes and pleons of the megalopa and the first juvenile crab stage. Additionally, we performed an enrichment analysis on the key physiological and morphological changes, including energy metabolism, morphological changes of pleon,and ventral nerve cord fusion. Specifically, our analysis focused on how the Hox genes, abdominal-A(abd-A) and abdominal-B (abd-B), regulate the carcinization. Researching the regulatory mechanisms of carcinization process may help us better understand the developmental biology ofE.sinensis.

2 MATERIAL AND METHOD

2.1 Sample preparation and morphological observation of pleons

Megalopaes (M) and the first juvenile crab stage(J) ofE.sinensiswere collected from a farm in Nantong, Jiangsu province, China. For molecular studies, the cephalothoraxes and pleons were frozen in liquid nitrogen immediately. Three biological replicates were performed in cephalothoraxes of megalopaes (MC) and the first juvenile crab stage(JC), and in pleons of megalopaes (MP) and the first juvenile crab stage (JP) respectively, receiving a total of 12 samples.

For morphological observations, the external characteristics (pleopods and pleon shape) of exuviaes from megalopae (n=20) and whole body from the first juvenile crab stage (n=20), were photographed using Olympus SZ51 microscopes.

2.2 cDNA library construction and sequencing

Total RNA was isolated with the Trizol Reagent(Invitrogen, USA) according to the manufacturer's instructions. RNA quality and concentration were determined using a Qubit®RNA Assay Kit (Life Technologies, CA, USA) and a NanoDrop 6000(Agilent Technologies, Palo Alto, CA, USA). mRNA poly (A) was separated using oligo (dT) beads, and then digested into short fragments with a TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA,USA). Double strands cDNA were synthesized with random hexamer primers. The cDNA fragment and sequencing adaptor were ligated with T4 DNA ligase according to Illumina manufacturer’s protocol. The ligated products were amplified to create a cDNA library. Each cDNA library was sequenced by the Illumina HiSeq™ 2000 platform.

2.3 Transcriptome assembly and annotation

Clean reads were guaranteed by removing adapter,filtering low-quality sequences (＜Q20) and sequences shorter than 50 bp. The clean reads were assembled using Trinity software (http://trinityrnaseq.sourceforge.net/) as described for de novo transcriptome assembly without a reference genome (Garber et al., 2011).

Gene functional annotations were performed basing on public databases. Each transcript was compared using Blastx algorithm (E≤1e-5) with the National Center for Biotechnology Information (NCBI) nonredundant protein sequence (Nr) database (http://www.ncbi.nlm.nih.gov/) and clustered into a unique sequence (unigene). Next, Hmmscan (HMMER 3.0)was used to search for Protein family (Pfam) domain signatures in unigenes. Moreover, animal transcription factors (TFs) database AnimalTFDB version 2.0(http://bioinfo.life.hust.edu.cn/AnimalTFDB/) was used to identify TFs. Gene Ontology (GO) (http://www.geneontology.org/) annotations and functional classification of all unigenes were obtained by using Blast2GO. The biochemical pathways were predicted by Kyoto Encyclopedia of Genes and Genomes(KEGG) (http://www.genome.jp/kegg/). Furthermore,NCBI nucleotide sequences database (Nt), SwissProt database, and eukaryotic Orhtolog Group database(KOG) were also used to annotate unigenes.

2.4 Gene expression analysis

To detect the expression levels of each unigene, all clean reads in each sample were first aligned with Unigene database using Bowtie (Langmead, 2010)before performing the estimation of expression level by combined RNA-Seq by expectation maximization(RSEM) (Li and Dewey, 2011). The fragments per kilobase of transcript sequence per million mapped reads (FPKM) (Trapnell et al., 2010)value was used to estimate the unigene expression abundance.

2.5 Identification and validation of diff erentially expressed genes

Diff erentially expressed genes (DEGs) between two transcriptomes were identified by the DESeq program (Love et al., 2014). DEGs were screened outaccording to the diff erence of expression amount fold change (FC≥2) and significantP-value (Padj＜0.05).GO enrichment analysis of DEGs was conducted by using the topGO software. KOG statistical classification and KEGG annotation of DEGs were also performed. The significance of enriched pathway was analyzed using the method of Fisher’s exact test(Tavazoie et al., 1999). Functional domain analysis was performed with ExPASy PROSITE (http:// www.expasy.ch/tools/scanprosite/).

Table 1 Summary of assembly statistics in E. sinensis

Total RNA were isolated from MC, MP, JC, and JP as described above. Approximately 2-μg RNA were used to the reverse transcription of cDNA. The fluorescent quantitative real-time PCR (qRT-PCR)was used to validate the expressed transcripts of DEGs. DEGs in Hedgehog and Hippo pathways included abd-A, abd-B, p53, and Dachs. Meanwhile,the randomly selected genes included hypothetical protein T11_4388 (HP), prophenoloxidase-activating factor (PPAF), fatty acid binding protein (FABP),crustin, hemocyanin subunit 6 (HS) and cuticle proprotein proCP2.2 (CP). Gene specific primers were designed according to the Illumina sequencing data and the housekeeping gene β-actin was used as an internal control (Supplementary Table S1). The qRT-PCR was carried out with 2×SYBR Select Master Mix (Applied Biosystems, Austin, TX, USA)as the manufacturer’s instructions in a 7500 Real-time PCR System (Applied Biosystems, Carlsbad, CA USA). All experiments were performed in triplicate.Gene expression level relative to control was determined by the 2-ΔΔCtmethod (Livak and Schmittgen, 2001). All data were shown as mean±SD.Statistical analysis was performed using IBM SPSS 20. The data obtained from qRT-PCR analysis were subjected to One-way Analysis of Variance (ANOVA)followed by Dunnett’s 2-sided post-hoc test.Significance was set atP＜0.05.

3 RESULT

3.1 External characteristics of pleons

After the M molted to become the J stage, the pleon was shorter, compressed, and folded beneath the cephalothorax by morphometry. The fifth pair of pleopods was smaller in M and completely reduced in J stage, and only four pairs of pleopods were observed in the pleon (Fig.1). However, there are no external sexual characteristics between male and female crabs,in M and J stage.

3.2 Transcriptome sequencing and read assembly

To cover the transcriptome of cephalothoraxes and pleons in megalopaes and the first juvenile crab stage(Fig.1), 12 cDNA libraries (three biological replicates)ofE.sinensiswere sequenced using Illumina Hiseq 2500 platform. The average raw reads for MC, MP,JC, and JP were 55 921 961, 52 237 985, 45 869 354,and 47 598 691, respectively. The cycleQ30 percentage of all samples was more than 91.69% (Table 1). A total of 267 158 transcripts were assembled, with a mean length of 1 123 bp using Trinity and Corset software.From these transcripts, 131 138 unigenes were assembled, with an average length of 958 bp(Supplementary Fig.S1).

Fig.1 Megalopaes and the first juvenile crab stage of E. sinensis

3.3 Functional annotation

Seven diff erent databases, including Nr, Nt,KEGG, SwissProt, Pfam, GO, and KOG were used to annotate the function of 131 138 unigenes.Among the unigenes, 63 303 were annotated with at least one database, and 2 592 were annotated in all seven databases. A summary of functional annotation of contigs is showed in Supplementary Table S2.

GO enrichment analysis (Supplementary Fig.S2)showed that 45 736 unigenes were associated with 56 GO terms, 26 of which were involved in biological processes (35 167 unigenes), 10 in molecular functions (36 721 unigenes) and 20 in cellular components (26 892 unigenes). The most prominent terms for biological processes, molecular functions,and cellular components are presented in Fig.2. KOG analysis showed that 15 974 unigenes were assigned to 25 KOG terms. The top four KOG terms were general function prediction (2 156 unigenes);translation, ribosomal structure and biogenesis(2 138 unigenes); posttranslational modification,protein turnover, chaperones (1 944 unigenes); and signal transduction mechanisms (1 762 unigenes)(Supplementary Fig.S3). KEGG analysis revealed that 7 658 unigenes assigned to 5 categories which were comprised of 229 KEGG pathways. The most prominent KEGG pathways were related to ribosome(610 unigenes), carbon metabolism (234 unigenes),phagosome (193 unigenes), oxidative phosphorylation(172 unigenes), purine metabolism (157 unigenes),protein processing in endoplasmic reticulum(153 unigenes) and biosynthesis of amino acids(152 unigenes). Based on the Pfam database,1 040 unigenes were assigned to TFs, and the main TF families were zf-C2H2 (442 unigenes), ZBTB(113 unigenes), Homeobox (79 unigenes), and bHLH(66 unigenes).

3.4 The analysis of diff erentially expressed genes

The analysis of four transcriptome data revealed 27 023 DEGs. The comparative analysis identified 9 309 DEGs (5 538 up-regulated and 3 771 downregulated) between MC and MP, 9 090 (6 650 upregulated and 2 440 down-regulated) between JC and JP, 12 760 (1 670 up-regulated and 11 090 downregulated) between MC and JC, and 16 907 (4 564 up-regulated and 12 343 down-regulated) between MP and JP. Some DEGs from four pairwise comparisons are listed in Supplementary Table S3.The distribution of DEGs was illustrated in volcano plots (Supplementary Fig.S4). Based on the GO analysis, the DEGs were mainly assigned to three functional groups: biological process, cellular component, and molecular function (Supplementary Fig.S5). In each comparison, top 10 GO terms of the most up-regulated and down-regulated DEGs are listed in Table 2. The most prominent terms for biological processes were cellular process and metabolic process; for molecular functions were binding and catalytic activity; for cellular components were membrane and cell. The top 20 KEGG pathway of the four pairwise comparisons are shown in Supplementary Fig.S6.

3.5 Up-regulated DEGs from MP to JP

A homeobox TF and a fork-head TF, which are related to axial patterning, were identified among 91 dramatically up-regulated TFs during carcinization in JP compared with MP (Supplementary Table S4). In addition, a number of multifunctional TFs were also up-regulated, including the E26 transforming-specific protein (ETS), AT-rich interaction domain (ARID),E2F TFs, MYB TFs, P53 TFs, the basic leucine zipper protein of TF family (TF-bZIP), the thanatosassociated protein (THAP), basic helix-loophelix protein (bHLH) and LPS-induced tumor necrosis factor alpha factor (zf-LITAF-like). They regulate biological processes such as cell cycle,cellular diff erentiation, cell proliferation and apoptosis. Additionally, a Doublesex/MAB-3 domain(DM) was significantly up-regulated in JP.

Table 2 The top 10 GO terms with significant DEGs from four pair-wise comparisons

Table 3 The top 20 KEGG pathways with significant DEGs in MP-to-JP

The most up-regulated top 20 KEGG pathways in MP-to-JP comparison are listed in Table 3. Regarding the energy metabolism pathways, unigenes involved in 5′-AMP-activated protein kinase (AMPK) signaling pathway were up-regulated (Fig.3a). In the Hedgehog pathway (morphogenesis pathway), PKA, casein kinase 1 (CK1) and glycogen synthase kinase 3 beta(GSK3β) were up-regulated in JP (Fig.4a). Meanwhile,in the Hippo pathway, discs overgrown (Dco),approximated (App), MOB kinase activator 1 (Mats)and p53 were up-regulated in JP (Fig.4b).

3.6 Down-regulated DEGs from MP to JP

The most down-regulated top 20 KEGG pathways in MP-to-JP comparison are also listed in Table 3.Unigenes involved in another energy metabolism pathways oxidative phosphorylation pathway were dramatically down-regulated (Fig.3b). In the Hedgehog pathway, patched (Ptc) was down-regulated in JP (Fig.4a). As to the Hippo pathway, Dachs and keren (Krn), were down-regulated in JP (Fig.4b).Unigenes involved in apoptosis related pathways,including parkinson’s disease, huntington’s disease,alzheimer’s disease in humans and apoptosis pathways were dramatically down-regulated in JP compared with MP (Fig.5).

3.7 The DEGs of abd-A and abd-B

Analysis of the DEGs revealed that abd-A was higher in MP (P＜0.01) than in MC, and abd-B was higher in JP (P＜0.01) than in JC (Supplementary Table S3).

3.8 Validation of Illumina sequencing results by qRT-PCR

The expression levels of key DEGs in Hedgehog and Hippo pathways and other selected genes shown in qRT-PCR showed a high consistency with the results of RNA-Seq (Fig.6).

Fig.3 Putative energy metabolism pathways in MP-to-JP comparisons

Fig.4 Putative morphological change pathways in MP-to-JP comparison

Fig.5 Putative the apoptosis and apoptosis-related pathways in MP-to-JP comparisons

4 DISCUSSION

Crabs undergo morphological changes after molting to the first juvenile crab stage. Compared to the megalopae with five pairs of pleopos (Montú et al., 1996), the pleopods were degenerate, only four pairs of pleopods in J stage. The development of pleopods inE.sinensisis similar to that in the Japanese mitten crab,Eriocheirjaponicus(Lee et al.,1994) and the mud crabScyllaparamamosain(Jiang et al., 2020). High-throughput sequencing is a newly developed technology for discovering and profiling genes effi ciently (Reuter et al., 2015). A previous study used the suppression subtractive hybridization(SSH) method to indetify the DEGs during the carcinization of the Chinese mitten crab assembled only 325 unigenes (Li et al., 2011), but with highthroughput sequencing, 16 907 diff erent expression genes were assembled between MP and JP. From a comparison with SSH, the RNA-Seq generated more data. The most prominent KEGG pathways in the development of M to J stage (MC-to-JC and MP-to-JP) ofE.sinensiswere similar to the metamorphosis research in the prawnMacrobrachiumrosenbergii(Ventura et al., 2013), which also undergoes morphological change from planktonic larval stages to the benthic post-larval stage (Abdu et al., 1998).

4.1 Energy metabolism during carcinization

Fig.6 qRT-PCR validation of the selected DEGs in RNA-Seq

In the oxidative phosphorylation pathway, NADH dehydrogenase, cytochrome c reductase, cytochrome c oxidase and F-type ATPase work together to generate energy (Dimroth et al., 1999; Heikal et al.,2014). The AMPK pathway is a canonical route regulating energy homeostasis of the whole body(López et al., 2016). In this study, the oxidative phosphorylation pathway had the most downregulated genes, while genes related to AMPK signaling pathway were up-regulated, in MP compared with JP. This result suggests that the energy metabolism pathway have changed in the pleon from the M to J stage, the main pathway was the AMPK pathway rather than the oxidative phosphorylation pathway in J stage, and the mechanism of energy utilization have changed. The pleon of megalopae crabs remains extended for better swimming ability(Atkins, 1954), while the pleon of the juvenile mitten crab is tucked under the cephalothorax in an adaptation to a benthic lifestyle, now protecting reproductive organs for mating and carrying eggs in adult female crabs (Števčić, 1971; Herborg et al., 2006; Davie et al., 2015). It is likely that the changed energy metabolism pathways of the pleon coincide with the adjustment of swimming ability, as describe in the bryozoanBugulaneritina(Wong et al., 2014), in which the level of energy production decreased, when the larvae ceased swimming.

4.2 Morphological changes of pleon in J stage

TFs have a leading role in the initiation of gene expression (Spitz and Furlong, 2012). In this study,91 TFs were dramatically up-regulated in JP compared with MP, among which 13 zf-C2H2 TFs were identified. The Homeobox TF and fork-head TF play an important role in the development of organisms,such asStrongylocentrotuspurpuratus(Howard-Ashby et al., 2006) andXenopuslaevis(Rohl and Knöchel, 2005). As a class of zinc finger motifs, zf-C2H2 TFs execute important regulatory functions in the development ofStrongylocentrotuspurpuratus(Materna et al., 2006). Therefore, we deduced that zf-C2H2 TFs as well as Homeobox TF and fork-head TF may be involved in the pleon development from M to J stage expecially in the morphological change by regulating target gene expression. Futhermore, a DM was also up-regulated in MP-to-JP comparisons. DM domain has been shown to be involved in sex-specific diff erentiation ofCaenorhabditiselegansand the medaka,Oryziaslatipes(Lints and Emmons, 2002;Nagahama, 2005). In the brachyuran decapods, the diff erentiation of gonopores and pleopods take place simultaneously from the first to the third juvenile crab stage (Lee et al., 1994). Therefore, we propose that the up-regulation of DM domain might correlate to the diff erentiation of gonopores or pleopods in the early juveniles.

Hedgehog signaling pathway and Hippo signaling pathway-fly that up-regulated in JP compared with MP, are implicated in the regulation of metamorphosis of many metazoan groups (Ingham and McMahon,2001; Pan, 2007). The Hedgehog pathway is involved in organizing growth and patterning in many developmental processes (Mohler, 1988) and the Hippo pathway coordinates with cell death,diff erentiation and inhibition of cell proliferation (Yu and Guan, 2013). Our transcriptomic analysis showed that the DEGs in MP-to-JP comparison were assigned to Hedgehog pathway, which was similar to the previous study of the carcinization ofE.sinensis(Song et al., 2015). Additionally, we found the Hippo pathway might take part in the carcinization of the Chinese mitten crab in this research.

The Hox gene abd-B plays an important role in appendage development. For example, it suppresses the development oflepidopteran proleg in the posterior abdomen and suppresses the development of appendages in genus ofDrosophila(Estrada and Sánchez-Herrero, 2001; Tomita and Kikuchi, 2009).Additionally, the interaction between abd-B and Hedgehog pathway regulates the chick hindgut andDrosophilagenital disc development (Roberts et al.,1995; Sánchez et al., 2001). Analysis of the DEGs showed that the abd-B was higher in JP (P＜0.01) than in JC. Hence, we deduced that abd-B might participate in the degeneration of the fifth pair of pleopods during carcinization by regulating the Hedgehog signaling pathway.

4.3 The fusion of ventral nerve cord during carcinization

Apoptosis or programmed cell death (PCD) is a conserved process among species (Vaux et al., 1994).We found that the unigenes involved in apoptosis related pathways were dramatically down-regulated in JP compared with MP. In carcinization, the ventral nerve cord fuses with the thoracic ganglion to form the thoracic ganglion in the thorax (Davie et al.,2015). In the Chinese mitten crab,E.sinensis, the fusion of ventral nerve cord and thoracic ganglion begins in M stage and completes in J stage (Song et al., 2017). Apoptosis-related pathway including,parkinson’s disease, huntington’s disease, and alzheimer’s disease in humans, as well as apoptosis pathway are related to central nervous system development and neurodegeneration (Schmidt et al.,1997; Federico et al., 2012). Furthermore, abdominal neuroblasts stop dividing at mid-larval stages due to the activation of apoptosis, which is mediated by the Hox gene abd-A, inDrosophila(Bello et al., 2003).Analysis of the transcriptome data, we found that the abd-A was higher in MP (P＜0.01) than in MC, which suggested the abd-A might involve in the fusion of ventral nerve cord during carcinization.

5 CONCLUSION

In this study, we reported the transcriptomic changes during megalopa molting to the first juvenile crab stage ofE.sinensis. Our results reveal major transcriptomic changes of pleons in terms of the oxidative phosphorylation and the AMPK pathways of the energy metabolism, the TFs, the Hedgehog and Hippo pathway of the morphological changes, and apoptosis-related pathways during the fusion of ventral nerve cord. By studying the key Hox genes abd-B and abd-A in carcinization, we found that the abd-B might participate in the degeneration of pleopod by regulating the Hedgehog signalling pathway-fly and the abd-A might regulate the ganglion fusion.

6 DATA AVAILABILITY STATEMENT

All the original Illumina HiSeq sequencing data are available from the NCBI Sequence Read Archive(SRA) database (PRJNA644959). Other data that support the current study are available from the corresponding author on reasonable request.