Insights into Sexual Precocity of Female Oriental River Prawn Macrobrachium nipponense through Transcriptome Analysis

Background The oriental river prawn (Macrobrachium nipponense) is the most prevalent aquaculture species in China. The sexual precocity in this species has received considerable attention in recent years because more and more individuals matured at a small size, which devalues the commercial production. In this study, we developed deep-coverage transcriptomic sequencing data for the ovaries of sexually precocious and normal sexually mature M. nipponense using next-generation RNA sequencing technology and attempted to provide the first insight into the molecular regulatory mechanism of sexual precocity in this species. Results A total of 63,336 unigenes were produced from the ovarian cDNA libraries of sexually precocious and normal sexually mature M. nipponense using Illumina HiSeq 2500 platform. Through BLASTX searches against the NR, STRING, Pfam, Swissprot and KEGG databases, 15,134 unigenes were annotated, accounting for 23.89% of the total unigenes. 5,195 and 3,227 matched unigenes were categorized by GO and COG analysis respectively. 15,908 unigenes were consequently mapped into 332 KEGG pathways, and many reproduction-related pathways and genes were identified. Moreover, 26,008 SSRs were identified from 18,133 unigenes. 80,529 and 80,516 SNPs were yielded from ovarian libraries of sexually precocious and normal sexually mature prawn, respectively, and 29,851 potential SNPs between these two groups were also predicted. After comparing the ovarian libraries of sexually precocious and normal sexually mature prawn, 549 differentially expressed genes (DEGs) and 9 key DEGs that may be related to sexual precocity of M. nipponense were identified. 20 DEGs were selected for validation by quantitative real-time PCR (QPCR) and 19 DEGs show consistent expression between QPCR and RNAseq-based differential expression analysis datasets. Conclusion This is the first report on the large-scale RNA sequencing of ovaries of sexually precocious and normal sexually mature M. nipponense. The annotated transcriptome data will provide fundamental support for future research into the reproduction biology of M. nipponense. The large number of candidate SNPs and SSRs detected in this study could be used as genetic markers for population genetics and functional genomics in this species. More importantly, many DEGs, especially nine key DEGs between sexually precocious and normal sexually mature prawns were identified, which will dramatically improve understanding of molecular regulatory mechanism of sexual precocity of this species.


Results
A total of 63,336 unigenes were produced from the ovarian cDNA libraries of sexually precocious and normal sexually mature M. nipponense using Illumina HiSeq 2500 platform. Through BLASTX searches against the NR, STRING, Pfam, Swissprot and KEGG databases, 15,134 unigenes were annotated, accounting for 23.89% of the total unigenes. 5,195 and 3,227 matched unigenes were categorized by GO and COG analysis respectively. 15,908 unigenes were consequently mapped into 332 KEGG pathways, and many reproduction-related pathways and genes were identified. Moreover, 26,008 SSRs were identified from 18,133 unigenes. 80,529 and 80,516 SNPs were yielded from ovarian libraries of sexually precocious and normal sexually mature prawn, respectively, and 29,851 potential SNPs between these two groups were also predicted. After comparing the ovarian libraries of sexually precocious and normal sexually mature prawn, 549 differentially expressed genes (DEGs) and 9 key DEGs that may be related to sexual precocity of M. nipponense were identified. 20 DEGs were selected for validation by quantitative real-time PCR (QPCR) and 19 DEGs show consistent expression between QPCR and RNAseqbased differential expression analysis datasets.

Introduction
The oriental river prawn Macrobrachium nipponense, a member of the Palaemonidae family of decapod crustaceans, is a commercial freshwater prawn species that naturally distributed throughout most freshwaters and low-salinity estuarine regions in Japan, Korea, Vietnam, Myanmar, and China [1,2]. It is considered as an important fishery resource and widely farmed in China due to its flavor, high nutritive value and excellent adaptability, with an average annual production value of over 23 million USD [3]. However, sexual precocity has become prevalent with the expansion of production scale. Sexual precocity refers to early male and female gonad development, which has been observed in many crustaceans. Sexually precocious M. nipponense has a low value due to low growth rate, poor survival, and short life span [4][5][6], which seriously restricts the sustainable development of this species. The adverse effect of sexual precocity on female M. nipponense is particularly prominent. Sexually precocious female M. nipponenseis far less than normal sexually mature female prawn in both weight and size. Some female prawns even become sexually mature at 0.5 grams, but fecundity is only 200-500 eggs. Eggs obtained from these females often result in poor-quality offspring. Thus, understanding the mechanisms that regulate ovarian maturation in female M. nipponense and controlling sexual precocity of this prawn is crucial to improving the production of this species.
The ovary is a multifunctional organ that plays a key role in reproduction and secretion of hormones for regulation of growth and development in female prawns [7]. Ovarian maturation in prawn is a complex process controlled by several factors, such as endocrine control, nutrition and environmental factors [8][9][10][11]. However, the molecular mechanisms involved in stimulating ovarian development in prawn are still unclear. Till now, some reproduction-and ovary development-related genes have been identified from ovaries in M. nipponense, such as SUMO-conjugating enzyme (Ubc9) gene [12], heat shock protein 90 (HSP90) gene [13] and Gonadotropin-releasing hormone receptor (GnRHR) gene [14], etc. However, the genomic or transcriptomic level resources that come from M. nipponense ovary remain limited. So far, only one study has reported sequenced transcriptome from ovary of M. nipponense, which sequenced non-normalized EST library using traditional Sanger sequencing [15]. This recorded cDNA library contains 3256 ESTs that can greatly accelerate the research and discovery of new genes and molecular markers on M. nipponense ovary. However, the underlying mechanism of sexual precocity of this female prawn has not been fully revealed, especially at the molecular level, including genes and pathways. In a word, the lack of genomic and transcriptomic information of M. nipponense ovary poses an obstacle to identify genes and construct regulatory networks associated with sexual precocity of this prawn. sequencing methods such as Solexa/Illumina RNA-seq and Digital gene expression (DGE) [18], have opened a new avenue into transcriptome characterization and gene-expression profiling for various species, and rapidly dominated transcriptome studies because the higheraccuracy, higher-speed and lower-cost than the first-generation sequencing technology (Sanger sequencing). The RNA-Seq, a technique based on sequencing the poly-A RNA fraction, is a powerful tool to study complex transcriptomes because it allows for not only characterizing isoforms from known genes but also discovering novel or predicted coding genes [19]. It gives a general view of gene expression especially in these species lack of a fully sequenced and assembled genome such as M. nipponense. Some cDNA libraries and transcriptome-level datasets have been obtained from M. nipponense by RNA-Seq to lay a foundation for functional genomics approaches used for improving the aquaculture performance of this species [2,20,21]. Based on these transcriptome studies, there are about 81,411 expressed sequence tags (ESTs) from M. nipponense in the public databases up to date. However, there have been no transcriptome studies regarding the ovary of sexually precocious M. nipponense was reported until now.
In the present study, we performed high-throughput sequencing of the ovaries of sexually precocious and normal sexually mature M. nipponense using Illumina RNA-Seq to generate a transcriptome database that will enlarge the public EST database for this species and help support future studies. The identification of differentially expressed genes and pathways in the ovary of these two types of prawn will help build a more complete understanding of the regulatory mechanisms associated with sexual precocity. In addition, the simple sequence repeats (SSRs) and single nucleotide polymorphisms (SNPs) reported in this transcriptome study are also potentially useful for population genetics and functional genomics studies in this species.

Sample preparation and RNA extraction
There were two groups of female experimental prawn, one group was sexually precocious M. nipponense (MNOP) (2.5-3.5 cm, 0.5-1.5 g) which has grown about 90 days from hatching to sexual maturity, another group was normal sexually mature M. nipponense (MNON) (4.5-5.5 cm, 2.5-4.5 g) which took about one year to reach sexual maturity after hatching. Each group had ten samples which were all obtained from aquaculture base in Zhejiang institute of freshwater fisheries, Huzhou, China, and their ovaries were all at the stage V (maturation stage) based on the external morphology, color, gonado-somatic index (GSI), and histological features [13,15]. After prawns were anesthetized and finally sacrificed by placing them in an ice bath, ovaries were dissected and snap frozen in liquid nitrogen and stored at-80°C until use. The study was approved by the Institutional Animal Care and Use Ethics Committee of Agriculture Ministry Key Laboratory of Healthy Freshwater Aquaculture, Zhejiang Institute of Freshwater Fisheries (Huzhou, China). Total RNA was extracted using Trizol Reagent (Invitrogen, Carlsbard, CA, USA). The integrity and purity of the total RNA were detected by an Agilent-2100 bioanalyzer (Agilent Technologies, Inc., Santa Clara CA, USA) and its concentration was determined using the NanoDrop 2000 spectrophotometer (Nano-Drop, Thermo Scientific, Wilmington, DE, USA). Ten RNA samples from each group were pooled in equal amounts, generating two mixed RNA samples from MNOP and MNON respectively, which were then used to prepare two separate RNA-seq transcriptome libraries.
Transcriptome library preparation and Illumina sequencing RNA purification, reverse transcription, library construction and sequencing were performed at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China) according to the manufacturer's instructions (Illumina, San Diego, CA). The two RNA-seq transcriptome libraries were prepared using Illumina TruSeqTM RNA sample preparation Kit (San Diego, CA). Poly(A) mRNA was purified from total RNA using oligo-dT-attached magnetic beads and then fragmented by fragmentation buffer. Taking these short fragments as templates, double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). Then the synthesized cDNA was subjected to end-repair, phosphorylation and 'A' base addition according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 200-300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (New England Biolabs, Boston, MA) for 15 PCR cycles. After quantified by TBS380, two RNAseq libraries were sequenced in single lane on an Illumina Hiseq 2500 sequencer (Illumina, San Diego, CA) for 2×100bp paired-end reads [22,23].

Sequencing data assembly and functional annotation
Image data output from the sequencer was transformed into raw reads by base calling, and stored in FASTQ format. The raw reads were filtered into the high quality clean reads after removing adaptor sequences, empty reads, and low quality sequences (reads with unknown sequences 'N' or less than 20bp) using SeqPrep (https://github.com/jstjohn/SeqPrep)and Sickle (https://github.com/najoshi/sickle). The clean reads were then assembled de novo into transcripts using the Trinity program (http://trinityrnaseq.sourceforge.net/, version number: trinityrnaseq-r2013-02-25) [24]. To eliminate redundant sequences, transcripts were clustered based on sequence similarities and the longest transcript in each cluster represented the final unique sequence (unigene). All the unigenes were searched against the protein databases such as non-redundant (NR) protein sequence database, Kyoto Encyclopedia of Genes and Genomes (KEGG) database, Search Tool for the Retrieval of Interacting Genes (String) database, Swissprot and Pfam databases to identify the proteins with the highly sequence similarity to the given unigenes and retrieve their function annotations, using BLASTX with a typical cut-off E-value <1e -5 . Biological processes, molecular functions and cellular components were described by GO annotations of unigenes using BLAST2GO (http://www.blast2go.com/ b2ghome) [25]. Clusters of Orthologous Groups (COG) (http://www.ncbi.nlm.nih.gov/COG/) and KEGG (http://www.genome.jp/kegg/) were used to predict possible functional classifications and metabolic pathways, respectively [26,27].

Molecular Marker Detection
MSATCOMMANDER V. 0.8.2 (http://code.google.com/p/msatcommander/) [28] was used to search for microsatellites (simple sequence repeats, SSRs). The minimum repeat lengths considered were as follows: ten repeats for mono-SSRs; six repeats for di-SSRs; and four repeats for tri-, tetra-, penta-, and hexa-SSRs. SSR-containing unigenes were annotated based on BLAST similarity searches. The location of SSRs was estimated based on open reading frames (ORFs) and untranslated regions (UTRs) within the unigenes that were identified using Trinity [24]. All the consensus assembly sequences generated from the two transcriptome libraries were employed as reference sequences to detect potential single nucleotide polymorphisms (SNPs) in MNON and MNOP using Samtools V. 0.1.18(http://samtools.sourceforge.net) [29] and VarScan v.2.2.7(http://varscan.sourceforge.net/) [30]. In addition, to identify SNPs between MNON and MNOP groups, all sequences generated from the MNON transcriptome library were employed as reference sequences and SNPs in MNOP were detected using Samtools V. 0.1.18 and VarScan v.2.2.7.

Identification and annotation of differentially expressed genes
We used the RNASeq by Expectation Maximization (RSEM) (http://deweylab.biostat.wisc.edu/ rsem/) [31] to calculate unigene abundances from RNA-Seq data. FPKM (fragments per kilobase of exon model per million mapped reads) was used to quantify the gene expression levels. FPKM value can be directly used for the comparison of gene expression level between samples. FPKM is calculated as follows: FPKM ¼ total exon Fragments mapped reads ðMillionsÞ Â exon length ðkbÞ In this formula, total exon Fragments is the number of reads that aligned to a specific unigene, mapped reads (Millions) is the total number of reads that aligned to all unigenes, exon Length (kb) is the length of the unigene.
The detection of differentially expressed genes (DEGs) between ovaries of MNON and MNOP was performed using edgeR package (Empirical analysis of Digital Gene Expression in R) (http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) [32] with the rigorous algorithm method. Following the procedures described in the edgeR documentation, read count tables were loaded into R, normalized using the default method for edgeR (trimmed mean of M values, or TMM), and then tested for differential expression using the Fisher exact Test method. 0.1, a typical value for the common biological coefficient of variation (square root-dispersion) for datasets arising from well-controlled experiments, was chosen as a nominal dispersion value for the data without biological replicates. The P-values reported by edgeR were adjusted for multiple hypothesis testing using the False Discovery Rate (FDR) control method [33], and the ratio of FPKMs was converted to the fold-change in the expression of each gene in two samples simultaneously. FDR < 0.05 and the absolute value of logFC (Fold change) ! 1 were set as the threshold for significantly differential expression. The function of DEGs was annotated by the GO and KEGG pathway enrichment analyses at a P 0.05.

Quantitative real-time PCR (QPCR) validation
Twenty candidate DEGs (Thirteen annotated DEGs and seven unannotated DEGs) were randomly selected for validation by QPCR. Total RNA was extracted from ovaries of MNON and MNOP using Trizol Reagent (Invitrogen, Carlsbard, CA, USA), and the RNA samples were reverse-transcribed using PrimeScript 1 RT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Dalian, China). Oligonucleotide primers (S1 Table) were designed with the Primer Premier 5.0 software and synthesized by Invitrogen (Shanghai, China). QPCR was performed using SYBR 1 Premix Ex Taq™ II (Tli RNaseH Plus) (TaKaRa, Dalian, China) on a CFX96 TM Real-Time PCR Detection System (Bio-Rad, USA). Each PCR reaction (25 μL) contained 12.5 μL SYBR 1 Premix Ex Taq II, 0.4 μM of each primer and appropriately diluted cDNA. The thermal cycler program was 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Six biological and three technical replicates were performed respectively. Beta-actin (β-actin) was amplified in parallel as an internal control for normalization of QPCR data due to its stable expression in M. nipponense after the expressions of three reference genes including β-actin, elongation factor 1-alpha (Ef1α) and ribosomal protein L8 (RpL8) were compared in different tissues of MNOP and MNON by the method of Qian et al (2014) [34]. The relative changes in gene expression levels were calculated using the 2 −ΔΔCT method [35,36]. Statistical analyses were performed with SPSS 20.0 software (SPSS Inc., Chicago, IL, USA). Results were presented as the means±S.D. (n = 6) for each group. One-way analysis of variance (ANOVA) followed by the post hoc test were carried out to determine whether the differences between groups were significant (P < 0.05).

Transcriptome sequencing and assembly
Two cDNA libraries were constructed from ovaries of MNOP and MNON. The Illumina Hiseq 2500 sequencing generated 47,254,012 raw reads from the MNOP and 56,793,358 raw reads from the MNON. The raw sequencing data were deposited in the Short Read Archive (SRA) of the NCBI with accession numbers of SRP063589. After filtering out the adaptor primers and empty and low-quality reads, 42,115,012 and 50,789,578 clean reads were generated from MNOP and MNON libraries, respectively (Table 1). After de novo assembly of clean reads using Trinity, a total 63,336 unigenes were identified. The lengths of all unigenes were distributed as 893.88 bp, 28,661 bp, 201 bp and 1846 bp from average, longest, shortest and N50 length, respectively ( Fig 1A). Of these, 31,441 (49.64%) were 1-400 bp; 9,273 (14.64%) were 401-600 bp; 4,634 (7.32%) were 601-800 bp; 3,021 (4.77%) were 801-1000 bp and 2,187(3.45%) were 1001-1200 bp ( Fig 1B).

Discussion
Little is known about the genes and biological mechanisms controlling sexual precocity in crustaceans. In this study, by profiling the ovary transcriptomes of sexually precocious and normal sexually mature M. nipponense, we expand the limited amount of sequence data that are available for M. nipponense in the public EST database and identified differentially expressed genes potentially influencing the sexual precocity of this prawn.
In this study, the RNA sequencing produced total 63,336 unigenes with average length was 893.88 bp from the ovaries of MNOP and MNON. This is considerably higher than the sum of unigenes (1,514) and average length of unigenes (734) obtained from ovary of M. nipponense using traditional method of expressed sequence tag (EST) reported by Wu et al. (2009). For one thing, the results show that high-throughput RNA-Seq technology has great superiority in   Wu et al. (2009). To the best of our knowledge, this is the first comprehensive ovary transcriptome study of not only normal sexually mature but also sexually precocious M. nipponense. Therefore, this transcriptome dataset provides a useful resource for future analyses of genes related to reproduction, gonadal development and sexual precocity in M. nipponense. A total of 15,134 (23.89%) unigenes were annotated successfully in the present study, whereas the rest 48,202 (76.11%) unigenes had no annotations which will make a meaningful contribution to the knowledge of M. nipponense by further characterizing their functions.  BLAST results against the Nr database showed that the M. nipponense sequences had the highest number of matches with the gene sequences from Z. nevadensis, followed by D. pulex. Z. nevadensis was the latest arthropod of the class insecta and D. pulex was the first arthropod of the class crustacea to have their whole genome sequenced [43,44]. However, only about 11.92%    transcriptome research in L. vannamei was the most extensive of all decapod crustaceans of the suborder Natantia in recent years [45][46][47][48][49][50], therefore, it is not surprising that L. vannamei was the only decapod crustaceans of the suborder Natantia in the top twenty species that match to the sequences of M. nipponense. The results of GO and KEGG analysis extend our knowledge on the gene function classification and potential involvement of genes in cellular metabolic pathways [51,52]. In this study, the M. nipponense ovarian sequences fall into GO categories with a roughly similar distribution to that of some other arthropods, such as Drosophila melanogaster, Parhyale hawaiensis, D. pulex [53] and Oncopeltus fasciatus [54], in which most ovarian sequences were enriched in functional category of "catalytic activity", suggesting that arthropod ovaries contain a large diversity of enzyme genes involved in various molecular function. As might be expected, many reproduction-related pathways and genes were found in the KEGG analysis, which maybe play Fig 11. QPCR analysis of twenty selected DEGs. X-axis: The 20 differentially-expressed unigenes or genes, Y-axis: Relative expression level of each unigenes or genes. Beta-actin is as internal control gene. * indicates significantly differential expression (P < 0.05), ** indicates most significantly differential expression (P < 0.01).
doi:10.1371/journal.pone.0157173.g011 an important role in development of oocytes and ovaries in crustaceans. However, further studies are also needed to understand the molecular functions of some putative reproductionrelated genes in these pathways in M. nipponense.
The focus of our project was to mine the differentially expressed genes (DEGs) between MNOP and MNON. In the present study, we identified 549 DEGs which would give clues on the molecular mechanism of sexual precocity determination. In these DEGs, the number of down-regulated unigenes (337) was more than the number of up-regulated unigenes (212), as well as the number of MNON-specifically expressed unigenes (187) was more than the number of MNOP-specifically expressed unigenes (90), indicating that the lack or low expression of quite a few genes in the ovarian development of M. nipponense is one of the important causes of sexual precocity of this species. Besides, in these DEGs, the number of unannotated sequences (363) was much larger than the annotated ones (186), which made it impossible to exploit all DEGs, some potentially useful genetic information relate to sexual precocity of M. nipponense may be missed. From these 186 annotated DEGs, we identified nine genes of interest that may be play an important role in sexual precocity of M. nipponense. They were VTG, CTSL, CST6, INSR, IGF1R, COX2, GPX, SOD1 and FASN genes.
VTG, CTSL, CST6 are 'ovary development'-related genes. The expressions of these three genes were significantly up-regulated in MNOP compared with MNON in the present study. VTG is the precursor of vitellin (VN), the major yolk protein that provides nutrition during embryonic development. VTG is synthesized extraovary or intraovary, transported by the hemolymph and finally sequestered into the growing oocytes through vitellogenin receptor (VTGR) mediated endocytosis [55,56]. In oviparous animals, oocyte development and ovarian maturation depend on massive production and accumulation of VTG. In vertebrates, estrogen hormones play a key role in regulating the synthesis of VTG through the enhancement of the transcription of the VTG gene [57,58]. This process is mediated by an E-ER-HSP90-VTG regulation pathway which widely existed in vertebrates [56][57][58][59][60][61]. During vitellogenesis, estrogen (E) binds to its receptor (ER) form a hormone-receptor complex and then acts on estrogen responsive elements (ERE) located at the upstream of the VTG DNA, which leads to the activation or enhancement of the VTG gene transcription [61]. However, HSP90 binds to estrogen receptor (ER) in the absence of estrogen (E) to increase the activity of estrogen hormone receptor complex to transcribe target genes [59,60,62,63]. In our transcriptomic data, HSP90 genes were identified. Although an ER gene has not been found yet, we identified estrogen related receptor (ERR) gene in our transcriptome database (S5 Table). However, HSP90 and ERR genes had no significantly differences between MNOP and MNON. So, whether estrogen (E) caused the significant up-regulation of VTG gene expression in MNOP pass through the E-ER-HSP90-VTG regulation pathway need to be confirmed by more studies. After all, the regulatory mechanism of E-ER-HSP90-VTG pathway for VTG synthesis in crustaceans is still controversial [13,37]. Cathepsins are lysosomal cysteine proteases, but cystatins are natural tight-binding reversible inhibitors of cysteine proteases. Based on sequence homology and conserved amino acid motifs, the group of cathepsins can be divided into three subgroups: cathepsin L (CTSL), cathepsin B (CTSB) and cathepsin F (CTSF) [64]. CTSL plays an important role during yolk formation and processing in vertebrates, for example, CTSL can degrade VTG into smaller yolk proteins which are stored in oocyte and will be the primary source of nutrients and energy for the developing embryo [65]. Current findings in Metapenaeus ensis [66] and M. nipponense [67] suggest that the ovary is the unique non-digestive organ showing CTSL expression and CTSL in the ovary in crustacean is also functionally important in reproductive physiology. On the contrary, as a cysteine protease inhibitor, CST6 plays an important role in regulating and protecting against uncontrolled proteolysis by cathepsins L (CTSL) and V (CTSV) and the asparaginyl endopeptidase legumain (LGMN) [68]. In the present study, the significant up-regulation of CST6 gene expression in MNOP ovary means that the large amounts of CST6 produced and inhibited the hydrolysis of VTG by CTSL. Although CTSL also produced abundantly because of the significant up-regulation of CTSL gene expression, the inhibition of CST6 probably overwhelmed the hydrolysis of CTSL. In sum, the magnificent VTG production due to the significant up-regulation of VTG gene expression, together with the inhibition of the hydrolysis of VTG by CST6, led to the accumulation of large amounts of VTG in the ovary, and finally resulted in the sexual precocity of M. nipponense.
INSR, IGF1R, COX2 genes were three DEGs in 'ovarian steroidogenesis' pathways (Table 3). Compared with MNON, INSR and IGF1R gene expressions were significantly up-regulated, while COX2 gene expression was significantly down-regulated in MNOP. INSR and IGF1R are transmembrane receptors that are activated by hormones called insulin and insulin-like growth factor 1 (IGF-1) respectively. Insulin is a key peptide hormone that regulates growth, metabolism and reproduction in vertebrates [69], However, insulin-like peptides (ILPs), the counterpart of vertebrate insulin, belong to the insulin-like superfamily and were widely found in invertebrate species [70][71][72]. The insulin-like superfamily includes two major subgroups, ILPs and the insulin-like growth factors (IGFs) [72]. IGF-I is a member of this superfamily, its function overlaps with that of insulin/ILPs [73]. In vertebrates, insulin and IGF-1 banded to the INSR and IGF1R respectively, and then they can stimulate DNA synthesis and promote the proliferation of ovarian theca-interstitial cell through insulin signaling pathway [74]. Insulin signaling pathway is also present in invertebrates, and is closely related to their reproduction. For instance, insulin-like growth factor signaling pathway in Drosophila was considered important for the regulation of vitellogenesis, oogenesis, and reproductive diapause [75]. The interaction and crosstalk between ILP signaling networks in Aedes aegypti was thought to be relevant to the regulation of yolk protein precursor, juvenile hormone, ecdysteroids, and nutrients [76]. In the present study, the significantly up-regulation of INSR and IGF1R gene expressions in MNOP means that a great deal of ILPs and IGF-1 generated in the sexually precocious prawn, which promoted vitellogenesis, oogenesis and proliferation of ovarian theca-interstitial cell, leading to the rapid ovarian maturation of M. nipponense. Cyclooxygenase-2 (COX2) is an inducible enzyme that catalyzes the formation of prostaglandins (PGs) from arachidonic acid (AA). In crustaceans, PGs have important significance in the regulation of female reproductive maturation. In Procambarus paeninsulanus [77] and Oziothelphusa senex senex [78], there has a positive correlation between the production of PGs (PGE 2 and PGF 2α ) and the ovarian maturation. Furthermore, the amount of PGs in ovaries of crustaceans varies along with developing ovary stages. For instance, the amounts of PGE 2 and PGF 2α were highest in ovaries at stage I and continued to decrease until stage IV in Metanephrops japonicus [79]. In this study, COX2 gene expression was significantly down-regulated in ovary at stage V of MNOP compared with MNON, suggesting the decrease of PGs in ovary at stage V of sexually precocious M. nipponense. We believe that the expression of COX2 gene in immature ovaries (stageI-VI ovaries) of sexually precocious M. nipponense was higher than that in mature ovary (stage V ovary) as M. japonicus, which promoted the synthesis of the PGs in large quantities in immature ovaries, and finally resulted in the rapid ovarian maturation. However, due to the reduced need for PGs in mature ovary, the expression of COX2 gene decreased sharply in the stage V ovary of sexually precocious M. nipponense, even significantly lower than that in the stage V ovary of normal sexually mature prawn. But the real expression levels of COX2 gene in ovaries at different stages of sexually precocious and normal sexually mature M. nipponense require further investigation.
Aerobic metabolism in animals is associated with the generation of free radicals or reactive oxygen species (ROS) that include the hydroxyl radicals, superoxide anion, hydrogen peroxide, and nitricoxide. ROS have been recently recognized as key signaling molecules involved in reproductive biology including oocyte maturation, fertilization, and embryo development [80][81][82]. However, because of the double-edged sword property of ROS, excessive ROS can also lead to serious damage to major cellular components including protein, lipid and DNA [83]. To minimize the oxidative damage caused by ROS, cells possess a range of enzymatic scavenger systems including superoxide dismutase (SOD), catalase (CAT) and glutathione peroxidase (GPX). In this study, SOD and GPX belong to 'peroxisome' and 'glutathione metabolism' pathways, respectively, while they are important antioxidant enzymes in antioxidant system. SOD first catalyzes the dismutation of superoxide radicals to H 2 O 2 , which is further metabolized to H 2 O and O 2 by CAT and GPX [84]. GPX, a selenoprotein, plays an important role in detoxifying lipids and hydrogen peroxide, with the concomitant oxidation of glutathione [85]. In this study, compared with MNON, SOD and GPX gene expressions in MNOP were significantly down-regulated and up-regulated respectively, while CAT gene expression had no significant difference between these two groups. There must be a precise balance in ROS production and recycling under the regulation of these antioxidant enzymes for avoiding the oxidative damage to cells in ovary of M. nipponense. The Low expression of SOD gene and the over expression of GPX gene disrupted this balance and caused overproduction of ROS. The excess ROS promoted the mature of the ovary as signaling molecules, and finally leading to the sexual precocity of M. nipponense.
It was reported that lipids (including all kinds of fatty acids) are the essential nutrients which stored in the oocytes of female broodstocks to nourish subsequent embryonic development during ovarian maturation of prawn [86]. Moreover, fatty acids also have effects on growth and reproductive success in female prawns through stimulating ovarian maturation and oocyte differentiation [87,88]. For example, C-20 polyunsaturated fatty acids (PUFAs), include arachidonic acid (AA), eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA), serve as prostanoid precursors and convert into a variety of prostanoids via enzymatic catalysis. One of the more prominent roles of prostanoids in crustaceans is the regulation of female reproductive maturation [89]. Fatty acid synthase (FASN), a key enzyme in the process of the de novo fatty acid synthesis in 'fatty acid biosynthesis' pathway, catalyzes all of the reaction steps for synthesizing saturated fatty acids (SFA) from acetyl-CoA and malonyl-CoA in an NADPH-dependent manner [90]. Then, various unsaturated fatty acids (USFA) are synthesized from SFA under catalysis of the fatty acid desaturases. In this study, various fatty acids and prostanoids in ovary of MNOP should be richer than those in ovary of MNON because of the up-regulation of FASN gene expression in ovary of MNOP compared with MNON, which may be one of the key factors of sexual precocity of M. nipponense.
Conventional Quantitative real-time PCR was frequently used to confirm the data obtained from high-throughput sequencing. Nineteen out of twenty selected DEGs from RNA-Seq were confirmed by QPCR, which indicates the high reliability of RNA-Seq for detecting DEGs. The differential expressions of above nine key genes that may be related to sexual precocity of M. nipponense were also verified with QPCR (Fig 11). Although these nine genes are highly related to reproduction of crustaceans, their involvement in sexual precocity of crustaceans has not yet been reported. Our findings in this study firstly revealed the potential relationship between these nine genes and sexual precocity of M. nipponense. Nevertheless, conclusion on this point cannot be drawn completely until the regulation pathways of these nine genes in sexual precocity of M. nipponense are fully understood.
SSRs and SNPs are useful molecular markers that have been widely applied into genetic mapping, population genetic analysis and breeding [35]. In the present study, we identified large volumes of SSRs and SNPs which were a lot higher than those in previous reports on M. nipponense [2,20,21]. This is a powerful supplement to SSR and SNP databases available for M. nipponense. Moreover, further in-depth research of the potential SNPs between MNON and MNOP groups predicted in this study will also provide certain help for the understanding of the sexually precocious mechanism of this prawn.

Conclusion
This is the first comprehensive transcriptome dataset for the ovaries of sexually precocious and normal sexually mature M. nipponense to be reported. A total of 63,336 unigenes were obtained and a number of reproduction-related unigenes and metabolic pathways were identified. Many potential SSRs and SNPs were predicted and can be used for subsequent marker development, genetic linkage and gene localization analysis. Moreover, 549 DEGs and nine key DEGs that are potentially involved in sexual precocity of M. nipponense were identified for the first time and are worthy of further investigation. In sum, these findings will serve as the important foundation for insight into the regulatory mechanism of sexual precocity of M. nipponense, and may also provide some references for future molecular research in freshwater prawn.
Supporting Information S1