Journal of Oceanology and Limnology   2019, Vol. 37 issue(1): 223-234     PDF       
http://dx.doi.org/10.1007/s00343-019-7254-6
Institute of Oceanology, Chinese Academy of Sciences
0

Article Information

SONG Chengwen, LIU Lei, HUI Min, LIU Yuan, LIU Hourong, CUI Zhaoxia
Primary molecular basis of androgenic gland endocrine sex regulation revealed by transcriptome analysis in Eriocheir sinensis
Journal of Oceanology and Limnology, 37(1): 223-234
http://dx.doi.org/10.1007/s00343-019-7254-6

Article History

Received Sep. 8, 2017
accepted in principle Oct. 15, 2017
accepted for publication Jan. 10, 2018
Primary molecular basis of androgenic gland endocrine sex regulation revealed by transcriptome analysis in Eriocheir sinensis
SONG Chengwen1,2, LIU Lei1, HUI Min1, LIU Yuan1,2, LIU Hourong1,3, CUI Zhaoxia1,2     
1 CAS Key Laboratory of Experimental Marine Biology, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China;
2 Laboratory for Marine Biology and Biotechnology, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China;
3 University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: In crustaceans, the male sexual differentiation and maintenance are specially regulated by androgenic gland (AG). However, little is known about the genes involved in the regulation process. RNA-Seq was performed on AG with ejaculatory duct (AG_ED) and ejaculatory duct (ED) as control in Eriocheir sinensis, one of the most important economic and fishery crabs with typically sex dimorphism. A total of 925 unigenes were identified as differentially expressed genes (DEGs) and the expression of nine genes randomly selected was confirmed by qRT-PCR. 667 unigenes were up-regulated in AG_ED, being supposed to be AG preferential genes. Among them, the full length of insulin-like androgenic gland factor (IAG) cDNA named as Es-IAG was obtained as a logo gene of AG, which together with the genes insulin-like receptor (INR), and single insulin binding domain protein (SIBD), might constitute the sex regulation pathway. Several sex related genes were identified, and their function will have to be investigated. Also, the identification of juvenile hormone epoxide hydrolase 1 (JHEH1), ecdysteroid 22-hydroxylase (DIB) and ecdysone receptor (ECR) preliminarily clarified the molecular regulation mechanism of eyestalk-AG-testis axis, which plays important roles in molting and reproduction. The results will enhance our understanding for the molecular basis of the AG involved in male sex regulation in crabs.
Keywords: androgenic gland    comparative transcriptome    Eriocheir sinensis    insulin-like androgenic gland factor (IAG)    sex regulation    
1 INTRODUCTION

The androgenic gland (AG) has been specially found in crustaceans as an endocrine organ involving in the regulation of the male sexual characteristics differentiation by secretion of hormones since decades ago (Charniaux-Cotton, 1955; Katakura, 1989). In 2002, Khalaila et al. proposed an eyestalk-AG-testis endocrine axis for male sex differentiation based on the relationship of AG, sinus gland and male reproductive system in decapod crustaceans. Although the followed studies have been focused on the morphology, anatomy and physiology of AG (Guan, 2003; Manor et al., 2004; Aflalo et al., 2006; Tropea et al., 2011), the embedded mechanism of sex regulation by AG is still vague.

At molecular level, the insulin-like androgenic gland hormone (IAG) gene as a logo gene of AG has been identified in many decapod crustaceans, such as in crabs Callinectes sapidus, Portunus pelagicus and Scylla paramamosain (Sroyraya et al., 2010; Chung et al., 2011; Huang et al., 2014). Silencing IAG expression in some species can lead to demasculinization (Ventura et al., 2009; Rosen et al., 2010), and even result in allmale monosex populations that are the progeny of fully sexually reversed males (neo-females) (Ventura and Sagi, 2012). However, more novel genes and pathways potentially involved in the sex-determination process of IAG remain to be characterized, which are necessary to uncover the molecular mechanism of sex differentiation and benefit the application of biotechnological manipulation for economic crustaceans.

The decapod crustaceans include many commercially valuable lobster, crayfish, shrimp and crab species. The Chinese mitten crab Eriocheir sinensis is one of the important aquaculture crabs in China. The production of E. sinensis has reached 8 million tons per year recorded by China Fisheries Yearbook. The culture of male population can make the cultured species reach market size at a faster rate, while female crabs are of higher economic values than males. At the same time, the culture of monosex population can minimize the invasion of E. sinensis in Europe and North America (Herborg et al., 2003; Dittel and Epifanio, 2009). Therefore, it is important to gain a better understanding of AG regulated sex differentiation in this species.

Transcriptome analysis provides the foundation for studying gene expression and pathway prediction (Zhang et al., 2014; Yang et al., 2015). The transcriptome analyses for AG have been performed in the oriental river prawn Macrobrachium nipponense (Jin et al., 2013) and the Eastern spiny lobster Sagmariasus verreauxi (Chandler et al., 2015), and many candidate novel sex-related genes have been identified. For example, 47 sex-determination related gene families were identified from the M. nipponense androgenic gland transcriptome and 40 candidate novel genes were found, which are extremely high expressed in the androgenic compared to other sex glands (Jin et al., 2013). However, the present study is the first instance to identify genes exclusively expressed in the androgenic gland and to predict their function and regulation according to the eyestalkAG-testis endocrine axis by comparative transcriptome analysis in crab. Genes with function of sex determination, reproduction and endocrine regulation were identified and their putative interactions were deduced accordingly. The results are expected to enhance our understanding for the molecular basis of AG- male sex regulation in crabs.

2 MATERIAL AND METHOD 2.1 Experimental design and animal sampling

A total of 60 healthy mature males E. sinensis at intermolt stage (mean weight=50 g) were obtained from a commercial market in Qingdao. Individuals were immediately brought to laboratory and acclimated in aquaria with aerated water at room temperature for 72 h before sampling (Zhang et al., 2015, 2017). Under a microscope, a quite small pair of AGs was observed on the surface of ED. Considering the simplicity of the experiment and to obtain enough RNA-Seq samples easily, two kinds of tissues were dissected: (1) androgenic gland with ejaculatory duct (AG_ED), and (2) the ejaculatory duct (ED) without AG. After dissection, the samples were immediately frozen in liquid nitrogen for RNA extraction. Finally, AG_EDs from 30 individuals were mixed for one sequencing sample, while EDs from the other 30 individuals were mixed for the other one. By comparative transcriptomic analysis between AG_ED and ED, the genes preferentially expressed in AG were expected to be identified.

2.2 Transcriptome library construction and Illumina sequencing

The total RNA was extracted using the TRIzol kit (Invitrogen, USA) according to the manufacturer's protocol. RNA concentration and integrity were measured using Qubit (Life Technologies) and the RNA Nano 6000 (Agilent Technologies), respectively. A total amount of 3 μg RNA per sample was used for preparation. Sequencing libraries were generated using NEBNext RNA Library Prep Kit for Illumina (NEB) according to the manufacturer's protocol (http://international.neb.com/protocols/2014/12/01/protocol-for-use-with-nebnext-poly-a-mrna-magneticisolation-module-neb-e7490) and index codes were added to assign sequences to each sample. After that, the libraries were sequenced on an Illumina Hiseq 2500 platform and paired-end reads were generated.

2.3 Transcriptome assembly, gene annotation and classification

All raw sequences were deposited to NCBI Short Read Archive (SRA) database (http://www.ncbi.nlm.nih.gov/Traces/sra/). After removing adaptor sequences, ambiguous 'N' nucleotides and low quality sequences, the remaining clean reads were assembled using Trinity software (Grabherr et al., 2011). Gene annotation was analyzed using an established approach as described by Cui et al. (2013). All the unigenes were compared with the NCBI (National Center for Biotechnology Information) Nr, Nt and Swiss-prot (http://www.ebi.ac.uk/uniprot/) databases for gene annotation. Gene encoding protein domains were identified by searching against Protein Family (Pfam) database (http://pfam.janelia.org/). After that, GO (Gene Ontology; http://www.geneontology.org/) and KOG (euKaryotic Ortholog Group; http://www.ncbi.nlm.nih.gov/COG/) were used to obtain the potential gene functions. GO enrichment analysis was performed using the GO seq R packages (Young et al., 2010). Using Enzyme Commission (EC) terms, the potential metabolic pathway for the unigenes was constructed from KEGG (Kyoto Encyclopedia of Genes and Genomes; http://www.genome.jp/kegg/). Based on the annotation and published literatures (Liu et al., 2015), the genes related to sex differentiation or reproduction were manually checked.

2.4 AG preferential gene identification and enrichment

The RSEM software was used to estimate gene expression level (Li and Dewey, 2011) and then FPKM (expected number of Fragments Per Kilobase of transcript sequence per Million base pairs sequences) of each gene was calculated accordingly (Trapnell et al., 2010). Prior to DEGs analysis, the read counts were adjusted by edge R program package through one scaling normalized factor. DEGs were then identified using the DEGseq (Anders and Huber, 2010) R package. The P-values were adjusted (Storey et al., 2003) and the corrected q value < 0.005 & |log2(fold change)| > 1 (fold change > 2) was set as the threshold for significantly differential expression. From all DEGs, the up-regulated genes in AG_ED were selected and supposed to be AG preferential genes. GO enrichment and KEGG pathway analysis were also implemented for these AG preferential genes as above.

2.5 Quantitative reverse-transcription PCR (qRTPCR) confirmation

The qRT-PCR was performed used new AG_ED and ED samples to verify the accuracy of RNA-Seq and the expression level of key DEGs. Total RNA was extracted and the first-strand cDNA was synthesized by using PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Japan). Then, the cDNA was diluted 100 times by nucleic acid free water. SYBR Green Real-time RT-PCR was performed in an ABI 7500 Sequence Detection System (Applied Biosystems). Nine pairs of gene-specific primers (Table 1) were used to amplify the partial cDNA sequences, respectively. The β-actin gene was amplified in parallel as an internal control. There were three biological and three technical replicates respectively. The 2-ΔΔCt method was used to analyze the expression level of different genes (Livak and Schmittgen, 2001). Simply, the fold change in target gene relative to the β-actin endogenous control is determined by: Fold change= 2-ΔΔCt, where the Ct value is defined as the cycle number at which the fluorescence signal crosses the threshold; ΔΔCt=(CtTargetCtActin)AG_ED– (CtTargetCtActin)ED.

Table 1 Genes and primers used in the quantitative realtime PCR analysis
2.6 IAG sequence analysis and its regulation mechanism prediction

Based on Blastx algorithm, one sequence with annotation of insulin-like androgenic gland factor [Callinectes sapidus] (E-value=2.97E-22) from DEGs was found and named Es-IAG. The open reading frame (ORF) of this sequence was identified using the translate tool (http://web.expasy.org/translate/) and BLAST programs in NCBI (http://www.ncbi.nlm.nih.gov/BLAST). The structure characteristics of the deduced amino acids were analyzed by comparing with IAGs in other crabs (Sroyraya et al., 2010; Chung et al., 2011; Huang et al., 2014). The molecular weight (MW) and isoelectric point (pI) of deduced amino acid sequences were analyzed using the Compute pI/Mw tool (http://web.expasy.org/compute_pi/). The sequence identity between Es-IAG and other decapod IAG genes were calculated based on the aligned amino acid sequences using BioEdit 7.0.

Considering that AG is an endocrine organ and IAG is a logo gene of AG, the IAG, endocrine and sex related genes were screened from AG preferential genes based on KEGG pathway analysis.

3 RESULT AND DISCUSSION 3.1 Gene annotation

All raw reads were deposited into the SRA database under the accession numbers SRR2170964 (AG_ED) and SRR2170970 (ED). A total of 71 490 unigenes were assembled, among which 21 139 genes (29.56%) could be annotated in at least one database and 14 817 (20.7%) matched in Nr database (E-value < 1E-05) (Table 2). The large proportion of genes without annotation might be attributed to limited number of crustacean sequences in public databases. The top three crustacean species with the most number of sequence hits were Daphnia pulex (1 053, 7.11%), Penaeus monodon (213, 1.44%) and S. paramamosain (130, 0.88%).

Table 2 Genes annotation information for AG_ED and ED of E. sinensis

By GO analysis, a total of 16 831 unigenes (23.54%) were assigned to the biological process, cellular component and molecular function categories (Supplementary Fig. 1). By searching against KOG database, 8 677 (12.13%) unigenes were classified into 25 groups (Supplementary Fig. 2). According to KEGG analysis, 6 785 (9.5%) unigenes were assigned to 32 groups (Fig. 1). These genes were involved in more than 250 pathways, among which many pathways were related to signal transduction or endocrine regulation, such as 'GnRH signaling pathway', 'Insulin signaling pathway', 'Insulin secretion', 'Oxytocin signaling pathway' and 'Estrogen signaling pathway'.

Fig.1 Functional distribution of all unigenes in AG and ED transcriptomes of E. sinensis based on KEGG analysis
3.2 Sex determination and reproduction related genes

The AG and ED are both male-specific reproductive organs. From the assembled transcriptome, a series of crucial sex determination or germ line genes were identified (Supplementary Table 1). These genes were divided into two categories of male-related and female-related.

3.2.1 Reproduction

A number of Cytochrome P450 genes, Vitellogenin (Vg) and Vitellogenin receptor (VgR) were also identified from the transcriptome, which had been reported to show sex-biased expression and regulation for different metabolic processes in our previous study (Liu et al., 2015). As reported by Ventura et al. (2017), several CYP450s participate in the biosynthesis of the active form of the molt hormone, 20-hydroxyecdysone (20HE), including CYP302A1 and CYP18A1, which were high expressed in AG and ED in E. sinensis. However, it is the first time to reveal Vg, together with its receptor, expression in male reproductive organ.

3.2.2 Sex determination

Many of them are key switch genes in the sex determination pathways of animals (Hansen and Pilgrim, 1999; Smith et al., 2003; Kopp, 2012), such as members of Doublesex- and mab-3-related transcription factor (DMRT), Feminization (Fem), SOX, Forkhead box protein (FOXO) and Transformer (Tra). However, each of these genes will have to be tested on their function in further research.

3.3 AG preferential genes isolation and enrichment

The focus of our study is to mine AG preferential genes and deeply identify possible genes and pathways involved in the sex regulation function of AG. By DEG analysis, 925 unigenes were revealed to be significantly differentially expressed between AG_ ED and ED. Of which, 667 unigenes were supposed to be AG preferential genes with up-regulated expression in AG_ED. All nine tested genes including IAG, INR, SIBD and so on, showed a consistent expression between qRT-PCR and RNA-Seq, which verified the reliability of RNA-Seq (Fig. 2).

Fig.2 Comparison of gene expression results revealed by qRT-PCR and RNA-Seq "Fold change" means the value of gene expression of AG_ED: ED.

GO enriched analysis showed that the most representative function of AG preferential genes was 'single organism metabolic process' (142 unigenes, 32.87% of the genes with GO terms), followed by 'oxidation-reduction process' (77, 17.82%) and 'oxidoreductase activity' (70, 16.2%) (Supplementary Fig. 3). By KEGG pathway analysis, these genes were assigned to 180 different pathways, of which the genes involved in 'Carbon metabolism, Biosynthesis of amino acids', 'Oxidative phosphorylation', 'Glycolysis/Gluconeogenesis' and 'Citrate cycle' were the most abundant. Considering the specific sex regulation function of AG, endocrine and sex regulation related genes from AG preferential genes and their involving pathways were further identified.

3.4 Endocrine sex regulation related genes

Since the silence of IAG has been reported to inhibit development of male secondary sex characteristics and even result in all male population in M. rosenbergii (Ventura and Sagi, 2012), the putative endocrine sex regulation of IAG was established accordingly in this study (Table 3).

Table 3 Putative endocrine sex regulation related genes identified from AG preferential genes in E. sinensis

Insulin-like hormone related genes In all AG preferential genes, IAG showed the highest expression level with an FPKM value 15 821.59 in AG_ED, while the value was 3.12 in ED with a fold change 5 070 (Table 1). This result was confirmed by qRTPCR and obviously highlighted the biosynthesis of IAG in AG. The full length of Es-IAG cDNA was 1 455 bp and the ORF was 456 bp encoding 151 amino acids (Fig. 3a). The molecular weight of the predicted Es-IAG protein was 16.59 kDa with a theoretical isoelectric point of 5.30. When compared with other decapod crustaceans (Manor et al., 2007; Sroyraya et al., 2010; Chung et al., 2011; Mareddy et al., 2011; Ma et al., 2013; Huang et al., 2014), Es-IAG was also composed of a signal peptide, a B chain, a C peptide and an A chain which might insure the function performing of the proteins. However, its amino acid sequences showed low sequence similarities with other species from 11.7% (with Cherax quadricarinatus) to 43.3% (with C. sapidus), and the low similarity of IAG is also found in other crabs or shrimps (Fig. 3b).

Fig.3 Sequence analysis of Es-IAG a. nucleotide and its deduced structure. Nucleotide and deduced amino acid residues are numbered on the right. The start and stop codons are boxed and an ATTAAA polyadenylation signal in the 3′-UTR is underlined. The predicted cleavage sites (RxxR, RFRR) are shaded in gray. Six cysteine residues are highlighted in bold letters. N-glycosylation sites with the sequence of NxS/T are highlighted with blue circles. Asterisk denotes the stop codon; b. sequence identity matrix of decapod IAG genes based on aligned amino acid sequences. The species and the GenBank accession numbers are as follow: Callinectes sapidus IAG1(CsIAG1, AEI72263), C. sapidus IAG2 (CsIAG2, AHM93481), Cherax destructor (CdIAG, ACD91988), C. quadricarinatus (CqIAG, ABH07705), Fenneropenaeus chinensis IAG1(FcIAG1, AFU60548), F. chinensis IAG2 (FcIAG2, AFU60549), Jasus edwardsii (JeIAG, AIM55892), Litopenaeus vannamei (LvIAG, AIR09497), Macrobrachium lar (MlIAG, BAJ78349), M. nipponense IAG1 (MnIAG1, AGB56976), M. nipponense IAG2 (MnIAG2, AHA33389), M. rosenbergii (MrIAG, ACJ38227), M. vollenhovenii (MvIAG, AHZ34725), Marsupenaeus japonicus (MjIAG, BAK20460), Palaemon pacificus (PpacIAG, BAJ84109), P. paucidens (PpauIAG, BAJ84108), Penaeus monodon (PmIAG, ADA67878), Portunus pelagicus (PpeIAG, ADK46885), Sagmariasus verreauxi (SvIAG, AHY99679) and Scylla paramamosain (SpIAG, AIF30295).

Two insulin related genes besides Es-IAG, insulinlike receptor (INR) and single insulin binding domain proteins (SIBD) were also up-regulated in AG_ED with high expression level. In C. quadricarinatus, an insulin-like growth factor-binding protein (IGFBP), is identified and found to bind the gender-specific androgenic hormone IAG (Rosen et al., 2013). The binding complex then moves in the circulatory system and subsequently regulate the signaling by acting as competitors to the designated receptors, leading to the release of insulin-like growth factor (IGF) at the target site (Hwa et al., 1999; Kelley et al., 2002). As reported, EsSIBD is a potential member of the IGFBP family with a close relationship with IGFBP7 and involved in the endocrine system (Gai et al., 2010). Therefore, upon Es-IAG with SIBD binding, INR acting as a tyrosine kinase might initiate a phosphorylation cascade, and further regulation is maintained through triggering specific signaling pathway which promotes the sex regulation process.

Endocrine related genes A total of 24 genes involved in endocrine system were identified from AG preferential genes (Table 3). By KEGG analysis, we found that nine of them (LDH, PDHA, PDHB, PYG, PK, PCK, PGAM, CALM, CREB2) were involved in Glucagon signaling pathway, while five genes (PYG, HK, PCK, CALM, MKNK) were in Insulin signaling pathway (Supplementary Fig. 4) and four genes (PLD1_2, CALM, PRKCD, CREB2) were in GnRH signaling pathway (Supplementary Fig. 5). In insulin signaling pathway, insulin binding to its receptor results in the tyrosine phosphorylation of the insulin receptor substrates (IRS) by the insulin receptor tyrosine kinase (INSR). As a member of insulin family, the role of IAG might similar with insulin and triggered the signaling pathway. The main function of this pathway is maintaining glucose homeostasis, and it also involved in proliferation and differentiation.

It is interesting that GnRH signaling pathway can regulate the gonadotrophins gene expression and secretion, and thereby mediate control of the reproduction in many vertebrates (Bliss et al., 2010). The presence of GnRH-like receptors in the oriental river prawn M. nipponense may suggest the function of GnRH signaling pathway in regulating ovarian development in the crustacean (Du et al., 2015). In the male blue swimming crab, P. pelagicus, lGnRH-IIIlike peptide possibly enhances germ cell proliferation and maturation in the testes and increases sperm production (Senarai et al., 2016). At same time, Siangcham et al. (2013) found that gonadotropinreleasing hormones can induce increases in AG and testicular size, IAG production and spermatogenesis in the giant freshwater prawn, M. rosenbergii. Calmodulin can promote the exocytotic gonadotropin secretion as well as the expression of many genes including those for gonadotropin subunits (Armstrong et al., 2009). JDP2 is a member of the AP-1 family which is partially regulated by c-Jun, while c-Jun is the downstream gene in GnRH signaling pathway. Transcript factor CREB2 is also involved in this pathway (Supplementary Fig. 5). Combined with the above researches and our results, IAG might play important roles in the sex regulation and it may interact with GnRH signaling pathway.

Sex related genes Seven genes were found to be associated with sex differentiation or reproduction. Two of them (Spermatogenesis-associated protein 1, SPATA1; Spermatogenesis-associated protein 13, SPATA13) are relevant to spermatogenesis (Qiao et al., 2015). Sox17 is also involved in spermatogenesis in the frog, mouse and scallop (Kanai et al., 1996; Hudson et al., 1997; Liang et al., 2017). LIM (Male reproductive-related LIM protein) can encode protein related to male reproduction (He et al., 2012), while BOKA (Bcl-2-related ovarian killer protein homolog A) plays a role in gonad development (Hsu et al., 1997). Most of these genes are highly expressed in testis and are candidate reproduction-related genes. So, these genes might be the downstream target genes of the IAG signaling pathway.

3.5 The sex regulation of eyestalk for androgenic gland

Other endocrine regulation hormone related genes such as juvenile hormone epoxide hydrolase 1 (JHEH1), ecdysteroid 22-hydroxylase (DIB), and ecdysone receptor (ECR), were also found upregulated in AG_ED. JHEH1 is well known to switch on the degradation of juvenile hormone (JH), which is synthesized primarily by the corpora allata in insect. In crustaceans, the methyl farnesoate (MF) secreted by mandibular organs (MO), can affect endocrine regulation and reproduction in a manner similar to JH (Olmstead and Leblanc, 2002). Previous studies have proposed the effect of MF on AG and testis in crustaceans (Laufer et al., 1993; Nagaraju, 2007). In this study, the high expression of JHEH1 in AG_ED probably provides a genetic proof for MF regulatory role on AG (Fig. 4).

Fig.4 Diagrammatic representation of the putative molecular basis by eyestalk-androgenic gland-testis endocrine axis for sex regulation in E. sinensis The genes of IAG, JHEH, ECR and DIB are identified from AG preferential genes in the present study. XO-SG: X-organ/sinus gland complex; MOIH: mandibular organ-inhibiting hormone; MIH: molt-inhibiting hormone; MF: methyl farnesoate; EC: ecdysone; +: indicates positive influence; -: indicates negative regulation.

Meanwhile, DIB is related to ecdysone (EC) biosynthesis and ECR is compulsory for the ecdysteroid action (Durica et al., 1999; Hopkins, 2009; Nakagawa and Henrich, 2009), whose expressions have been reported to be molt cycledependent in many organs, such as ovary and testis (Chung et al., 1998; Durica et al., 2002; Kim et al., 2005; Asazuma et al., 2007). The found of ECR imply the recognition of a molting hormone (ecdysteroid) in the AG. In eyestalk-ablated male Litopenaeus vannamei, the Lv-IAG expression is up-regulated and synchronized to the molt cycle (Vázquez-Islas et al., 2015). These results support a reproductive endocrine axis as an eyestalk-androgenic gland linkage by ecdysteroid. In common, both the synthesis of MF and EC are regulated by X-organ/sinus gland (XOSG) of eyestalk through secretions of mandibular organ-inhibiting hormone (MOIH) and moltinhibiting hormone (MIH) (Nagaraju, 2007; Chang and Mykles, 2011). As reported by Rotllant et al. (2017), the AG itself is regulated by the XO-SG through secreting a multitude of neurohormones including gonad/vitellogenesis inhibiting hormone (GIH/VIH) and MIH, which are believed to regulate IAG expression directly and have an ongoing inhibitory effect during sexual development. Based on these results, we further find the gene information to explain the regulation of eyestalk for AG (Fig. 4). Together, by regulation of eyestalk-borne hormones for AG, the latter gland controls the differentiation and functioning of male sexual characteristics through IAG pathway, which might partly contribute to explain the endocrine axis-like relationship between eyestalk, AG, and testis at molecular level.

4 CONCLUSION

Our study represents the first AG related transcriptome analysis in crab. It successfully isolates dozens of sex differentiation or reproduction related genes, such as IAG with full length and other male- or female-related genes. A series of AG preferential genes are identified, including IAG, endocrine and sex related genes. The primary molecular basis of the eyestalk-androgenic gland-testis is revealed, which advance our understanding of AG involving endocrine sex regulation pathways in crustaceans. The results will facilitate further research on the mechanisms of sex determination in crabs.

5 DATA AVAILABILITY STATEMENT

Authors can confirm that all relevant data are included in the article and its supporting information files. All raw sequences were deposited to NCBI Short Read Archive (SRA) database (http://www.ncbi.nlm.nih.gov/Traces/sra/) under accession numbers SRR2170964 (AG_ED) and SRR2170970 (ED).

Electronic supplementary material

Supplementary material (Supplementary Figs.1–5 and Supplementary Table 1) is available in the online version of this article at https://doi.org/10.1007/s00343-019-7254-6.

References
Aflalo E D, Hoang T T T, Nguyen V H, Lam Q, Nguyen D M, Trinh Q S, Raviv S, Sagi A. 2006. A novel two-step procedure for mass production of all-male populations of the giant freshwater prawn Macrobrachium rosenbergii. Aquaculture, 256(1-4): 468-478. DOI:10.1016/j.aquaculture.2006.01.035
Anders S, Huber W. 2010. Differential expression analysis for sequence count data. Genome Biology, 11(10): R106. DOI:10.1186/gb-2010-11-10-r106
Armstrong S P, Caunt C J, Fowkes R C, Tsaneva-Atanasova K, McArdle C A. 2009. Pulsatile and sustained gonadotropinreleasing hormone (GnRH) receptor signaling:does the Ca2+/NFAT signaling pathway decode GnRH pulse frequency?. Journal of Biological Chemistry, 284(51): 35746-35757. DOI:10.1074/jbc.M109.063917
Asazuma H, Nagata S, Kono M, Nagasawa H. 2007. Molecular cloning and expression analysis of ecdysone receptor and retinoid X receptor from the kuruma prawn, Marsupenaeus japonicus. Comparative Biochemistry and Physiology Part B:Biochemistry and Molecular Biology, 148(2): 139-150.
Bliss S P, Navratil A M, Xie J J, Roberson M S. 2010. GnRH signaling, the gonadotrope and endocrine control of fertility. Frontiers in Neuroendocrinology, 31(3): 322-340. DOI:10.1016/j.yfrne.2010.04.002
Chandler J C, Aizen J, Elizur A, Battaglene S C, Ventura T. 2015. Male sexual development and the androgenic gland:novel insights through the de novo assembled transcriptome of the Eastern spiny lobster, Sagmariasus verreauxi. Sexual Development, 9(6): 338-354. DOI:10.1159/000443943
Chang E S, Mykles D L. 2011. Regulation of crustacean molting:a review and our perspectives. General and Comparative Endocrinology, 172(3): 323-330.
Charniaux-Cotton H. 1955. Le déterminisme hormonal des caractères sexuels d'Orchestia gammarella (Crustacé Amphipode). Comptes Rendus de l'Academie des Sciences, 240: 1487-1489.
Chung A C K, Durica D S, Hopkins P M. 1998. Tissue-specific patterns and steady-state concentrations of ecdysteroid receptor and retinoid-X-receptor mRNA during the molt cycle of the fiddler crab, Uca pugilator. General and Comparative Endocrinology, 109(3): 375-389. DOI:10.1006/gcen.1997.7046
Chung J S, Manor R, Sagi A. 2011. Cloning of an insulin-like androgenic gland factor (IAG) from the blue crab, Callinectes sapidus:implications for eyestalk regulation of IAG expression. General and Comparative Endocrinology, 173(1): 4-10.
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 challenged Eriocheir sinensis larvae for immune gene identification and SNP development. PLoS One, 8(12): e82156. DOI:10.1371/journal.pone.0082156
Dittel A I, Epifanio C E. 2009. Invasion biology of the Chinese mitten crab Eriochier sinensis:a brief review. Journal of Experimental Marine Biology and Ecology, 374(2): 79-92. DOI:10.1016/j.jembe.2009.04.012
Du Y X, Ma K Y, Qiu G F. 2015. Discovery of the genes in putative GnRH signaling pathway with focus on characterization of GnRH-like receptor transcripts in the brain and ovary of the oriental river prawn Macrobrachium nipponense. Aquaculture, 442: 1-11. DOI:10.1016/j.aquaculture.2015.02.016
Durica D S, Chung A C K, Hopkins P M. 1999. Characterization of EcR and RXR gene homologs and receptor expression during the molt cycle in the crab, Uca pugilator. Integrative and Comparative Biology, 39(4): 758-773.
Durica D S, Wu X H, Anilkumar G, Hopkins P M, Chung A C K. 2002. Characterization of crab EcR and RXR homologs and expression during limb regeneration and oocyte maturation. Molecular and Cellular Endocrinology, 189(1-2): 59-76. DOI:10.1016/S0303-7207(01)00740-7
Gai Y C, Wang L L, Song L S, Zhao J M, Qiu L M, Xing K Z. 2010. A putative endocrine factor SIBD (single insulin binding domain protein) involved in immune response of Chinese mitten crab Eriocheir sinensis. Fish & Shellfish Immunology, 28(1): 10-17.
Grabherr M G, Haas B J, Yassour M, Levin J Z, Thompson D A, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q D, Chen Z H, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren B W, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology, 29(7): 644-652. DOI:10.1038/nbt.1883
Guan W B. 2003. Studies on the structure and function of male and female reproductive system and spermatophore artifical transfer of the mud crab Scylla serrata (Forskal). Xiamen University, Xiamen. (in Chinese with English abstract)
Hansen D, Pilgrim D. 1999. Sex and the single worm:sex determination in the nematode C. elegans. Mechanisms of Development, 83(1-2): 3-15. DOI:10.1016/S0925-4773(99)00024-6
He L, Wang Q, Jin X K, Wang Y, Chen L L, Liu L H, Wang Y. 2012. Transcriptome profiling of testis during sexual maturation stages in Eriocheir sinensis using Illumina sequencing. PLoS One, 7(3): e33735. DOI:10.1371/journal.pone.0033735
Herborg L M, Rushton S P, Clare A S, Bentley M G. 2003. Spread of the Chinese mitten crab (Eriocheir sinensis H.Milne Edwards) in Continental Europe:analysis of a historical data set. Hydrobiologia, 503(1-3): 21-28.
Hopkins P M. 2009. Crustacean ecdysteroids and their receptors. In: Smagghe G ed. Ecdysone: Structures and Functions. Springer, Dordrecht, p.73-97.
Hsu S Y, Kaipia A, McGee E, Lomeli M, Hsueh A J W. 1997. Bok is a pro-apoptotic Bcl-2 protein with restricted expression in reproductive tissues and heterodimerizes with selective anti-apoptotic Bcl-2 family members. Proceedings of the National Academy of Sciences of the United States of America, 94(23): 12401-12406. DOI:10.1073/pnas.94.23.12401
Huang X S, Ye H H, Huang H Y, Yang Y N, Gong J. 2014. An insulin-like androgenic gland hormone gene in the mud crab, Scylla paramamosain, extensively expressed and involved in the processes of growth and female reproduction. General and Comparative Endocrinology, 204: 229-238. DOI:10.1016/j.ygcen.2014.06.002
Hudson C, Clements D, Friday R V, Stott D, Woodland H R. 1997. Xsox17α and -β mediate endoderm formation in Xenopus. Cell, 91(3): 397-405. DOI:10.1016/S0092-8674(00)80423-7
Hwa V, Oh Y, Rosenfeld R G. 1999. The insulin-like growth factor-binding protein (IGFBP) superfamily. Endocrine Reviews, 20(6): 761-787.
Jin S B, Fu H T, Zhou Q, Sun S M, Jiang S F, Xiong Y W, Gong Y S, Qiao H, Jiang W Y. 2013. Transcriptome analysis of androgenic gland for discovery of novel genes from the oriental river prawn, Macrobrachium nipponense, using Illumina Hiseq 2000. PLoS One, 8(10): e76840.
Kanai Y, Kanai-Azuma M, Noce T, Saido T C, Shiroishi T, Hayashi Y, Yazaki K. 1996. Identification of two Sox17 messenger RNA isoforms, with and without the high mobility group box region, and their differential expression in mouse spermatogenesis. Journal of Cell Biology, 133(3): 667-681. DOI:10.1083/jcb.133.3.667
Katakura Y. 1989. Endocrine and genetic control of sex differentiation in the malacostracan crustacea. Invertebrate Reproduction & Development, 16(1-3): 177-182.
Kelley K, Schmidt K E, Berg L, Sak K, Galima M M, Gillespie C, Balogh L, Hawayek A, Reyes J A, Jamison M. 2002. Comparative endocrinology of the insulin-like growth factor-binding protein. Journal of Endocrinology, 175: 3-18. DOI:10.1677/joe.0.1750003
Khalaila I, Manor R, Weil S, Granot Y, Keller R, Sagi A. 2002. The eyestalk- androgenic gland-testis endocrine axis in the crayfish Cherax quadricarinatus. General and Comparative Endocrinology, 127(2): 147-156. DOI:10.1016/S0016-6480(02)00031-X
Kim H W, Chang E S, Mykles D L. 2005. Three calpains and ecdysone receptor in the land crab Gecarcinus lateralis:sequences, expression and effects of elevated ecdysteroid induced by eyestalk ablation. Journal of Experimental Biology, 208(16): 3177-3197. DOI:10.1242/jeb.01754
Kopp A. 2012. Dmrt genes in the development and evolution of sexual dimorphism. Trends in Genetics, 28(4): 175-184. DOI:10.1016/j.tig.2012.02.002
Laufer H, Ahl J S B, Sagi A. 1993. The role of juvenile hormones in crustacean reproduction. American Zoologist, 33(3): 365-374. DOI:10.1093/icb/33.3.365
Li B, Dewey C N. 2011. RSEM:accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12: 323. DOI:10.1186/1471-2105-12-323
Liang S S, Zhang Z F, Yang D D, Chen Y Y, Qin Z K. 2017. Different expression of sox17 gene during gametogenesis between scallop Chlamys farreri and vertebrates. Gene Expression Patterns, 25-26: 102-108. DOI:10.1016/j.gep.2017.06.009
Liu Y, Hui M, Cui Z X, Luo D L, Song C W, Li Y D, Liu L. 2015. Comparative transcriptome analysis reveals sexbiased gene expression in juvenile Chinese mitten crab Eriocheir sinensis. PLoS One, 10(7): e0133068. DOI:10.1371/journal.pone.0133068
Livak K J, Schmittgen T D. 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods, 25(4): 402-408. DOI:10.1006/meth.2001.1262
Ma K Y, Lin J Y, Guo S Z, Chen Y, Li J L, Qiu G F. 2013. Molecular characterization and expression analysis of an insulin-like gene from the androgenic gland of the oriental river prawn, Macrobrachium nipponense. General and Comparative Endocrinology, 185: 90-96. DOI:10.1016/j.ygcen.2013.01.018
Manor R, Aflalo E D, Segall C, Weil S, Azulay D, Ventura T, Sagi A. 2004. Androgenic gland implantation promotes growth and inhibits vitellogenesis in Cherax quadricarinatus females held in individual compartments. Invertebrate Reproduction & Development, 45(2): 151-159.
Manor R, Weil S, Oren S, Glazer L, Aflalo E D, Ventura T, Chalifa-Caspi V, Lapidot M, Sagi A. 2007. Insulin and gender:an insulin-like gene expressed exclusively in the androgenic gland of the male crayfish Cherax quadricarinatus. General and Comparative Endocrinology, 150(2): 326-336. DOI:10.1016/j.ygcen.2006.09.006
Mareddy V R, Rosen O, Thaggard H B, Manor R, Kuballa A V, Aflalo E D, Sagi A, Paterson B, Elizur A. 2011. Isolation and characterization of the complete cDNA sequence encoding a putative insulin-like peptide from the androgenic gland of Penaeus monodon. Aquaculture, 318(3-4): 364-370. DOI:10.1016/j.aquaculture.2011.05.027
Nagaraju G P C. 2007. Is methyl farnesoate a crustacean hormone?. Aquaculture, 272(1-4): 39-54. DOI:10.1016/j.aquaculture.2007.05.014
Nakagawa Y, Henrich V C. 2009. Arthropod nuclear receptors and their role in molting. FEBS Journal, 276(21): 6128-6157. DOI:10.1111/j.1742-4658.2009.07347.x
Olmstead A W, Leblanc G A. 2002. Juvenoid hormone methyl farnesoate is a sex determinant in the crustacean Daphnia magna. Journal of Experimental Zoology, 293(7): 736-739. DOI:10.1002/(ISSN)1097-010X
Qiao H, Xiong Y W, Jiang S F, Fu H T, Sun S M, Jin S B, Gong Y S, Zhang W Y. 2015. Gene expression profile analysis of testis and ovary of oriental river prawn, Macrobrachium nipponense, reveals candidate reproduction-related genes. Genetics and Molecular Research, 14(1): 2041-2054. DOI:10.4238/2015.March.20.14
Rosen O, Manor R, Weil S, Gafni O, Linial A, Aflalo E D, Ventura T, Sagi A. 2010. A sexual shift induced by silencing of a single insulin-like gene in crayfish:ovarian upregulation and testicular degeneration. PLoS One, 5(12): e15281. DOI:10.1371/journal.pone.0015281
Rosen O, Weil S, Manor R, Roth Z, Khalaila I, Sagi A. 2013. A crayfish insulin-like-binding protein:another piece in the androgenic gland insulin-like hormone puzzle is revealed. Journal of Biological Chemistry, 288(31): 22289-22298. DOI:10.1074/jbc.M113.484279
Rotllant G, Nguyen T V, Sbragaglia V, Rahi L, Dudley K J, Hurwood D, Ventura T, Company J B, Chand V, Aguzzi J, Mather P B. 2017. Sex and tissue specific gene expression patterns identified following de novo transcriptomic analysis of the Norway lobster, Nephrops norvegicus. BMC Genomics, 18: 622. DOI:10.1186/s12864-017-3981-2
Senarai T, Saetan J, Tamtin M, Weerachatyanukul W, Sobhon P, Sretarugsa P. 2016. Presence of gonadotropin-releasing hormone-like peptide in the central nervous system and reproductive organs of the male blue swimming crab, Portunus pelagicus, and its effect on spermatogenesis. Cell and Tissue Research, 365(2): 265-277. DOI:10.1007/s00441-016-2375-0
Siangcham T, Tinikul Y, Poljaroen J, Sroyraya M, Changklungmoa N, Phoungpetchara I, Kankuan W, Sumpownon C, Wanichanon C, Hanna P J, Sobhon P. 2013. The effects of serotonin, dopamine, gonadotropinreleasing hormones, and corazonin, on the androgenic gland of the giant freshwater prawn, Macrobrachium rosenbergii. General and Comparative Endocrinology, 193: 10-18. DOI:10.1016/j.ygcen.2013.06.028
Smith C A, Katz M, Sinclair A H. 2003. DMRT1 is upregulated in the gonads during female-to-male sex reversal in ZW chicken embryos. Biology of Reproduction, 68(2): 560-570. DOI:10.1095/biolreprod.102.007294
Sroyraya M, Chotwiwatthanakun C, Stewart M J, Soonklang N, Kornthong N, Phoungpetchara I, Hanna P J, Sobhon P. 2010. Bilateral eyestalk ablation of the blue swimmer crab, Portunus pelagicus, produces hypertrophy of the androgenic gland and an increase of cells producing insulin-like androgenic gland hormone. Tissue and Cell, 42(5): 293-300. DOI:10.1016/j.tice.2010.07.003
Storey J D. 2003. The positive false discovery rate:a Bayesian interpretation and the q-value. Annals of Statistics, 31(6): 2013-2035. DOI:10.1214/aos/1074290335
Trapnell C, Williams B A, Pertea G, Mortazavi A, Kwan G, van Baren M J, Salzberg S L, Wold B J, Pachter L. 2010. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology, 28(5): 511-515. DOI:10.1038/nbt.1621
Tropea C, Hermida G N, López Greco L S. 2011. Effects of androgenic gland ablation on growth and reproductive parameters of Cherax quadricarinatus males(Parastacidae, Decapoda). General and Comparative Endocrinology, 174(2): 211-218. DOI:10.1016/j.ygcen.2011.08.023
Vázquez-Islas G, Guerrero-Tortolero D A, Garza-Torres R, Álvarez-Ruiz P, Mejía-Ruiz H, Campos-Ramos R. 2015. Quantitative analysis of hypertrophy and hyperactivity in the androgenic gland of eyestalk-ablated male Pacific white shrimp Litopenaeus vannamei during molt stages. Aquaculture, 439: 7-13. DOI:10.1016/j.aquaculture.2015.01.015
Ventura T, Bose U, Fitzgibbon Q P, Smith G G, Shaw P N, Cummins S F, Elizur A. 2017. CYP450s analysis across spiny lobster metamorphosis identifies a long sought missing link in crustacean development. The Journal of Steroid Biochemistry and Molecular Biology, 171: 262-269. DOI:10.1016/j.jsbmb.2017.04.007
Ventura T, Manor R, Aflalo E D, Weil S, Raviv S, Glazer L, Sagi A. 2009. Temporal silencing of an androgenic glandspecific insulin-like gene affecting phenotypical gender differences and spermatogenesis. Endocrinology, 150(3): 1278-1286. DOI:10.1210/en.2008-0906
Ventura T, Sagi A. 2012. The insulin-like androgenic gland hormone in crustaceans:from a single gene silencing to a wide array of sexual manipulation-based biotechnologies. Biotechnology Advances, 30(6): 1543-1550. DOI:10.1016/j.biotechadv.2012.04.008
Yang M L, Wu Y, Jin S, Hou J Y, Mao Y J, Liu W B, Shen Y C, Wu L F. 2015. Flower bud transcriptome analysis of Sapium sebiferum (Linn.) Roxb. and primary investigation of drought induced flowering:pathway construction and G-Quadruplex prediction based on transcriptome. PLoS One, 10(3): e0118479. DOI:10.1371/journal.pone.0118479
Young M D, Wakefield M J, Smyth G K, Oshlack A. 2010. Gene ontology analysis for RNA-seq:accounting for selection bias. Genome Biology, 11(2): R14. DOI:10.1186/gb-2010-11-2-r14
Zhang F, Lv J J, Liu P, Gao B Q, Li J, Chen P. 2015. Cloning and expression of chitinase under low salinity stress during molting in Portunus trituberculatus. Oceanologia et Limnologia Sinica, 46(4): 948-957. (in Chinese with English abstract)
Zhang F, Lv J J, Liu P, Gao B Q, Li J. 2017. Cloning and expression analysis of the cDNA of PtCht3 in Portunus trituberculatus. Progress in Fishery Sciences, 38(2): 167-176. (in Chinese with English abstract)
Zhang Y N, Xia Y H, Zhu J Y, Li S Y, Dong S L. 2014. Putative pathway of sex pheromone biosynthesis and degradation by expression patterns of genes identified from female pheromone gland and adult antenna of Sesamia inferens(Walker). Journal of Chemical Ecology, 40(5): 439-451. DOI:10.1007/s10886-014-0433-1