Adult M. japonicus were purchased from a local fish market in Okayama Prefecture, Japan. For tissue-specific gene expression analysis, the brain, eyestalk, thoracic ganglia, heart, hepatopancreas, intestine, and gonad were dissected from both three male (20.7 g average body weitht) and three female prawns (23.0 g average body weight; 1.0% average gonadosomatic index, GSI). The prawns are determined to be in intermolt (C0–C1) and early premolt (D0) stages thorugh the observation of the setal development of the pleopods using a method modified from the previous report (57). Ovarian developmental stages of the female prawns were determined as previtellogenic by histlogical analysis (20, 35). Tissues were stored in RNAlater solution (Thermo Fisher Scientific, MA, USA) at −20°C until further use. For RNA-sequencing, the ovary of two intermolt female prawns in previtellogenic stage (26.4 g body weight; 1.0% GSI) and early exogenous vitellogenic stage (56.3 g body weight; 2.5% GSI) were dissected out and stored as described above. Total RNA was isolated using the illustra RNAspin mini RNA isolation kit (GE Healthcare Bio-Sciences AB, Uppsala, Sweden).

The concentration of the two total RNA samples isolated from the ovary was measured using the Qubit RNA BR assay kit (Thermo Fisher Scientific). The cDNA libraries were constructed with 1 μg of the total RNA using the NEBNext ultra directional RNA library prep kit for Illumina, NEBNext Poly(A) mRNA Magnetic Isolation Module, and NEBNext multiplex oligos for Illumina (index primers set 1; New England BioLabs, MA, USA). All protocols were performed according to the manufacturer's instructions with minor modifications, namely that the fragmentation of RNA was performed by 94°C followed by incubation for 7.5 min. The final library fragment size was estimated to be 200–870 bp (average of 520 bp) using the Agilent high sensitivity DNA kit (Agilent Technologies, CA, USA). The library was sequenced using the MiSeq with Reagent Kit v3 (Illumina, CA, USA) in the paired-end mode with a read length of 300 bases.

Bases with a quality score (QV < 20) were trimmed from the 5′ and 3′ ends of each read, and reads containing ≥ 30% of low quality bases (QV < 14) with <25 bp were removed before assembling. The preprocessed reads were assembled using the Trinity platform (58). The resultant contigs were analyzed using the Basic Local Alignment Search Tool + (BLAST+; version 2.3.0) to perform a homology search against the National Center for Biotechnology Information non-redundant (NCBI-nr) protein database (ver. 170504) with an E-value cutoff of 1 × 10⁻⁵. The above-mentioned operations were executed using the DNA Data Bank of Japan (DDBJ) Read Annotation Pipeline and the supercomputer at the Research Organization of Information and Systems (ROIS), National Institute of Genetics (NIG), Japan. The BLAST output contained a maximum of 30 hits for each sequence which were used to assign the functional Gene Ontology (GO) terms to the protein sequences and for further GO Slim analysis on Blast2GO (59, 60). The signal peptide was predicted using SignalP 4.1 Server (61).

Poly (A)+ RNA was prepared from 50 μg of the total RNA as described above using the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England BioLabs). First-strand cDNA was synthesized from the purified Poly(A)+ RNA using the Transcriptor First Strand cDNA Synthesis kit (Roche Diagnostics, Mannheim, Germany) and an anchored-oligo(dT)₁₈ primer. This first-strand cDNA was purified with AMPure XP (Beckman Coulter, IN, USA) and tailed with poly(A) using terminal transferase (Roche Diagnostics). The following PCRs were performed to obtain the correct open reading frame (ORF) of each target gene using the synthesized cDNA.

For 5′ -RACE, PrimeSTAR HS DNA polymerase or TaKaRa LA Taq DNA polymerase (Takara Bio, Shiga, Japan) were used with the cDNA on conventional PCR programs as per a previously described method (35). Adapter1 and adapter2, as forward primers, and bursA-R01 and -R02, as reverse primers, were used for the bursicon A subunit. Adapter1, adapter2, bursA-R01, and -R02 were used for the bursicon B subunit. Adapter1, adapter2, ilp-R01, -R07, -R10, R17, and R18 were used for the insulin-like peptide (Supplementary Table 1).

All PCR products were subcloned into the pGEM-T easy vector (Promega, WI, USA) after the addition of an adenine nucleotide at the 3′ ends. All plasmids were then sequenced on the 3730xl DNA analyzer (Applied Biosystems, CA, USA).

Expression plasmids for the M. japonicus CHH-like peptide and M. japonicus neuroparsin-like peptide (Maj-pCHH-B and Maj-NPLP, respectively), were prepared as per a previously described method (62). Both Maj-pCHH-B and Maj-NPLP cDNA fragments were amplified by PCR using the pchhbexF1/R1 and nplexF/R primer pairs, respectively (Supplementary Table 1). Each PCR product was mixed with the pET-44a(+) plasmid (Novagen, WI, USA), digested using Sma I and EcoR I (New England BioLabs), purified with AMPure XP, and then ligated. The thrombin protease cleavage site of the recombinant Maj-pCHH-B (rMaj-pCHH-B) expression plasmid was modified to a tobacco etch virus (TEV) protease cleavage site using the pchhbexF2/R2 primer pair (Supplementary Figure 1). Additionally, the thrombin protease cleavage site of the recombinant Maj-NPLP (rMaj-NPLP) expression plasmid was modified to a human rhinovirus 3C (HRV 3C) protease cleavage site using the nplexF2/R2 primer pair as per a previously described method (62). These modifications accompanied the substitution of N-terminal residues from Gln to Gly in rMaj-pCHH-B and from Ala to Gly in rMaj-NPLP (Supplementary Figure 2).

The transformation and culture of E. coli strain BL21(DE3) STAR (Thermo Fisher Scientific) with expression plasmids, induction of recombinant protein overexpression, and preparation of the soluble fraction of the cell lysate were performed as per previously described methods (62, 63).

The soluble fraction of the cell lysate containing the recombinant fusion protein, His-Nus-His-tagged rMaj-pCHH-B, was purified using the Ni Sepharose 6 Fast Flow resin (GE Healthcare). While the recombinant fusion protein was captured with the resin, the affinity tags were cleaved from rMaj-pCHH-B by protease digestion in a buffer containing 25 mM Tris-HCl (pH 7.4), 0.1 mM EDTA, 0.4 M urea, and ProTEV protease (Promega) at 20°C for 24 h. The untagged rMaj-pCHH-B was washed out from the resin and further purified by reverse phase high-performance liquid chromatography (RP-HPLC) on a Capcell Pak C18 SG300 column (150 × 6 mm; Shiseido, Tokyo, Japan) using the following program: a 2-min hold at 5% acetonitrile (MeCN) in 0.05% trifluoroacetic acid (TFA), a 5-min linear gradient of 5–25% MeCN in 0.05% TFA, a 15-min gradient of 29– 37% MeCN in 0.05% TFA, a 1.2-min gradient of 37– 85% MeCN in 0.05% TFA, and a 5-min hold at 85% MeCN in 0.05% TFA at a flow rate of 0.8 mL/min.

The soluble fraction containing recombinant His-Nus-His-tagged rMaj-NPLP was incubated with the Ni Sepharose 6 resin at 4°C for 20 h. The resin was then washed with washing buffer (20 mM phosphate buffer, 0.2 M NaCl, 50 mM imidazole, pH 7.4) and equilibrated with washing buffer without imidazole. For the cleavage of the tags, HRV 3C protease (Takara Bio) was added, and the resin slurry was incubated at 4°C for 3 days. Untagged rMaj-NPLP was eluted from the resin and further purified by RP-HPLC on the aforementioned column using the following program: a 1-min hold at 5% MeCN in 0.05% TFA, a 4-min linear gradient of 5–21% MeCN in 0.05% TFA, a 15-min gradient of 21– 29% MeCN in 0.05% TFA, a 3.25-min gradient of 29–85% MeCN in 0.05% TFA, and a 5-min hold at 85% MeCN in 0.05% TFA at a flow rate of 0.8 mL/min.

The mass spectra of the purified recombinant peptides were measured on an Agilent 6,520 Accurate-Mass Quadrupole-TOF mass spectrometer with an electrospray ionization (ESI) interface (Agilent Technologies) as we have previously described (62).

The effect of one of the CHH family of peptides (Pej-SGP-I) on the mRNA expression of the putative hormones was assessed using our ex-vivo ovarian incubation system (28, 29). The same system was also used to assess the effects of rMaj-pCHH-B and rMaj-NPLP on vitellogenin (Maj-VG) expression. The recombinant Pej-SGP-I (rPej-SGP-I) was prepared as per previously described methods (62–64). Adult female prawns (22.5 g average body weight; 0.9% average gonadosomatic index, GSI) acted as donors of the ovary. As shown in Supplementary Figure 3, the abdominal part of the ovary, where left and right ovarian lobes were sticking to each other, was dissected out and divided into two lobes, and one lobe was incubated as control lobe (medium only), whereas the other received the hormone treatment (experimental). The adjacent part of ovary was kept as initial sample (without incubation) and as sample for histological analysis to determine the vitellogenic stage. Only a single sample set of the ovary (initial, control, and experimental) was prepared from one prawn and counted as n = 1; total 30 prawns were used for the experiment in Figure 5 (6 for each graph), 28 for Figure 6A, 24 for Figure 6B, and 12 for Figure 6C. All prawns used had immature ovary and were in intermolt to early premolt (C0, C1, D0, and D1) stages.

Following incubation, the ovarian tissue fragments were immersed in RNAlater solution and stored at −20°C. Total RNA extraction was then performed as described above.

Quantitative real-time reverse transcriptase PCR (qRT-PCR) was used for the quantification of the putative hormone genes, Maj-VG, and arginine kinase (Maj-AK). The sequences of primers and TaqMan probes used in this study have been listed in Supplementary Table 2. The concentrations of total RNA prepared from various tissues were quantified using the Qubit RNA BR assay kit, and for each sample 8 ng RNA was used for qRT-PCR. The qRT-PCR reactions were carried out using the iTaq universal probes one-step kit (Bio-Rad, CA, USA), and the same methods were used for the real-time monitoring of the fluorescence signal on the CFX96 real-time PCR detection system (Bio-Rad) as has been previously described (35). For the quantification of gene expression levels, relative standard curve method was used. DNA templates were amplified with respective gene-specific primer sets (Supplementary Table 2) so that they include qRT-PCR amplicon sequence of respective target genes. RNA standards were synthesized by in vitro transcription using in vitro Transcription T7 Kit (Takara Bio) with the DNA templates. The synthesized RNA standards were purified using the illustra RNAspin mini RNA isolation kit and quantified using the Qubit RNA BR assay kit as described above. The standard curves were generated using each RNA standard ranging from 40 ng to 0.4 pg prepared by serial 10-fold dilutions, and arbitrary values ranging from 40,000 to 0.4 were assigned correspondingly. Relative gene expression levels were determined based on the threshold cycles using the standard curves. Each standard had almost the same length (371–426 nt, Supplementary Table 2), and similar amplification efficiencies were achieved (95.9–99.8%) in the qRT-PCR.

For the incubated ovary samples, the relative expression of the target genes were standardized to those of Maj-AK, and expressed as a percent change from the initial (Figure 5) or 0 nM group (Figure 6) samples (35) (Supplementary Figure 3). In contrast, for tissue-specific gene expression analysis (Figure 4), the relative expression of the target genes were standardized to the total RNA input in the qRT-PCR to account for variations in the Maj-AK expression between tissues.

Gene expression levels have been represented as the mean ± standard error mean (SEM). Statistical differences in gene expression levels were analyzed using the Wilcoxon signed rank test or one-way analysis of variance (ANOVA) followed by the Dunnett's post hoc test in the GraphPad Prism version 4.0 for Windows (GraphPad Software, CA, USA).

De novo assembly produced 98,509 contigs with a total size of ~91.8 Mbp. The N50 of the transcriptome was 1,676 bp long, and the longest contig was 15,684 bp long. Redundant contigs were eliminated based on a sequence similarity. Among the remaining 47,026 contigs, 24,033 (51.1%) had a BLAST hit, and 13,839 (29.4%) were assigned at least one GO term (Supplementary File). Following bioinformatics analysis on the ovarian transcriptome, several transcripts encoding putative hormones were identified as described below.

Two contigs putatively encoding bursicon, a heterodimer composed of burs-α and burs–β, were identified in the ovarian transcriptome. Additional cDNA cloning and homology analysis revealed that these were the α and β subunits; Maj-burs-α and Maj-burs-β, respectively. The Maj-burs-α precursor contains a signal peptide (19 amino acid residues) which, once cleaved, gives a mature peptide comprising of 121 amino acid residues. Additionally, the Maj-burs-β precursor also contains a signal peptide (30 amino acid residues) and once processed results in a mature peptide composed of 115 residues (Figure 1). Mature Maj-burs-α and Maj-burs-β both showed ≥ 50% amino acid sequence identity with the known bursicon subunits, and they both contain 11 Cys residues, 9 of which are conserved in the cystine knot-like domain (smart00041 in Conserved Protein Domain Family, NCBI) of the transforming growth factor-β superfamily (Supplementary Figure 4).

The ovarian transcriptome contained one contig encoding a precursor of Maj-NPLP which contains a 24 amino acid-long signal peptide and 78 amino acid-long mature peptide (Figure 1). We identified 12 conserved Cys residues, which were a characteristic of the crustacean neuroparsin family, in the mature Maj-NPLP peptide. Maj-NPLP shares 46% amino acid sequence identity with an neuroparsin in the greasyback shrimp Metapenaeus ensis (MeNPLP) which has been hypothesized to regulate ovarian maturation (65). Additionally, Maj-NPLP shares 32% amino acid sequence identity with the ovary ecdysteroidogenic hormone (OEH) in Aedes aegypti, the yellow fever mosquito (66) (Supplementary Figure 5).

We found one contig encoding the precursor of a CHH-like peptide in the ovarian transcriptome. Additional cDNA cloning revealed that this precursor consisted of a 25 amino acid signal peptide, a 22 amino acid CHH precursor-related peptide (CPRP), an RXRR cleavage signal, and a 74 amino acid mature hormone (Figure 1). Considering the existence of CPRP and the absence of a single Gly residue 5 amino acids downstream of the first Cys residue, we considered this peptide to belong to the type-I subfamily of the CHH superfamily. On the other hand, there was no C-terminal amidation signal, a characteristic feature of type-I subfamily precursors. Additionally, there was a single-residue insertion seven amino acids downstream of the third Cys residue. Based on the primary structure nomenclature of the M. japonicus CHH-family peptides (25, 67), we categorized this peptide as a putative CHH-B (Maj-pCHH-B). Its mature peptide sequence showed 37–46% amino acid identity with the Marsupenaeus CHH-family peptides (Supplementary Figure 6). Phylogenetic analyses of the CHH superfamily characterized in the XOSG of M. japonicus so far (22, 23, 25, 26, 67) suggests that Maj-pCHH-B is diverse and most likely from the typical type-I and type-II subfamilies (Figure 2).

A putative insulin-like peptide precursor (Maj-ILP1) was identified in the ovarian transcriptome of M. japonicas. The full-length ORF of Maj-ILP1 was obtained by several rounds of 5′ -RACE. The precursor comprised a signal peptide composed of 20 residues, a B-chain consisting of 30 amino acid residues, an RXRR cleavage signal, a C-peptide composed of 131 amino acid residues, the other RXRR cleavage signal, and an A-chain composed of 24 amino acid residues (Figure 1). The primary structure of Maj-ILP1 was distinct from IAG in the same species (Maj-IAG); mature peptide (deduced B- and A- chains) of Maj-ILP1 shared only 28.6% amino acid sequence identity with that of Maj-IAG (Supplementary Figure 7). Phylogenetic analysis of known insulin/relaxin family shows that Maj-ILP1 is part of the insulin/insulin-related peptide group and not the IAG or relaxin groups (Figure 3).

Based on BLAST analysis, several contigs encoding the G-protein coupled receptor family, the receptor guanylyl cyclase family, and the insulin receptor family were found in the present ovarian transcriptome (Table 1). A full-length ORF of contig N13763 was obtained by additional cDNA cloning, and its seven-transmembrane domain showed 62% amino acid sequence identity to the neuropeptide Y (NPY) receptor in the Nevada dampwood termite Zootermopsis nevadensis. Additional cDNA cloning also revealed that the ORF of contig N28645 was composed of an N-terminal leucine-rich repeat domain and a seven-transmembrane domain that was highly similar to the transmembrane domain of the relaxin-family peptide receptors (cd15137 in Conserved Protein Domain Family). Furthermore, contigs N06137 and N35023 were the most similar to the receptor for OEH in A. aegypti (66) and to IAGR in the Chinese white shrimp Fenneropenaeus chinensis (53), respectively.

Tissue-specific expression levels of the putative hormone genes were examined by qRT-PCR analysis (Figure 4). We found that Maj-pCHH-B was mainly expressed in the nervous system and in the intestine of both male and female prawns. Expression in the ovary was low compared to the nervous system, and we did not detect any significant sexual dimorphism in the expression pattern. Additionally, Maj-NPLP was expressed primarily in the heart. Apparent expression was also observed in the thoracic ganglia and in the intestine. Maj-burs-α and Maj-burs-β are primarily expressed in the thoracic ganglia of both sexes and in the ovary. The levels of Maj-burs-α and -β expression in the ovary were significantly higher compared to those in the testis. Maj-ILP1 displayed a gender-specific expression pattern as it was expressed primarily in the ovary. Clear expression was also observed in the nervous system, while no expression was observed in the hepatopancreas of both sexes.

To investigate the potential involvement of the putative hormones in vitellogenesis, the effects of Pej-SGP-I, a VIH of M. japonicus (28), on hormone gene expression were first examined (Figure 5). In an ex-vivo ovarian incubation system, we found that Maj-pCHH-B mRNA level was significantly reduced by 50 nM rPej-SGP-I. Maj-NPLP expression was also reduced by rPej-SGP-I, but not changed (p = 0.063). Maj-burs-α, Maj-burs-β, and Maj-ILP1 mRNA levels were not affected by rPej-SGP-I.

Based on the above results, we further investigated potential involvement of Maj-NPLP and Maj-pCHH-B in vitellogenesis using recombinant peptides. rMaj-NPLP and rMaj-pCHH-B were expressed both as Nus-tagged fusion proteins; they were mostly recovered in the soluble fraction of cell lysates and successfully purified using a Ni-sepharose resin and subsequent HPLC. The deconvoluted mass spectra of untagged and HPLC-purified recombinant peptides have been shown in Supplementary Figure 8. Electrospray ionization (ESI) mass spectrum of rMaj-NPLP revealed a molecular mass of 8,357.1 which agreed with the calculated value (8,368.7) minus 12 Da, suggesting the presence of 6 disulfide bonds in the structure. Similarly, the observed molecular mass of rMaj-pCHH-B was 8,391.7 which was similar to the calculated value (8,397.6) minus 6 Da, indicating the presence of 3 disulfide bonds in the structure.

The effects of the recombinant peptides on Maj-VG expression were assessed in the ovarian incubation system. Although rMaj-NPLP did not affect Maj-VG expression at lower doses, Maj-VG expression increased at higher doses with significant changes observed in the ovary fragments treated with 20 nM rMaj-NPLP (Figure 6A). In comparison, rMaj-pCHH-B treatment did not affect the Maj-VG expression up to 200 nM (Figure 6B). When 200 nM rMaj-pCHH-B was co-incubated with 50 nM rPej-SGP-I, rMaj-pCHH-B acted neither cooperatively nor antagonistically on the vitellogenesis-inhibiting activity of rPej-SGP-I (Figure 6C).