Comparative transcriptome analysis explores maternal to zygotic transition during Eriocheir sinensis early embryogenesis

Accepted Manuscript Comparative transcriptome analysis explores maternal to zygotic transition during Eriocheir sinensis early embryogenesis

Junhao Ning, Chengwen Song, Danli Luo, Yuan Liu, Hourong Liu, Zhaoxia Cui PII: DOI: Reference:

S0378-1119(18)31070-9 doi:10.1016/j.gene.2018.10.036 GENE 43292

To appear in:

Gene

Received date: Revised date: Accepted date:

21 April 2018 27 September 2018 11 October 2018

Please cite this article as: Junhao Ning, Chengwen Song, Danli Luo, Yuan Liu, Hourong Liu, Zhaoxia Cui , Comparative transcriptome analysis explores maternal to zygotic transition during Eriocheir sinensis early embryogenesis. Gene (2018), doi:10.1016/ j.gene.2018.10.036

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Comparative transcriptome analysis explores maternal to zygotic transition during Eriocheir sinensis early embryogenesis

Junhao Ninga,c,d, Chengwen Songa,b,c, Danli Luoa,c,d, Yuan Liua,b,c, Hourong Liua,c,d,

SC

RI

PT

Zhaoxia Cuia,b,c*

a. CAS Key Laboratory of Experimental Marine Biology, Institute of Oceanology,

NU

Chinese Academy of Sciences, Qingdao 266071, China

MA

b. Laboratory for Marine Biology and Biotechnology, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China

PT E

D

c. Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China

AC

CE

d. University of Chinese Academy of Sciences, Beijing 100049, China

* Corresponding authors.

1

ACCEPTED MANUSCRIPT Address: Institute of Oceanology, Chinese Academy of Sciences, 7 Nanhai Road, Qingdao 266071, China

Tel: +86 532 82898509.

PT

Fax: +86 532 82898509.

SC

RI

E-mail address: [email protected] (Z. Cui).

NU

Abstract

MA

The maternal genome directs almost all aspects of early animal development. As development proceeds, the elimination of maternal gene products and zygotic genome

D

activation (ZGA) occur during the maternal to zygotic transition (MZT). To study the

PT E

molecular mechanisms regulating this developmental event in Eriocheir sinensis, RNA-Seq technology was applied to generate comprehensive information on

CE

transcriptome dynamics during early embryonic stages. In total, 32,088 annotated

AC

unigenes were obtained from the transcriptomes of fertilized eggs and embryos at the cleavage (2-4 cell) and blastula stage. A total of 566 maternal genes and 1165 zygotic genes were isolated, among which 103 and 266 genes were predicted conserved maternal transcripts (COMATs) and conserved zygotic transcripts (COZYTs), respectively. The COMATs performed housekeeping gene functions and may be essential for initiating early embryogenesis of the Bilateria. Furthermore, 87, 76 and 117 differentially expressed genes associated with the MZT, morphogenesis and 2

ACCEPTED MANUSCRIPT immunity were identified when compared the three transcriptomic datasets. We also unmask that the MZT takes place around the cleavage stage, when the genes involved in the clearance of maternal gene products and the ZGA were significantly up-regulated. Taken together, these datasets provide a valuable resource for

PT

understanding the mechanisms of early developmental events in E. sinensis, and

RI

facilitate further studies on molecular mechanisms of asynchronous development in

SC

crabs.

NU

Keywords: mitten crabs; embryogenesis; maternal to zygotic transition; transcriptomics

MA

1. Introduction

D

The Chinese mitten crab (Eriocheir sinensis) is one of the most commercially

PT E

important aquatic species in China. This crab has contributed obviously to the understanding of invertebrate development and functional genomics due to its

CE

available resources, such as typical zoeal larval and metamorphotic stages (Montú et al., 1996). Nevertheless, cannibalism owing to asynchronous development in

AC

crustacean has been restricting the development of E. sinensis industry (Lovrich, 1997; Marshall et al., 2005; Moksnesa, 1997). Asynchronous development of crab embryos will deteriorate sharply from late cleavage stage (32 cell) to early gastrula stage in natural conditions, which causes serious cannibalism and reduces the yield of juvenile crabs (Zhao et al., 1993). Although early embryonic development of E.

3

ACCEPTED MANUSCRIPT sinensis has been reported previously (Huang et al., 2011), little is known about the molecular mechanism of the asynchronous development.

Embryonic development is controlled by a complicated interaction between maternal and zygotic provisioned transcripts. Maternally deposited products (mRNAs

PT

and proteins) control virtually all aspects of early embryonic development prior to the

RI

zygotic genome activation (ZGA) (Korzh, 2009; Tadros and Lipshitz, 2009). With the

SC

development of animal embryos, the ‘maternal-to-zygotic transition’ (MZT) is

NU

triggered by the concomitant degradation of maternal products and the activation of zygotic transcription (Tadros and Lipshitz, 2009; Walser and Lipshitz, 2011).

MA

Meanwhile, Asynchronous cleavages has been reported to be associated with the formation of cell interphase during the MZT (Tadros and Lipshitz, 2009). In recent

D

years, the MZT has been extensively studied in model organisms, including

PT E

echinoderms, nematodes, insects, fishes, amphibians and mammals, but the timing and scale of this event varies from 2-cell to the gastrula (Tadros and Lipshitz, 2009).

CE

For example, the first cohort of zygotic genome is activated at the 10th cell division in

AC

zebrafish, while at 13th cell division in Drosophila (Kimmel et al., 1995). Although several researches have tried to reveal the regulating mechanisms of maternal transcripts degradation as well as zygotic gene activation, most analyses were focused on single transcript or small subset of transcripts. Accordingly, studies on the molecular mechanism of the MZT may be conducive to further elucidate E. sinensis embryonic development and reduce asynchronous rate. 4

ACCEPTED MANUSCRIPT Advances in high-throughput sequencing technologies have had a tremendous impact on genomics (Green et al., 2010), transcriptomics (Sultan et al., 2008), and stem cell biology (Tang et al., 2010). RNA-sequencing (RNA-Seq) has offered a high-effective and convenient tool for analysis of gene expression, identification of

PT

novel genes and differentially expressed genes (Garber et al., 2011). Studies on

RI

transcriptomes of early embryogenesis have been reported in different invertebrates,

SC

such as D. melanogaster (De Renzis et al., 2007), Caenorhabditis elegans (Baugh et al., 2003), Oncopeltus fasciatus (Ewen-Campen et al., 2011), and Lymnaea stagnalis

NU

(Liu et al., 2014). An extensive analysis of transcriptome dynamics during

MA

invertebrate embryogenesis would provide insights into the molecular regulation mechanism of the MZT. Here, we apply RNA-seq to gain a comprehensive

D

understanding of transcriptional processes occurring from the fertilized egg to the

PT E

blastula. Our focus is on the major developmental events during the MZT, such as axis formation, segmentation and immunity. The work builds a valuable resource for

AC

CE

E. sinensis developmental biology and functional genomics.

2. Materials and Methods 2.1. Sample preparation The adult mitten crabs were collected during their reproductive period from a farm in Panjin, Liaoning Province, China. Females and males were placed in one container filled with filtered seawater (salinity 20) for mating at 15 °C. Then the fertilized eggs 5

ACCEPTED MANUSCRIPT (Fe, collected immediately after spawning) and embryos from the cleavage stage (Cs, 2-4 cell) and blastula stage (Bs) were collected from healthy ovigerous crabs. The developmental stages of embryos were monitored under a dissecting microscope, and embryo samples at each stage were collected from five identical crabs. All samples

PT

were collected separately in 1.5 mL tubes, then immediately frozen in liquid nitrogen

NU

2.2. cDNA library construction and sequencing

SC

RI

until RNA extraction.

MA

Total RNA was extracted from frozen Fs, Cs and Bs embryos using Trizol Reagent (Invitrogen, USA). Residual DNA was removed from the extracted RNA samples

D

with RNase-free DNase I. RNA quality and concentration were determined using a

PT E

Qubit®RNA Assay Kit with a Qubit®2.0 Fluorimeter (Life Technologies, CA, USA) and a NanoDrop 6000 Assay Kit with an Agilent Bioanalyzer 2100 system (Agilent

CE

Technologies, Palo Alto, CA, USA). The equal amounts of RNA samples from five

AC

replicate crabs were pooled as input material for RNA-Seq. mRNA was purified using magnetic oligo (dT) beads, and then digested into short fragments with a TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA). The short mRNA fragments, mixed with random hexamer primers, were used as templates to synthesize first-strand cDNA. Then double-stranded cDNA was prepared and sequencing adaptors were ligated with T4 DNA ligase according to Illumina manufacturer’s protocol. The ligated products were amplified to create a final cDNA library after purification with 6

ACCEPTED MANUSCRIPT AMPureXP beads (Beckman Coulter, High Wycombe, United Kingdom). Then each cDNA library was sequenced by the Illumina HiSeqTM 2000 platform and paired-end

PT

reads were generated.

2.3. Transcriptome assembly and annotation

RI

All raw reads were deposited in the NCBI Short Read Archive (SRA) database

SC

(http://www.ncbi.nlm.nih.gov/Traces/sra/). After removing adapter sequences,

NU

high-quality clean reads were obtained by filtering the low-quality sequences (
MA

using Trinity software (http://trinityrnaseq.sourceforge.net/) as described for de novo

D

transcriptome assembly without a reference genome (Garber et al., 2011). The Fe, Cs

PT E

and Bs datasets were analyzed separately first, and then assembled and reanalyzed.

Gene annotations were performed by sequence alignment with public databases.

CE

Each transcript was searched against the National Center for Biotechnology

AC

Information (NCBI) non-redundant protein sequence (Nr) database (http://www.ncbi.nlm.nih.gov/) and clustered into a unique sequence (unigene). After that, GO (Gene Ontology; http://www.geneontology.org/) annotations were assigned to classify the potential functions of all unigenes using Blast2GO. The biochemical pathway were generated by Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/), and non-supervised Orthologous Groups (eggNOG)

7

ACCEPTED MANUSCRIPT (http://eggnog.embl.de/) were performed to predict and classify potential functions of the unigenes based on known orthologous gene products.

PT

2.4. Gene expression analysis To detect gene expression levels, clean reads in each sample were first mapped to

RI

the assembled transcriptome to obtain the read number of each gene, and then the

SC

RPKM (reads per kilobase of exon model per million mapped reads) value was

NU

calculated to compare gene expression differences between different samples (Mortazavi et al., 2008). A RPKM threshold value of 0.1 was set to detect the

MA

presence of a unigene, which corresponded to a false discovery rate (FDR) of less

PT E

D

than 0.05

2.5. Differentially expressed genes identification, enrichment and pathway analysis

CE

Differentially expressed genes (DEGs) between two transcriptomes were identified

AC

by the DESeq program (http://www-huber.embl.de/users/anders/DESeq/) (Anders and Huber, 2010). Unigenes with a significant adjusted P-value (P<0.05) and a foldchange value >2 were considered to be diff erentially expressed between two sample groups (Anders and Huber, 2010). Classification of genes by time of increase or decrease in abundance is relevant to their regulation and function (Baugh et al., 2003). Four dynamic expression classes (M: maternal genes, Z: zygotic genes, MZ: 8

ACCEPTED MANUSCRIPT maternal-zygotic genes and T: transient genes) were divided by the defining features of each expression profile. After that, GO, eggNOG, KEGG Orthology (KO) and KEGG pathway enrichment analysis were performed to categorize the DEGs and detected potential pathways that might be involved in. In the GO and KEGG pathway

PT

enrichment analysis, a P-value of 0.05 (-Log10 (0.05) = 1.3) was set as the threshold

RI

for a significant diff erence between each two transcriptomes. The KEGG metabolic

SC

pathway of the unigenes was constructed online

(http://www.genome.jp/kegg/tool/map_pathway2.html), and the gene expression level

NU

was represented with different color. According to the public databases and previous

MA

literatures, crucial DEGs related to morphological change, hormone, nervous system

PT E

D

and nutrition metabolism were further manually checked.

2.6. Validation by quantitative real-time PCR

CE

The expression of eleven key DEGs was validated by quantitative real-time PCR (RT-qPCR) analysis. Primers designed according to the Illumina sequencing data are

AC

listed in Table 1. Total RNA was extracted from the same samples as that in Illumina sequencing, and first-strand cDNA was synthesized using a PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). The RT-qPCR was performed with an ABI 7500 Sequence Detection System (Applied Biosystems, Foster City, CA, USA). The β-actin gene was used as an internal control and all experiments were performed in triplicate. The reaction was carried out in a total volume of 10 μL, 9

ACCEPTED MANUSCRIPT containing 5 μL of 2× SYBR Premix Ex TaqTM II (TaKaRa), 0.2 μL of 50× ROX Reference Dye, 2 μL of diluted cDNA mix, 0.2 μL of each primer (10 mM) and 2.4 μL of Milli-Q water. Thermal profile for SYBR Green RT-qPCR was 95 ℃ for 30 s, followed by 40 cycles of 95 ℃ for 5 s and 60 ℃ for 35 s. The 2-∆∆CT method was used

PT

to analyze relative gene expression (Livak and Schmittgen, 2001). The results were

RI

subjected to one-way analysis of variance with SPSS 16.0, where P<0.05 indicated

NU

SC

statistically significant difference.

3. Results and Discussion

MA

3.1. Transcriptome sequencing and assembly

D

In total, 57,135,691, 36,954,342 and 32,813,420 raw reads were obtained from Fe,

PT E

Cs and Bs transcriptome, then deposited in NCBI Short Read Archive database with accession numbers SRR1182019, SRX1162680 and SRR2180019, respectively. After

CE

assembly by Trinity, 628,264 contigs were generated with an average length of 287,92 bp from the three transcriptomic datasets. Further assembly and BlastX searches

AC

against the Nr database resulted in 32,088 annotated unigenes, whose lengths ranged from 201 bp to 18,986 bp with an N50 length of 2,389 bp (Table 2). These data could greatly enrich the genetic resources for E. sinensis, especially for its developmental biology and functional genomics.

10

ACCEPTED MANUSCRIPT 3.2. Functional annotation of unigenes According to the KEGG analysis, 7,158 (22.30%) annotated sequences were mapped to 305 KEGG pathways. The eggNOG analysis showed that a total of 29,461 （91.81%）unigenes were categorized into 26 functional groups (Supplementary Fig.

PT

S1). Except the unigenes with ‘Functional unknown’ (14.16%) and ‘General function

RI

prediction only’ (17.19%), the dominant categories were ‘Signal transduction

SC

mechanisms’ (10.37%), followed by ‘Post-translational modification, protein turnover, chaperones’ (8.93%) and ‘Translation, ribosomal structure and biogenesis’

MA

NU

(7.07%).

D

3.3. Functional classification of differentially expressed gene

PT E

In total, 1,756 DEGs were obtained from the three transcriptomes. And the distribution of DEGs between two developmental stages was illustrated in volcano

CE

plots (Supplementary Fig. S2). Then we divided these DEGs into four categories by dynamic expression pattern (Baugh et al., 2003). 566 down-regulated genes from Fe

AC

to Bs (P ≤ 0.05) were named predicted maternal genes (M), 1165 up-regulated genes from Fe to Bs (P ≤ 0.05) were assigned to predicted zygotic genes (Z). Meanwhile, 2 predicted maternal-zygotic genes (MZ) and 23 predicted transient genes (T) were identified in our datasets. In order to verify whether “maternal-to-zygotic transition” occurs at a certain time between Fe and Bs, we compared the predicted maternal and zygotic genes of our 11

ACCEPTED MANUSCRIPT transcriptomes with five published embryonic transcriptomes (Table 3) (Azumi et al., 2007; Baugh et al., 2003; De Renzis et al., 2007; Nestorov et al., 2013). Through reciprocal tBLASTx analysis, we identified some quite conserved genes that have one or more homologous genes in these five species. According to the previous

PT

nomenclature (Liu et al., 2014), these conserved genes were named “conserved

RI

maternal transcripts” (COMATs) and “conserved zygotic transcripts” (COZYTs),

SC

respectively. 103 COMATs and 266 COZYTs were identified as the maternal-only and zygotic-only transcripts in the early transcriptomes of E. sinensis (Table 3), and

NU

the proportion of COMATs in E. sinensis was 18.20% (between 3.12 and 36.38%),

MA

implying maternally deposited genes were conserved across the invertebrates and may act as housekeeping genes involved in nucleotide binding functions, protein

D

degradation and activities associated with the cell cycle (Liu et al., 2014). However,

PT E

COZYTs accounted for 22.83% (well over 14.00%), indicating the zygotic transcripts were diverse across species owing to their different developmental regulation in

CE

embryogenesis.

AC

Based on the Gene Ontology (GO) analysis, the DEGs were mainly assigned to three functional groups ‘biological process’, ‘cellular component’ and ‘molecular function’ (Supplementary Fig. S3), and the most abundant GO terms in ‘biological process’ were ‘biosynthetic process’, ‘nucleic acid binging transcription factor activity’, and ‘translation’. To understand functional distribution of these DEGs, the KEGG enrichment analyses were performed and 594 DEGs were involved in 41 12

ACCEPTED MANUSCRIPT predicated metabolic pathways (Fig.1). The predicted maternal genes were mainly assigned to ‘signal transduction’ (24.41%), ‘translation’ (20.03%), ‘infectious diseases’ (17.34%), ‘carbohydrate metabolism’ (14.16%), ‘endocrine system’ (13.46%) and ‘cancers’ (12.46%) categories. The predicted zygotic genes were

PT

mainly allocated to ‘signal transduction’ (7.4%), ‘digestive system’ (6.23%),

SC

RI

‘endocrine system’ (6.06%), and ‘immune system’ (3.87%) categories.

NU

3.4. Validation of Illumina sequencing results by RT-qPCR

RT-qPCR was used to confirm the expression profiles of key maternal and zygotic

MA

genes during the MZT. As shown in Fig. 2, maternal genes, including ANK, CCNB,

D

BCDC and CCNA, exhibited high expression abundance in Fe, while zygotic genes,

PT E

Dicer1, RPBII, SODC, HB and FA11, were significantly up-regulated in Cs or Bs. Moreover, the transient gene NANOS reached the highest expression level in Cs,

CE

while the expression of maternal-zygotic gene CP decreased significantly in Cs when compared with that in Fe and Bs. RPBII and Dicer are the determinants for activation

AC

of the zygotic genome (Bettegowda et al., 2007), and their high expression after Fe indicated that the MZT virtually occurred from Cs to Bs. The RT-qPCR data was consistent with the results of RNA-Seq data, which not only validated the expression profile of the DEGs involved in the MZT, but also verified the accuracy and reliability of our transcriptome analysis.

13

ACCEPTED MANUSCRIPT 3.5. Candidate DEGs involved in the regulation of MZT According to the KEGG analysis, 87 DEGs involved in the MZT were isolated and divided into three classes, including 38 maternal protein degradation related genes, 23 maternal mRNA degradation related genes and 26 zygotic genome activation related

RI

PT

genes (Supplementary Table S1).

SC

3.5.1 Maternal protein degradation

NU

The degradation of destabilized maternal products is the first event of the MZT, which means most maternal products are important for oogenesis but not essential for

MA

embryogenesis. For example, in D. melanogaster, about 65% maternal mRNAs would

D

be eliminated during the MZT process (Lecuyer et al., 2007).

PT E

For the elimination of maternal proteins, autophagy and ubiquitin proteasome system are two key cellular regulatory systems (Wang and Klionsky, 2003), which are

CE

highly complex and mediated by lysosome and ubiquitin, respectively (Ciechanover,

AC

1998). In the current datasets, 27 DEGs were participated in autophagy, including lysosomal components cathepsin gene family (cathepsin C, cathepsin B, cathepsin L precursor, cathepsin L-like proteinase precursor, cathepsin L-like), glucosidase (glucosylceramidase-like, beta-glucuronidase, beta-N-acetylhexosaminidase) and nuclease (deoxyribonuclease II, lysosomal) (Supplementary Table S1). The above-mentioned DEGs were predicted zygotic genes and significantly up-regulated from Fe to Bs, which suggested the ZGA virtually takes place at Cs. Similar result has 14

ACCEPTED MANUSCRIPT been found in C. elegans indicating that the ZGA occurs at 2-cell stage (Baugh et al., 2003), whereas the ZGA occurs at fertilized egg stage in mouse and sea urchin (Aoki et al., 1997; Wei et al., 2006).

In addition, the DEGs associated with lysosomal enzyme transport like V-type

PT

H+-Transporting ATPase (V-type H+-Transporting ATPase subunit A, V-type

RI

H+-Transporting ATPase subunit B2, V-type proton ATPase 16 kDa proteolipid

SC

subunit), clathrin (clathrin heavy chain 1), receptor proteins (GABA (A)

NU

receptor-associated protein) were all up-regulated from Fe to Bs (Supplementary Table S1). However, the expression of sulfoglucosamine (arylsulfatase A-like,

MA

N-sulfoglucosamine sulfohydrolase) were significantly down-regulated. Besides, ubiquitin-proteasome associated DEGs including proteasome subunit (Proteasome

D

subunit beta type-7, Proteasome p44.5 subunit) and ubiquitin components (E3

PT E

ubiquitin-protein ligase 32, ADP-ribosylation factor) presented obviously up-regulated in Cs, which suggested that the zygotic genome have been activated at

CE

the cleavage stage and then might maintain or induce subsequent waves of the ZGA

AC

(Tadros and Lipshitz, 2009).

3.5.2 Maternal mRNAs degradation The elimination of maternal mRNAs is mainly accomplished by cis-elements and microRNAs (miRNAs) dependent degradation mechanism. For the cis-elements dependent degradation pathway, Smaug (SMG), a RNA-binding protein, is essential 15

ACCEPTED MANUSCRIPT for clearance of unstable maternal mRNAs (Tadros et al., 2007). Meanwhile, in D. melanogaster, Pumilio (PUM) as a cis-element is enriched in destabilized maternal transcripts (De Renzis et al., 2007). In this study, seven unigenes including zinc finger protein 36, C3H1 type-like 1-like, nanos homolog 2-like, hunchback Transcription

PT

factor, Translation initiation factor 4F, helicase subunit, cyclic AMP response

RI

element-binding protein a-like, involved in the SMG and PUM mediate degradation

SC

pathways were all up regulated in Cs. The target transcript cyclin B of SMG was significantly reduced from Cs to Bs, indicated that maternal mRNAs might begin to

NU

degrade in Cs. Additionally, Pan protein complex (PAB-dependent poly (A) -specific

MA

ribonuclease subunit 2, PAN3 poly (A) specific ribonuclease subunit homolog), Lsm protein complex (U6 snRNA-associated Sm-like protein LSm1-like, LSM4 homolog)

D

and RNA-degrading enzymes (2-phospho-D-glycerate hydro-lyase, chaperonin

PT E

GroEL) were significantly up regulated after the Fe stage.

An increasing evidences indicate that miRNAs are the mediators of maternal

CE

degradation pathway, such as miR-430 in zebrafish and miR-309 in D. melanogaster

AC

(Benoit et al., 2009; Giraldez et al., 2006). In the present study, RNase III Dicer, cleaving precursor mRNA to produce microRNA, was significantly up-regulated in Cs followed by down regulated in Bs. These results showed that maternal products have begun to degrade at the Cs stage in the Chinese mitten crab.

16

ACCEPTED MANUSCRIPT 3.5.3 Zygotic genome activation Though the scale and mechanism of the ZGA vary from species to species, several genes actively transcribed at the onset of the ZGA have been identified (Tadros and Lipshitz, 2009), and the advanced ZGA hypothesis has been reported in the model

PT

species. During early embryogenesis of E. sinensis, the unigenes participated in

RI

transcription and translation of zygotic gene, such as marker gene of zygotic

SC

activation (RNA polymerase II largest subunit) and transcription initiation factors (Transcription initiation factor TFIIF subunit beta, Transcription initiation factor

NU

TFIID TATA-box-binding protein), were all significantly up-regulated in Cs. While

MA

the genes related to “maternal clock” hypothesis (cyclin A-like, cyclin-A, G2 / mitotic-specific cyclin-A) were down-regulated from Fe to Bs. In addition, the genes

D

involved in “chromatin modification” hypothesis like methyltransferase (arginine

PT E

n-methyltransferase, betaine-homocysteine S-methyltransferase, methyltransferase, putative-like, methyltransferase-like protein 6, O-methyltransferase,

AC

CE

O-methyltransferase) were up-regulated from Cs to Bs (Supplementary Table S1).

3.6. Candidate DEGs involved in axis formation and segmentation As development proceeds, embryonic morphology changes slightly along the anterior–posterior (AP) and dorsal–ventral (DV) axis. In the present study, 21, 35 and 20 DEGs respectively involved in anterior–posterior axis formation, dorsal–ventral

17

ACCEPTED MANUSCRIPT axis formation and segmentation were discovered in E. sinensis (Supplementary Table S2).

For the AP axis patterning, bicoid (bcd) and nanos are two crucial transcription factors, and can promote the identification of anterior and posterior axis through

PT

regulating the expression of zygotic gene hunchback (hb) (Hulskamp et al., 1989;

RI

Irish et al., 1989; Lehmann and Nussleinvolhard, 1987, 1991; Nussleinvolhard et al.,

SC

1987; Struhl, 1989). However, bcd is only isolated in the higher diptera, such as

NU

Nasonia and Tribolium. Increasing evidences show that bcd is originated from Hox3 (Zen) gene, and the sequence and function of Hox3 are conserved in the short germ

MA

insects (Rosenberg et al., 2009). Hox3 and nanos were both detected in the present study, the former was up regulated from Fe to Bs, and the latter was significantly up

D

regulated in Cs followed by significantly down regulated, which indicats that cleavage

PT E

stage was a key period for crab AP axis formation. Similar work in D. melanogaster suggests that activated zygotic genes are responsible for rapid establishment of the

CE

body plan at the syncytium and blastula stages (De Renzis et al., 2007). The

AC

downstream Hox genes involved in AP axis construction were also detected and most of them were up regulated (Supplementary Table S2).

For the DV axis patterning, the BMP and Wnt signaling pathways are quite conserved in the establishment of DV axis (De Robertis, 2008). And Toll signaling pathway is mainly responsible for the formation of DV polarity in Drosophila (Anderson et al., 1985). According to our transcriptome datasets, 118, 23 and 39 18

ACCEPTED MANUSCRIPT unigenes involved in Wnt, BMP and Toll signaling pathways were detected, respectively. Of them, 12, 12 and 3 unigenes expressed differentially during the three stages (Supplementary Table S2), indicated that these DEGs may play an important role in crab DV axis formation. Furthermore, the position and expression levels of

PT

these genes in pathways were revealed in Fig. 4.

RI

In short-germ insects, the localized activation of serine protease cascade reaction in

SC

follicle cells would activate Spatzle gene, then activate the Toll signaling pathway in

NU

the ventral part of the perivitelline space (Moussian and Roth, 2005). In the present transcriptomes, most serine protease related genes were highly expressed in Fe and

MA

then down regulated obviously in Cs or Bs, which suggests they may play vital roles in early embryonic development of crab. As a sister clade of insects, crustacean

D

embryos formed by cell lineage rather than extracellular induced signaling pathways

PT E

(Price et al., 2010). Accordingly, the extracellular serine proteases (easter and snake) were predicted maternal genes. These embryos are not nursed by follicle cells, but

CE

attach to the pleopods of females by secreted egg stalks (Davis, 1964). Although the

AC

effect of egg stalks on the morphogenesis of embryo was still unclear, our observation provides a possibility that extracellular signals, such as Toll, might also participate in axis formation of crab. Furthermore, BMP 2/4 is expressed predominately in ventral and acts as a morphogen, which constructs the DV mesoderm in a concentration-dependent manner (Niehrs, 2010). In this study, BMP 2/4 was

19

ACCEPTED MANUSCRIPT up-regulated, but its antagonist (kielin/chordin-like protein-like isoform 2) was down-regulated.

As for body segmentation, gap genes are regional expression in embryo to keep continuously segmenting of larval cuticle (Mito et al., 2006). Additionally, gap genes

PT

could regulate downstream pair-rule genes and secondary pair-rule genes to maintain

RI

intrinsic polarity in each segment (Small et al., 1992). In this article, gap genes

SC

(Hunchback transcription factor, Huckebein protein), pair rule genes (Segmentation

NU

protein even-skipped, Hairy and enhancer of split, Hairy enhancer of split-like 3, Runt-related transcription factor) and secondary pair-rule genes (Protein

MA

odd-skipped, Segmentation protein Runt, Paired box protein, engrailed) were overall up-regulated, which suggests that body segmentation of E. sinensis have been

PT E

D

determined in early embryonic stages.

CE

3.7. Candidate DEGs involved in immunity

AC

Most aquatic animals undergo external fertilization, whose early embryos immersed in water full of potential pathogens. However, the developing embryos could protect themselves from pathogen attacks via a cross-generational immunity, which is depend on a series of maternal derived immune factors. This phenomenon has been discovered in crustacean Daphnia magna and Penaeus monodon (Little et al., 2003), insects Bombus terrestris (Sadd et al., 2005). This protective mechanism is achieved through maternal transcripts rather than its own transcriptional products 20

ACCEPTED MANUSCRIPT before the MZT. Similar results have been found in freshwater Hydra indicating that they can protect their embryos using maternally-encoded antimicrobial peptides (Fraune et al., 2010).

Invertebrates lack adaptive immune system and mainly depend on innate immunity

PT

to defend against microbial invasion (Kurata et al., 2006). Accordingly, many

RI

immune relevant genes have been studied separately in E. sinensis (Cui et al., 2013;

SC

Li et al., 2013). Here, 117 immune-related DEGs were detected in our datasets,

NU

including 23 predicted maternal genes, 93 predicted zygotic genes and 1 transient gene (Supplementary Table S3). These DEGs were mainly involved in pattern

MA

recognition, antioxidant system, antimicrobial peptide synthesis, phenoloxidase system, hemolymph clotting, phagocytosis, melanogenesis and cell adhesion.

D

Moreover, these immune processes could be further integrated into four signaling

PT E

pathways, including Toll, IMD, MAPK and JAK-STAT pathways according to the

CE

previous report (Li and Xiang, 2013) . Most pattern recognition receptors (PRRs), including lipopolysaccharide and β-1,

AC

3-glucan binding protein (LGBP), C-type lectin receptor protein (CLER), macrophage mannose receptor 1-like and Down Syndrome Cell Adhesion Molecule (DSCAM), were predicted maternal genes with a significant up-regulation in Fe. The PRRs are critical in the innate immune system and can initiate downstream signal transmission. For instance, lectin can identify beta 1, 3-glucan from fungi surface, LGBP can discriminate lipid polysaccharides on the surface of gram-negative bacteria 21

ACCEPTED MANUSCRIPT and virus, which could activate the prophenoloxidase system and enhance the synthesis of antimicrobial peptides (Medzhitov and Janeway, 2002). The DSCAM adhesion molecule could mediate the recognition of gram-positive and gram-negative bacteria, and activate phagocytosis (Hauton, 2012). The high expression of these

PT

PRRs in Fe showing that crab early embryos might have already had an ability to

RI

identify pathogens before the MZT.

SC

As for the downstream immune factors, most DEGs involved in antioxidant system,

NU

antimicrobial peptide synthesis, prophenoloxidase system, hemolymph clotting, phagocytosis and cell adhesion were predicted zygotic genes (Fig. 5), which suggests

MA

crab early embryos formed a relatively complete immune system after the MZT. Notably, several DEGs involved in phagocytosis were highly expressed before the

early embryos.

CE

4. Conclusion

PT E

D

MZT, suggesting that phagocytosis played important roles in immune defense of crab

Here we present the de novo assembly of the embryonic transcriptomes for E.

AC

sinensis, sequenced by the next-generation sequencing technologies. We discuss the vital events before or around the MZT in the Chinese mitten crab. A total of 566 maternally provided genes and 1165 zygotic genes were identified in the three transcriptomes. We uncover that the maternal to zygotic transition occurs at early embryogenesis of E. sinensis, and the 2-4 cell stage might be the onset of zygotic transcription. These results might extend the knowledge on the regulation mechanism 22

ACCEPTED MANUSCRIPT of the MZT, as well as contribute to the identification of candidate developmental regulators in E. sinensis, which are likely to advance our understanding of

PT

developmental mechanisms in other invertebrates.

RI

Conflict of interest disclosure

NU

SC

The authors state this research is free of conflicts of interests.

MA

Acknowledgments

We are grateful to all the laboratory members for technical advice and helpful

D

discussions. This research was supported by Natural Science Foundation of Shandong

PT E

Province [grant number ZR2017QD001]; and the Scientific and Technological Innovation Project of Qingdao National Laboratory for Marine Science and

CE

Technology [grant number 2015ASKJ02]. We also thank Shanghai Personal

AC

Biotechnology Co.,Ltd (Shanghai, China) for Illumina transcriptome sequencing.

References

Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Genome Biology 11(10): R106.

23

ACCEPTED MANUSCRIPT Anderson, K.V., Bokla, L., Nüssleinvolhard, C., 1985. Establishment of dorsal-ventral polarity in the Drosophila embryo: the induction of polarity by the Toll gene-product. Cell 42(3), 791–798.

Aoki, F., Worrad, D.M., Schultz, R.M., 1997. Regulation of transcriptional activity

PT

during the first and second cell cycles in the preimplantation mouse embryo.

RI

Developmental Biology 181(2), 296–307.

SC

Azumi, K., Sabau, S.V., Fujie, M., Usami, T., Koyanagi, R., Kawashima, T.,

NU

Fujiwara, S., Ogasawara, M., Satake, M., Nonaka, M., 2007. Gene expression

Biology 308(2), 572–582.

MA

profile during the life cycle of the urochordate Ciona intestinalis. Developmental

D

Baugh, L.R., Hill, A.A., Slonim, D.K., Brown, E.L., Hunter, C.P., 2003. Composition

PT E

and dynamics of the Caenorhabditis elegans early embryonic transcriptome. Development 130(5), 889–900.

CE

Benoit, B., He, C.H., Zhang, F., Votruba, S.M., Tadros, W., Westwood, J.T., Smibert,

AC

C.A., Lipshitz, H.D., Theurkauf, W.E., 2009. An essential role for the RNA-binding protein Smaug during the Drosophila maternal-to-zygotic transition. Development 136(6), 923–932.

Bettegowda, A., Lee, K.B., Smith, G.W., 2007. Cytoplasmic and nuclear determinants of the maternal-to-embryonic transition. Reproduction, Fertility & Development 20(1), 45–53. 24

ACCEPTED MANUSCRIPT Ciechanover, A., 1998. The ubiquitin-proteasome pathway: onprotein death and cell life. Embo Journal 17(24), 7151–7160.

Cui, Z.X., Li, X.H., Liu, Y., Song, C.W., Hui, M., Shi, G.H., Luo, D.L., Li, Y.D., 2013. Transcriptome Profiling Analysis on Whole Bodies of Microbial

PT

Challenged Eriocheir sinensis Larvae for Immune Gene Identification and SNP

RI

Development. Plos One 8(12), e81256.

SC

Davis, C.C., 1964. A study of the hatching process in aquatic invertebrates.

NU

Hydrobiologia 23(1-2), 253–266.

MA

De Renzis, S., Elemento, O., Tavazoie, S., Wieschaus, E.F., 2007. Unmasking activation of the zygotic genome using chromosomal deletions in the Drosophila

D

embryo. PLoS Biology 5(5), e117.

CE

185–195.

PT E

De Robertis, E.M., 2008. Evo-Devo: Variations on ancestral themes. Cell 132(2),

Ewen-Campen, B., Shaner, N., Panfilio, K.A., Suzuki, Y., Roth, S., Extavour, C.G.,

AC

2011. The maternal and early embryonic transcriptome of the milkweed bug Oncopeltus fasciatus. BMC Genomics 12(1), 61.

Fraune, S., Augustin, R., Anton-Erxleben, F., Wittlieb, J., Gelhaus, C., Klimovich, V.B., Samoilovich, M.P., Bosch, T.C.G., 2010. In an early branching metazoan, bacterial colonization of the embryo is controlled by maternal antimicrobial

25

ACCEPTED MANUSCRIPT peptides. Proceedings of the National Academy of Sciences 107(42), 18067–18072.

Garber, M., Grabherr, M.G., Guttman, M., Trapnell, C., 2011. Computational methods for transcriptome annotation and quantification using RNA-seq. Nature

PT

methods 8(6), 469–477.

RI

Giraldez, A.J., Mishima, Y., Rihel, J., Grocock, R.J., 2006. Zebrafish MiR-430

SC

promotes deadenylation and clearance of maternal mRNAs. Science 312(5770),

NU

75–79.

MA

Green, R.E., Krause, J., Briggs, A.W., Maricic, T., Stenzel, U., Kircher, M., Patterson, N., Li, H., 2010. A draft sequence of the Neandertal genome. Science 328(5979),

D

710.

PT E

Hauton, C., 2012. The scope of the crustacean immune system for disease control.

CE

Journal of Invertebrate Pathology 110(2), 251–260.

Hulskamp, M., Schroder, C., Pfeifle, C., Jackle, H., Tautz, D., 1989. Posterior

AC

segmentation of the Drosophila embryo in the absence of a maternal posterior organizer gene. Nature 338(6217), 629–632.

Irish, V., Lehmann, R., Akam, M., 1989. The Drosophila posterior-group gene nanos functions by repressing hunchback activity. Nature 338(6217), 646–648.

26

ACCEPTED MANUSCRIPT Kimmel, C.B., Ballard, W.W., Kimmel, S.R., Ullmann, B., Schilling, T.F., 1995. Stages of embryonic development of the zebrafish. Developmental Dynamics 203(3), 253–310.

Korzh, V., 2009. Before maternal-zygotic transition. There was morphogenetic

PT

function of nuclei. Zebrafish 6(3), 295–302.

RI

Kurata, S., Ariki, S., Kawabata, S., 2006. Recognition of pathogens and activation of

SC

immune responses in Drosophila and horseshoe crab innate immunity.

NU

Immunobiology 211(4), 237–249.

MA

Lecuyer, E., Yoshida, H., Parthasarathy, N., Alm, C., Babak, T., Cerovina, T., Hughes, T.R., Tomancak, P., Krause, H.M., 2007. Global analysis of mRNA

D

localization reveals a prominent role in organizing cellular architecture and

PT E

function. Cell 131(1), 174–187.

Lehmann, R., Nüssleinvolhard, C., 1987. Hunchback, a gene required for

CE

segmentation of an anterior and posterior region of the Drosophila embryo.

AC

Developmental Biology 119(2), 402–417.

Lehmann, R., Nüssleinvolhard, C., 1991. The Maternal gene nanos has a central role in posterior pattern-formation of the Drosophila embryo. Development 112(3), 679–691.

Li, F., Xiang, J., 2013. Signaling pathways regulating innate immune responses in shrimp. Fish & Shellfish Immunology 34(4), 973–980. 27

ACCEPTED MANUSCRIPT Li, X.H., Cui, Z.X., Liu, Y., Song, C.W., Shi, G.H., 2013. Transcriptome analysis and discovery of genes involved in immune pathways from hepatopancreas of microbial challenged mitten crab Eriocheir sinensis. Plos One 8(7), e68233.

Little, T.J., O'Connor, B., Colegrave, N., Watt, K., Read, A.F., 2003. Maternal

PT

transfer of strain-specific immunity in an invertebrate. Current Biology 13(6),

RI

489–492.

SC

Liu, M.M., Davey, J.W., Jackson, D.J., Blaxter, M.L., Davison, A., 2014. A

NU

conserved set of maternal genes? Insights from a molluscan transcriptome.

MA

International Journal of Developmental Biology 58(6-8), 501–511.

Livak, K.J., Schmittgen, T.D, 2001. Analysis of relative gene expression data using

D

real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25(4),

PT E

402-408.

Lovrich, G.A., 1997. Cannibalism in the snow crab, Chionoecetes opilio (O.

CE

Fabricius) (Brachyura Majidae), and its potential importance to recruitment.

AC

Journal of Experimental Marine Biology and Ecology 211(2), 225–245.

Marshall, S., Warburton, K., Paterson, B., Mann, D., 2005. Cannibalism in juvenile blue-swimmer crabs Portunus pelagicus (Linnaeus, 1766): effects of body size, moult stage and refuge availability. Applied Animal Behaviour Science 90(1), 65–82.

28

ACCEPTED MANUSCRIPT Medzhitov, R., Janeway, C.A., 2002. Decoding the patterns of self and nonself by the innate immune system. Science 296(5566), 298–300.

Mito, T., Okamoto, H., Shinahara, W., Shinmyo, Y., Miyawaki, K., Ohuchi, H., Noji, S., 2006. Kruppel acts as a gap gene regulating expression of hunchback and

PT

even-skipped in the intermediate germ cricket Gryllus bimaculatus.

RI

Developmental Biology 294(2), 471–481.

SC

Moksnesa, P.O., 1997. Cannibal–prey dynamics in young juveniles and postlarvae of

NU

the blue crab. Journal of Experimental Marine Biology & Ecology 215(2),

MA

157-187.

Montú, M., Anger, K., Bakker, C., 1996. Larval development of the Chinese mitten

D

crab Eriocheir sinensis H. Milne-Edwards (Decapoda: Grapsidae) reared in the

PT E

laboratory. Helgoländer Meeresuntersuchungen 50(2), 223–252.

Mortazavi, A., Williams, B.A., Mccue, K., Schaeffer, L., Wold, B., 2008. Mapping

CE

and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods 5(7),

AC

621–628.

Moussian, B., Roth, S., 2005. Dorsoventral axis formation in the Drosophila embryo—Shaping and transducing a morphogen gradient. Current Biology 15(21), R887–R899.

29

ACCEPTED MANUSCRIPT Nestorov, P., Battke, F., Levesque, M.P., Gerberding, M., 2013. The maternal transcriptome of the crustacean Parhyale hawaiensis is inherited asymmetrically to invariant cell lineages of the ectoderm and mesoderm. Plos One 8(2), e56049.

Niehrs, C., 2010. On growth and form: a Cartesian coordinate system of Wnt and

PT

BMP signaling specifies bilaterian body axes. Development 137(6), 845–857.

RI

Nüssleinvolhard, C., Frohnhofer, H.G., Lehmann, R., 1987. Determination of

SC

Anteroposterior Polarity in Drosophila. Science 238(4834), 1675–1681.

NU

Price, A.L., Modrell, M.S., Hannibal, R.L., Patel, N.H., 2010. Mesoderm and

MA

ectoderm lineages in the crustacean Parhyale hawaiensis display intra-germ layer compensation. Developmental Biology 341(1), 256–266.

D

Rosenberg, M.I., Lynch, J.A., Desplan, C., 2009. Heads and tails: evolution of

PT E

antero-posterior patterning in insects. BBA-Gene Regulatory Mechanisms

CE

1789(4), 333–342.

Sadd, B.M., Kleinlogel, Y., Schmid-Hempel, R., Schmid-Hempel, P., 2005.

AC

Trans-generational immune priming in a social insect. Biology Letters 1(4), 386–388.

Small, S., Blair, A., Levine, M., 1992. Regulation of even-skipped stripe-2 in the Drosophila embryo. Embo Journal 11(11), 4047–4057.

30

ACCEPTED MANUSCRIPT Struhl, G., 1989. Differing strategies for organizing anterior and posterior body pattern in Drosophila embryos. Nature 338(6218), 741–744.

Sultan, M., Schulz, M.H., Richard, H., Magen, A., 2008. A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome.

PT

Science 321(5981), 956–960.

RI

Tadros, W., Goldman, A.L., Babak, T., Menzies, F., Vardy, L., Orr-Weaver, T.,

SC

Hughes, T.R., Westwood, J.T., Smibert, C.A., Lipshitz, H.D., 2007. SMAUG is

NU

a major regulator of maternal mRNA destabilization in Drosophila and its translation is activated by the PAN GU kinase. Developmental cell 12(1),

MA

143–155.

D

Tadros, W., Lipshitz, H.D., 2009. The maternal-to-zygotic transition: a play in two

PT E

acts. Development 136(18), 3033–3042.

Tang, F., Barbacioru, C., Bao, S., Lee, C., Nordman, E., Wang, X., Lao, K., Surani,

CE

M.A., 2010. Tracing the derivation of embryonic stem cells from the inner cell

AC

mass by single-cell RNA-Seq analysis. Cell stem cell 6(5), 468–478.

Walser, C.B., Lipshitz, H.D., 2011. Transcript clearance during the maternal-to-zygotic transition. Current Opinion in Genetics & Development 21(4), 431–443.

Wang, C.W., Klionsky, D.J., 2003. The molecular mechanism of autophagy. Molecular Medicine 9(3-4), 65–76. 31

ACCEPTED MANUSCRIPT Wei, Z., Angerer, R.C., Angerer, L.M., 2006. A database of mRNA expression patterns for the sea urchin embryo. Developmental Biology 300(1), 476–484.

Huang, X.R., Zhuang, P., Zhang, L.Z., Qu, L., Yao, Z.F., Liu, T., Feng, G.P., Liu, J.Y., 2011. Embryonic development and the variation of some metabolism

PT

enzyme activity during embryonic development of Chinese mitten crab

RI

(Eriocheir sinensis). Journal of Fisheries of China 35(2), 192–199.

SC

Zhao, Y.L., Du, N.S., Lai, W., 1993. Effects of different gradient temperature on

NU

embryonic development of the chinese mitten-handed crab, Eriocheir

AC

CE

PT E

D

MA

sinensis(Crustacea, decapod). Zoological Research 14(1), 49–53.

32

ACCEPTED MANUSCRIPT Figure captions: Fig. 1. KEGG pathway enrichment analysis of DEGs in the E. sinensis embryo transcriptomes. A: Human diseases, B: Organismal systems, C: Cellular processes, D: Environmental information processing, E: Genetic information processing, F:

RI

PT

Metabollism.

SC

Fig. 2. RT-qPCR validation of the selected DEGs in Illumina sequencing. Each

NU

column represented the mean of triplicate assays within ± S.D. Significant differences

MA

of the same gene in different developmental stages (Fe, Cs and Bs) were indicated by

PT E

D

different letters, and the same letters indicate no significant difference.

Fig. 3. The putative degradation process of maternal products in E. sinensis embryo, which was constructed based on KEGG mapping. The predicted maternal genes were

CE

labeled with green, zygotic genes were labeled in red, transient genes were labeled

AC

with blue, genes with no expression differences were labeled in pink, respectively.

Fig. 4. The putative axis formation and segmentation in E. sinensis embryo. The predicted maternal genes were labeled with green, zygotic genes were labeled in red,

33

ACCEPTED MANUSCRIPT transient gene were labeled with blue, genes with no expression differences were labeled in pink, respectively.

PT

Fig. 5. Outline of the innate immune system in E. sinensis. It was constructed based on the immune system of shrimp (Li and Xiang, 2013). The predicted maternal genes

RI

were labeled with green, the predicted zygotic genes were labeled in red, the predicted

SC

transient gene were labeled with blue, respectively, and the pink color represents

AC

CE

PT E

D

MA

NU

genes with no expression differences in three stages.

34

ACCEPTED MANUSCRIPT

Table 1

Details of primers used for RT-qPCR Cluster

Gene ID

Primer

Sequence (5'-3')

ANK

EM

c99630_g3_i4

RT-ANK-F

GCTCCAGATCAACTCCGAAT

RT-ANK-R

GGATGTCACAGCAGACCACT

RT-CCNB-F

CCTTACACACTCCACAGACC

RT-CCNB-R

CATCAAACTTCCAACTTCCT

RT-BCDC-F

TTGGTTCCTTTAGTCCGTT

RT-BCDC-R

CTCTCAGTTCTTTGGTTCAC

CCNA

SM

M

c102500_g1_i3

c99718_g5_i1

RT-CCNA-F

FA11

CP

SZ

Z

MZ

CGTATGTATGGTGTAAGGGA

RT-DICER1-R

D

HB

LZ

c88186_g1_i2

c93257_g3_i2

c90946_g1_i2

c99607_g2_i2

c104244_g1_i2

TGCTAACTTGGATGCGAGAA

RT-DICER1-F

GGGGTGAAGATGGTGGAGTG

RT-RPBII-F

TTCGCCTTCTCCGCATCCCT

RT-RPBII-R

TCCCCCTTATCTTCCCCGCC

RT-SODC-F

CGGTGTGGGGTGTGAGGGAG

RT-SODC-R

CGCAGGTTTGTGGGAGGAGA

RT-HB-F

TCCCTACGCGTCGGACTTGC

RT-HB-R

CACTGGCACTGTCCTTCCTG

RT-FA11-F

GCTCTATTGGGTCCCTTTGG

RT-FA11-R

CTTCCGTCTGCTTCTTGCTG

RT-CP-F

GGGAGCCGACACCCAAGAAG

RT-CP-R

TGACGCAAGACCACAACAAT

PT E

SODC

EZ

c94391_g1_i1

CE

RPBII

EZ

AC

DICER1

ACATAACAGTGGTGAGGGCC

MA

RT-CCNA-R

RI

c103426_g1_i1

SC

BCDC

LM

NU

CCNB

PT

Gene

35

ACCEPTED MANUSCRIPT c102162_g3_i1

RT-NANOS-F

GGGGCAATATTTGAGGGTGT

RT-NANOS-R

AACTGTTCGGCAGGAGAGGA

CE

PT E

D

MA

NU

SC

RI

PT

T

AC

NANOS

36

ACCEPTED MANUSCRIPT

Table 2

Fe

Cs

Bs

Raw reads

57,135,691

36,954,342

32,813,420

Q20 percentage %

92.02

92.94

92.98

GC percentage %

52.12

52.83

53.42

Clean reads

46,892,436

31,071,869

Total clean base pairs (Gb)

11.62

6.84

628,264

MA

Total number of contigs Mean length of contigs (bp)

287.92

N50 of contigs (bp)

344

PT E

D

Total number of transcripts

270,813

Mean length of transcripts (bp)

601

N50 of transcripts (bp)

859 32,088

Mean unigene length (bp)

1,320

AC

CE

Total number of annotated unigenes

N50 of unigenes (bp) Database

2,389 Number

Percentage

GO

25,558

79.65%

KEGG

7,158

22.31%

eggNOG

29,461

91.81%

KO

6,688

20.84% 37

27,633,559

SC

NU

Assembly statistics

PT

Statistics

RI

Statistics of output sequencing in the E. sinensis embryo transcriptomes

7.68

ACCEPTED MANUSCRIPT GO, Gene Onotology; eggnog, evolutionary genealogy of genes: Non-supervised Orthologous

AC

CE

PT E

D

MA

NU

SC

RI

PT

Groups; KEGG, Kyoto Encyclopedia of Genes and Genomes; KO, KEGG Orthology.

38

ACCEPTED MANUSCRIPT Table 3

The conserved maternal and zygotic unigenes in E. sinensis and other invertebrates COMATs (%)

Zygotic-only

COZYTs (%)

Source

E. sinensis

566

103 (18.20%)

1165

266 (22.83%)

This study

C. elegans

953

76 (7.97%)

441

46 (10.43%)

Baugh et al., 2003

D. melanogaster

563

106 (18.83%)

300

42 (14.0%)

De Renzis et al., 2007

C. intestinalis

4041

126 (3.12%)

1646

P. hawaiensis

690

117 (16.95%)

568

L. stagnalis

481

175 (36.38%)

--

RI

PT

Maternal-only

Azumi et al., 2007

74 (13.03%)

Nestorov et al., 2013

SC

42 (2.55%)

--

Liu et al., 2014

NU

Species

MA

Maternal-only transcripts present in the oocyte or egg but not in developing embryos, zygotic-only transcripts generate after zygotic transcriptome activation; COMATs represent

D

conserved maternal transcripts, CONZYTs represent conserved zygotic transcripts; -- mean no

AC

CE

PT E

data.

39

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

PT

Fig. 1

40

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

PT

Fig.2

41

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

PT

Fig.3

42

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

PT

Fig.4

43

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

PT

Fig.5

44

ACCEPTED MANUSCRIPT Abbreviations: RT-qPCR, quantitative real-time PCR; DEG, differentially expressed gene; MZT, maternal to zygotic transition; ZGA, zygotic genome activation; COMATs, conserved maternal transcripts; COZYTs, conserved zygotic transcripts; Fe, fertilized

AC

CE

PT E

D

MA

NU

SC

RI

PT

egg; Cs, cleavage stage; Bs, blastula stage.

45

ACCEPTED MANUSCRIPT Highlights 1. It is the first crab transcriptome profile to investigate the maternal-to-zygotic transition.

PT

2. A total of 32,088 annotated unigenes were obtained through de novo

RI

transcriptome assembly.

SC

3. We unravel the MZT in E. sinensis occurs around cleavage stage (2-4 cell stage).

NU

4. The DEGs were participated in the MZT, morphogenesis and immunity,

AC

CE

PT E

D

MA

respectively.

46