- Research article
- Open Access
Transcriptome characterization of BPG axis and expression profiles of ovarian steroidogenesis-related genes in the Japanese sardine
BMC Genomics volume 21, Article number: 668 (2020)
The clupeoid fishes are ecologically and commercially important fish species worldwide that exhibit a high level of population fluctuation, accompanied by alteration of reproductive traits. However, knowledge about their reproductive physiology in order to understand mechanisms underlying such population dynamics is limited. The endocrine system along with the brain–pituitary–gonadal (BPG) axis is critical for regulating reproduction. The aims of this study were to provide transcript data and genes related to the BPG axis, and to characterize the expression profiles of ovarian steroidogenesis-related genes in the Japanese sardine (Sardinops melanostictus, Clupeidae).
RNA sequencing was performed using the sardine brain, pituitary, and gonad in both sexes. A total of 290,119 contigs were obtained and 115,173 non-redundant ORFs were annotated. The genes differentially expressed between ovary and testis were strongly associated with GO terms related to gamete production. The tissue-specific profile of the abundance of transcripts was characterized for the major regulators in the BPG axis, such as gonadotropin-releasing hormone, gonadotropin, and steroidogenic enzyme. By comparing between ovary and testis, out of eight different 17β-hydroxysteroid dehydrogenase (Hsd17b) genes identified, higher hsd17b7 expression was found in testis, whereas higher expression of hsd17b8, hsd17b10, hsd17b12a, and hsd17b12b was found in ovary. The cDNAs encoding key endocrine factors in the ovarian steroidogenic pathway were cloned, sequenced, and quantitatively assayed. In the pituitary, follicle-stimulating hormone beta peaked during vitellogenesis, while luteinizing hormone beta peaked at the completion of vitellogenesis. In the ovary, follicle-stimulating hormone receptor and luteinizing hormone receptor were upregulated from mid- to late phase of vitellogenesis. Furthermore, three steroidogenic enzyme genes (cyp11a1, cyp17a1, and cyp19a1a) gradually increased their expression during ovarian development, accompanying a rise in serum estradiol-17β, while 3β-hydroxysteroid dehydrogenase and steroidogenic acute regulatory protein did not change significantly.
This is the first report of deep RNA sequencing analysis of Japanese sardine, in which many key genes involved in the BPG axis were identified. Expression profiles of ovarian steroidogenesis-related genes provide a molecular basis of the physiological processes underlying ovarian development in the sardine. Our study will be a valuable resource for clarifying the molecular biology of clupeoid fishes.
The application of next-generation sequencing (NGS) technology has recently increased in the field of aquaculture and fisheries. For aquaculture fish species, the administration of exogenous hormones and environmental manipulation are conducted to control gonadal development in captivity. Since knowledge of reproductive physiology is helpful to develop effective treatments for such control, RNA sequencing (RNA-seq) has been applied for a variety of aquaculture species to reveal the physiological mechanisms regulating reproduction. As recent examples, RNA-seq analyses were conducted on carp, eel, perch, and salmonid species [1,2,3,4,5]. For important target species of fisheries, the application of large-scale genomic approaches has mostly been conducted to investigate population structure and genetic differentiation (e.g., for cod  and tuna species ). Accordingly, the current understanding of large-scale transcriptome data for fish reproduction is considered to be biased towards aquaculture species. Meanwhile, recent research on swordfish (Xiphias gladius) fisheries showed that RNA-seq analysis of ovarian development can be a useful tool in establishing a biological basis of reproduction for successful stock management .
Clupeoid fishes include ecologically and commercially important fish species, such as anchovies, herrings, and sardines. They exhibit dramatic and cyclic population dynamics in response to climate change at the multi-decadal scale [9, 10]. Since such population fluctuation is accompanied by alteration of reproductive traits , knowledge about the reproductive physiology is key to understanding the mechanisms of such population dynamics. The first draft genome of Atlantic herring (Clupea harengus) was published in 2016 by Martinez Barrio et al.  and the genome assembly has recently been improved by Kongsstovu et al. ; however, gene profiling related to reproduction has not been characterized.
As in other vertebrates, fish reproduction is regulated by the activation of endocrine factors in the brain–pituitary–gonadal (BPG) axis. The brain gonadotropin-releasing hormone (Gnrh) stimulates the secretion of pituitary gonadotropins (Gths), follicle-stimulating hormone (Fsh) and luteinizing hormone (Lh), which in turn act on the gonads to stimulate the production of sex steroid hormones, such as estrogens and androgens [14,15,16]. The Gnrhs are decapeptides, of which teleost fish commonly have two or three forms within the brain of individual species, that is, species-specific-type Gnrh (Gnrh1), chicken-II-type Gnrh (Gnrh2), and salmon-type Gnrh (Gnrh3) . The Gths are glycoprotein hormones composed of a common glycoprotein α-subunit (Cga) and hormone-specific β-subunit (Fshb or Lhb) .
During oogenesis, the substantial growth of the oocyte is marked by vitellogenesis (yolk protein accumulation), which occurs after previtellogenic growth with the appearance and accumulation of cortical alveoli . The progress of vitellogenesis is generally accompanied by a rapid elevation of circulating estradiol-17β (E2), which induces hepatic vitellogenin production . E2 is produced in oocyte-surrounding follicle cells which are composed of an inner layer of granulosa cells, a basement membrane, and an outer layer of thecal cells . In the follicle cells, the physiological effects of Gths are mediated by Gth receptors (Gthrs; Fshr and Lhr) . In the ovarian steroidogenic pathway, a series of enzymes, such as cytochrome P450 side-chain cleavage (cyp11a1), 17α-hydroxylase/17,20-lyase (cyp17a1), ovarian aromatase (cyp19a1a), and 3β-hydroxysteroid dehydrogenase (hsd3b), as well as steroidogenic acute regulatory protein (star), catalyze the conversion of cholesterol to E2 . Information on the molecular regulation of steroidogenic proteins (steroidogenic enzymes and Star) in the ovarian steroidogenic pathway is a prerequisite for understanding physiological mechanisms governing ovarian development. In clupeoid fishes, there is no information available on the expression of genes encoding ovarian steroidogenic proteins.
The objective of the present study was to provide transcript data and genes related to the BPG axis with the goal of establishing a molecular basis for the endocrine regulation of ovarian development using the Japanese sardine (Sardinops melanostictus, Clupeidae). The Japanese sardine has been a key species for studies of population dynamics because of its high level of population fluctuation . Although some evidence of the mechanisms of population dynamics has been presented from an ecological perspective [21,22,23,24,25], the available information on aspects of the endocrine regulation of ovarian development is still limited. RNA-seq was performed using sardine brain, pituitary, and gonad from both sexes to generate a high-quality transcriptome assembly. Furthermore, complementary DNAs (cDNAs) encoding key endocrine factors in the ovarian steroidogenic pathway were cloned and sequenced. We then analyzed their expression profiles during ovarian development by quantitative real-time PCR.
Assembly and functional annotation
Japanese sardines at age 3 (three fish for each sex) sampled in February were used for RNA-seq. The body size was 186–203 mm in body length (BL) and 91–112 g in body weight (BW). Females had vitellogenic ovaries, while males were at the spermiation stage.
Sequencing of the cDNA libraries using NextSeq500 yielded a total of 5–7 million pared-end reads (150 bp) per organ, which were deposited in the DNA Data Bank of Japan (DDBJ) Sequence Read Archive (DRA) under accession numbers DRA010129. The adapter and low-quality sequences were trimmed, and then the remaining reads from all libraries were assembled by Trinity into 290,119 contigs with an average length of 696 bp (Table 1). Open reading frames (ORFs) of more than 300 bp were extracted from the contigs and then the redundant protein sequences with up to 95% similarity were clustered into 115,173 distinct ORFs, which consisted of complete (47,886), internal (24,524), 5′ partial (30,603), and 3′ partial (12,160) ORFs (Table 1). Those translated ORFs were then used for sequence homology searches. Among a total of 115,173 protein sequences, 47,201 (41.0%), 46,335 (40.2%), 45,307 (39.3%), and 54,667 (47.5%) protein sequences showed significant similarities (e-value ≤1e−05) to protein sequences of Atlantic herring, zebrafish (Danio rerio), medaka (Oryzias latipes), and the NCBI non-redundant protein (nr) database, respectively (Fig. 1a). Of these, 24,037, 24,388, 22,118, and 38,651 protein sequences are allocated unique accession numbers for Atlantic herring, zebrafish, medaka, and the NCBI nr database, respectively (Fig. 1a). Overall, 888, 104, 80, and 7194 protein sequences had unique similarity to Atlantic herring, zebrafish, medaka, and the NCBI nr database, respectively, and 42,525 protein sequences shared similarity among the three species and the NCBI nr database (Fig. 1a). The cumulative relative frequency distribution for amino acid identities between sardine and the three fish species showed that, at 80% amino acid identity, the frequencies were 47.2% (Atlantic herring), 75.5% (zebrafish), and 80.9% (medaka) (Fig. 1b).
Additionally, a total of 9053 Gene Ontology (GO) numbers were assigned to 19,722 protein sequences (17.1%), and 828 enzyme commission (EC) numbers were assigned to 4588 protein sequences (4.0%). All of the contigs were registered with DDBJ as Transcriptome Shotgun Assembly (TSA) sequences under accession numbers ICPT01000001–ICPT01290119. The results of the Protein Basic Local Alignment Search Tool (BLASTP) are shown in Additional file 1. The statistical distribution of GO-slim categories is shown in Additional file 2: Figure S1.
Differentially expressed genes
The numbers of transcripts commonly and differentially expressed in the three tissues (brain, pituitary, and gonad) for female and male were analyzed (Fig. 1c). In both sexes, the transcripts with the highest tissue-specificity were found in the brain, and the highest overlap of transcripts between two tissues was found between brain and pituitary. The number of transcripts shared by all tissues was 16,959 in female, while it was 19,334 in male. The numbers of tissue-specific transcripts of brain and pituitary were similar between female and male, while in the gonad, this was approximately 1.4 times higher in testis than in ovary.
Comparative expression profiling between ovary and testis showed that 9746 and 15,804 transcripts were significantly expressed (false discovery rate: FDR < 0.05), and 61 and 33 GO terms were significantly enriched (FDR < 0.05) in ovary and testis, respectively (Additional file 3: Table S1 and S2). Among the top 10 GO terms enriched in the ovary, five germ cell-related terms, namely, positive regulation of acrosome reaction (GO:2000344), acrosin binding (GO:0032190), egg coat formation (GO:0035803), binding of sperm to zona pellucida (GO:0007339), and structural constituent of egg coat (GO:0035804), were found (Table 2). Among the top 10 GO terms enriched in the testis, spermatogenesis (GO:0007283) and six cilium-related terms, namely, ciliary basal body (GO:0036064), cilium movement (GO:0003341), intraciliary transport (GO:0042073), motile cilium (GO:0031514), ciliary transition zone (GO:0035869), and outer dynein arm assembly (GO:0036158) were found (Table 3).
Identification of reproduction-related genes in the BPG axis
In the RNA-seq analysis, three forms of Gnrh (herring-type Gnrh [gnrh1], chicken-II-type Gnrh [gnrh2], and salmon-type Gnrh [gnrh3]) and Gth subunits (fshb, lhb, and cga) were identified with their receptors, such as Gnrh receptors (gnrhrs), Fsh receptor (fshr), and Lh receptor (lhr). Of the gnrhrs identified, the deduced amino acid sequences of some transcripts showed high homology to three types of Gnrhr identified in medaka, that is, Gnrhr1, Gnrhr2, and Gnrhr3 (GenBank accession nos. NP_001098352, NP_001098392, and NP_001098393, respectively). The sequences of genes encoding steroidogenic proteins involved in ovarian steroidogenesis were also identified for cyp11a1, cyp17a1, cyp19a1a, hsd3b, and star. Additionally, eight types of 17β-hydroxysteroid dehydrogenase (hsd17bs; hsd17b3, hsd17b4, hsd17b7, hsd17b8, hsd17b10, hsd17b12a, hsd17b12b, and hsd17b14) were found. In addition to the above 24 genes, genes involved in the system of kisspeptin (Kiss)/Kiss receptor (kiss1, kiss2, and kiss1r) and thyroid-stimulating hormone (Tsh)/Tsh receptor (tshb and tshr), and genes of the other cytochromes P450 (cyp1a1, cyp1b1, and cyp19a1b) were selected to generate clustering of the differentially expressed reproduction-related genes in the BPG axis. For four receptor genes (fshr, lhr, gnrhr1, and gnrhr3), the 5′ and 3′ partial ORFs were found from the identified transcripts, while for other genes, the complete ORFs were found. A total of 36 transcripts (28 complete and 8 partial ORFs) were chosen to profile the abundance of reproduction-related transcripts in each tissue and they formed eight clusters (G1 to G8; Fig. 2). In the G1 cluster, the highest expression of fshb, lhb, and cga was observed in the pituitary. In the G2 cluster, six types of hsd17bs (hsd17b4, hsd17b7, hsd17b8, hsd17b10, hsd17b12a, and hsd17b12b) and hsd3b showed high expression in all tissues. In the G3 cluster, tshb, tshr, and cyp1b1 showed high expression in the pituitary, whereas their expression was low in the gonad. In the G4 cluster, fshr, lhr, cyp11a1, cyp17a1, hsd17b14, and star were expressed mainly in the gonad. In the G5 and G6 clusters, kiss1r, gnrh1, gnrh2, gnrhr2, and gnrhr3 were expressed mainly in the brain and pituitary. In the G7 cluster, cyp1a1 and cyp19a1a were expressed mainly in the brain and gonad, and high expression of cyp1a1 was observed in the ovary. In the G8 cluster, kiss1, kiss2, gnrh3, gnrhr1, and cyp19a1b were mainly expressed in the brain.
The transcription levels of hsd17bs were further compared between ovary and testis (Fig. 3). The expression of hsd17b7 was higher in testis than in ovary. In contrast, the expression of hsd17b8, hsd17b10, hsd17b12a, and hsd17b12b was higher in ovary than in testis. hsd17b3, hsd17b4, and hsd17b14 showed no difference in expression between ovary and gonad.
Sequence alignments of ovarian steroidogenesis-related genes
Genes that play roles in the ovarian steroidogenic pathway were identified based on the KEGG database. Of these, cDNAs encoding Gth subunits (fshb, lhb, and cga), Gthrs (fshr and lhr), four steroidogenic enzymes (cyp11a1, cyp17a1, cyp19a1a, and hsd3b), and Star (star) were cloned and sequenced using a PCR-based strategy (Table 4; Additional file 4).
Multiple clones exhibiting different amino acid sequences for Fshb were obtained from ten pituitaries, where some showed a shift of an N-linked glycosylation site (see detailed explanation in Additional file 5). For real-time PCR, gene-specific primers for fshb were designed to include nucleotide sequences conserved in all cDNA sequences determined. In contrast, the deduced amino acid sequences of Lhb and Cga were highly conserved among different individuals. The deduced amino acid sequences of Japanese sardine Fshb, Lhb, and Cga are 57.3, 79.9, and 80.7% similarity to those of Atlantic herring, respectively. The Japanese sardine Fshb possesses 48.9 and 38.3% sequence similarity with zebrafish and medaka Fshb, but Lhb and Cga share 58–75% sequence similarity between them (Table 4; Additional file 4: Figure S2).
The deduced amino acid sequences of both Japanese sardine Fshr and Lhr are 81.8 and 76.1% similarity to those of Atlantic herring, respectively, while they possess 62–77% sequence similarities to those of zebrafish and medaka (Table 4; Additional file 4: Figure S3 and S4).
The deduced amino acid sequences of Japanese sardine steroidogenic proteins encoded by cyp11a1, cyp17a1, cyp19a1a, hsd3b, and star share high sequence similarities with those of Atlantic herring (92–94%), while they possess 79–90% sequence similarities to those of zebrafish and medaka (Table 4; Additional file 4: Figure S5–S9).
Changes in ovarian development
To analyze the expression profiles of ovarian steroidogenesis-related genes using real-time PCR, female Japanese sardines at age 2 were sampled in September, October, December, and January, and fish at age 3 were sampled in September (Table 5). The body size was 157.0 ± 3.3 mm in BL and 41.7 ± 3.2 g in BW in September at age 2, whereas it was 184.5 ± 3.9 mm in BL and 78.8 ± 3.4 g in BW in September at age 3. The gonadosomatic index (GSI) was 0.44 ± 0.03 in September at age 2, and it increased to 2.63 ± 0.36 in January, but decreased to 0.48 ± 0.03 in September at age 3.
In September, all fish at ages 2 and 3 had immature ovaries, occupied by perinucleolus-stage oocytes (Fig. 4a). In October, most fish had ovaries containing cortical alveolus-stage oocytes, indicating the progress of previtellogenic growth (Fig. 4b). In December, some females were at the middle phase of vitellogenesis (Fig. 4c), while in January, most females were at the late phase of vitellogenesis (Fig. 4d).
Changes in ovarian steroidogenesis-related gene expressions
The expression of fshb was highest in December at age 2, while it was lowest in September at age 3 (Fig. 5). The expression of lhb was lowest in September at age 2, gradually increased, peaked in January and declined in September at age 3. The expression of cga did not change significantly throughout the experimental period.
The expression of fshr was lowest between October and December at age 2, increased in January, but decreased in September at age 3 (Fig. 6). A significant increase in lhr expression was found from December to January at age 2. The expression of cyp11a1, cyp17a1, and cyp19a1a was low in September at age 2 and then increased, peaked in January but decreased in September at age 3. The expression of hsd3b and star did not change significantly throughout the experimental period, although the level of star tended to be higher in January.
Changes in serum steroid hormone levels
Serum androstenedione (AD) levels were low in September and December at ages 2 and 3 (72.2 ± 9.5 to 109.3 ± 11.9 pg/ml), while they were the highest in January at age 2 (351.4 ± 17.2 pg/ml) (Fig. 7). Serum E2 levels were low in September at age 2 (82.2 ± 7.4 pg/ml), gradually increased from October, peaked in January (812.1 ± 100.4 pg/ml), and then declined in September at age 3 (169.3 ± 48.2 pg/ml). Serum testosterone (T) levels remained low throughout the experimental period (below 46.2 ± 11.1 pg/ml). Serum estrone (E1) levels were below the assay detection limit (< 15 pg/ml) in the serum of all fish.
Japanese sardine transcriptome in the BPG axis
In the present study, a total of 115,173 ORFs were predicted from transcripts in the BPG axis of Japanese sardine. Of these, about 40% of protein sequences have similarity to Atlantic herring, zebrafish, and medaka (Fig. 1a). In addition, the distribution of percent identity of amino acids indicated that the sardine protein sequence is more identical to Atlantic herring than to zebrafish and medaka (Fig. 1b). In the genomes of Atlantic herring, zebrafish, and medaka, 43,000–53,000 proteins have been annotated [13, 26,27,28]. In the present study, the use of these three species as a reference detected about half of the proteins, for example, 24,037 proteins of Atlantic herring out of 47,201 proteins of the sardine (Fig. 1a). Among those proteins, we succeeded in obtaining sequence information of many key genes involved in the BPG axis of Japanese sardine.
A Venn diagram indicated that the numbers of tissue-specific transcripts in the brain and pituitary were similar between female and male, while in the gonad this number was higher in testis than in ovary (Fig. 1c). Similar trends for genes differentially expressed between ovary and testis have been reported in some fish species, such as bluehead wrasse (Thalassoma bifasciatum) , Nile tilapia (Oreochromis niloticus) , and sharpsnout seabream (Diplodus puntazzo) , where many more genes were expressed in testis than in ovary. Male-biased gene expression may be a general characteristic feature of fish gonad. Additionally, in the wrasse, such sex-biased gene expression in gonad was more potent than that in brain , similar to the tissue-specific expression patterns observed in the sardine.
The genes differentially expressed between ovary and testis were further characterized by the enriched GO terms. The majority of top GO terms enriched in the ovary were related to oogenesis and sperm–egg recognition, while those in the testis were related to spermatogenesis and cilium function, which may be linked to formation of the sperm flagellum. Thus, many genes differentially expressed between ovary and testis were directly linked to gamete production in each gonad. These results will provide us with numerous candidate genes potentially responsible for regulating gametogenesis in future study of the sardine.
Differential expression of reproduction-related genes in the BPG axis
Many genes encoding the major regulators in the BPG axis were identified in the Japanese sardine, and their transcription levels showed tissue- and sex-specific expression patterns (Fig. 2).
As in other vertebrates, teleost Gths are synthesized in gonadotropic cells in the pituitary, in which transcripts of Gth subunits are abundant . Indeed, the highest expression of fshb, lhb, and cga was detected in the sardine pituitary. We also found high expression of tshb in the sardine pituitary. Tsh is composed of common Cga and hormone-specific β-subunit (Tshb) and produced in the pituitary like Gth, but it exhibits various physiological functions, such as in growth and metamorphosis including reproduction . Teleost Tshr has been detected in various tissues, such as gonad, brain, kidney, and gills [33, 34]. In the present study, high expression of sardine tshr was found in the pituitary.
As in other teleosts, transcripts encoding three forms of Gnrh (gnrh1, gnrh2, and gnrh3) were identified within the present annotated sequences, where Gnrh1 is the herring type that is found in fish belonging to the order Clupeiformes . In fish brain, different expression of the three gnrhs has been found in a variety of regions, while Gnrh1 is suggested to be the main regulator of Gth secretion because Gnrh1-expressing neurons innervate into the pituitary [16, 36]. In the present study, gnrh1 and gnrh2 were mainly expressed in the brain and pituitary, while gnrh3 was expressed mainly in the brain. Numerous Gnrhrs have been identified in teleost, but their classification has not been fully established . In medaka, three types of gnrhr (gnrhr1, gnrhr2, and gnrhr3) have been characterized in detail [38, 39]. We also found these three types in the Japanese sardine, showing that the gnrhr2 and gnrhr3 expression was mainly detected in the brain and pituitary, whereas gnrhr1 was in the brain. This suggests that, in the sardine, Gnrhr2 and Gnrhr3 are plausible candidates for mediating signal transduction of Gnrh-stimulated Gth response. While the Kiss/Kissr system has been proposed to be involved in the neuroendocrine control of the BPG axis in fish as recognized in mammals, the physiological roles in fish reproduction are still not well defined [40,41,42]. In the present study, the expression of kiss1 and kiss2 was mainly detected in the sardine brain.
Five genes (fshr, lhr, cyp11a1, cyp17a1, and star) expressed mainly in the ovary are well known to play major roles in the ovarian steroidogenic pathway . Aromatase also plays a central role in the biosynthesis of estrogen, and teleosts possess two types of aromatase gene: ovarian-type cyp19a1a and brain-type cyp19a1b . In the present study, cyp19a1a was expressed in the brain and gonad, whereas cyp19a1b was mainly expressed in the brain, suggesting that the sardine ovarian aromatase is encoded by cyp19a1a, as in other fishes. Hsd3b is another critical enzyme in ovarian steroidogenesis. In the present study, high expression of hsd3b was found in the gonad, while it was also expressed in the brain and pituitary. A study in stinging catfish (Heteropneustes fossilis) showed that hsd3b is widely expressed in the brain, suggesting its role in neurosteroidogenesis . Cytochrome P450 1A1 (cyp1a1) and 1B1 (cyp1b1) play roles in E2 metabolism in mammals, while cyp1a1 is suggested to be primarily responsible for this function in zebrafish . In the present study, cyp1a1 expression was mainly found in the brain and gonad, especially in the ovary, whereas cyp1b1 was expressed mainly in the brain and pituitary, suggesting a specific role of cyp1a1 in the ovary like in zebrafish. Finally, we also found high expression of six types of hsd17b (hsd17b4, hsd17b7, hsd17b8, hsd17b10, hsd17b12a, and hsd17b12b) in all three tissues. We further compared the transcription levels of hsd17bs, focusing on their sex-specific expression.
Hsd17b catalyzes the interconversion of 17-ketosteroids to 17β-hydroxysteroids, such as E1 to E2, and a total of 15 Hsd17b family genes have been identified in mammals [46, 47]. In fish, several different hsd17b transcripts have been isolated, showing differential tissue-specific expression and catalytic activities with various types of Hsd17b [48,49,50]. However, the specific function of each Hsd17b in ovarian steroidogenesis has not been clearly confirmed. The present study identified eight types of Hsd17b in the Japanese sardine. By comparing transcription levels between ovary and testis, higher expression of hsd17b8, hsd1710, hsd17b12a, and hsd17b12b was found in the ovary, whereas higher expression of hsd17b7 was found in the testis (Fig. 3). A study of the Japanese eel (Anguilla japonica) revealed that Hsd17b12a is one of the enzymes mediating the biosynthesis of 11-ketotestosterone, which is the main androgen in fish spermatogenesis . Therefore, it is interesting that hsd17b12a was found to be most highly expressed in the sardine ovary. In mammals, Hsd17b1 expressed in granulosa cells play a central role in estrogen biosynthesis, and Hsd17b7 also catalyzes the conversion of E1 to E2 in the ovary and placenta . However, the sequence of the hsd17b1 was not obtained in our data, and hsd17b7 was found but showed higher expression in testes than in ovary. Accordingly, while our results suggest sex-specific differences in the function of each Hsd17b in the sardine, it is still difficult to suggest which type would have a role in estrogen biosynthesis in this species. The present sequence data will allow future studies on the physiological roles of Hsd17bs in sardine reproduction, which should be useful for understanding the gonadal steroidogenic pathway unique to teleost fish.
Molecular cloning of ovarian steroidogenesis-related genes
A total of 10 genes that have critical roles in ovarian steroidogenesis were selected, and the cDNAs were cloned and sequenced. The deduced amino acid sequences of Japanese sardine Cga and Lhb share approximately 80% sequence similarity with their homologs in Atlantic herring, whereas Fshb shares only 57% similarity. The sequences and structures of Cga and Lhb are generally more conserved than those of Fshb in teleosts . Likewise, the sequences of Lhr are more conserved than those of Fshr among vertebrate species . However, the sequence similarities of sardine Lhr to those in zebrafish and medaka were lower than those of Fshr. In contrast, the sequences of steroidogenic proteins are highly conserved among teleosts . In accordance with this, we found high sequence similarities (> 80%) of steroidogenic proteins encoded by cyp11a1, cyp17a1, cyp19a1a, hsd3b, and star between sardine and other fish species.
For sardine Fshb, multiple fshb clones were isolated from transcripts of ten pituitaries (Additional file 5). This suggests that the sardine Fshb may have many amino acid sequence variants between individuals. In contrast, the existence of isoforms of Gth subunits has been reported in some fish species such as zebrafish and salmonids and this seems to be derived from tetraploidy of the fish species . Genomic analysis of Atlantic salmon (Salmo salar) revealed that its genome exhibits largely tetraploid due to a genome duplication, showing a high repeat content (58–60%) . Therefore, it is also important to consider the possibility of the existence of Fshb paralogues in the sardine. Another sequencing approach targeting the isoforms and variants might be required to confirm this association.
Expression profiles of ovarian steroidogenesis-related genes
The present study and other previous studies [55, 56] showed that serum E2 levels in the Japanese sardine gradually increased during ovarian development. This was accompanied by increased serum vitellogenin levels . In the present study, we focused on the ovarian steroidogenic pathway and analyzed the expression profiles of ovarian steroidogenesis-related genes in relation to the vitellogenic process (Fig. 8).
The physiological role of Fsh has been studied in detail in salmonids, showing that Fsh stimulates the production of E2 in follicle cells and induces vitellogenin uptake [15, 57, 58]. Conversely, the major physiological role of Lh is well recognized in many teleost species as a main regulator of the production of maturation-inducing steroid, which initiates the final phase of oocyte development, that is oocyte maturation and ovulation . In the present study, the peak of pituitary fshb expression was found during the middle phase of vitellogenesis, whereas lhb expression peaked at the completion of vitellogenesis. This difference may reflect a shift of the activation of Gth synthesis from Fsh to Lh in association with a shift from a vitellogenic process to oocyte maturation. Pituitary cga expression did not change throughout ovarian development, suggesting that cga expression is kept constant regardless of the activation of Gth synthesis.
The present study demonstrated that, in the sardine ovary, the expression of fshr and lhr was upregulated during the late phase of vitellogenesis. The size-dependent separation of ovarian follicles (oocytes surrounded by follicle cell layers) revealed that fshr expression was elevated during vitellogenesis, peaking at the completion of vitellogenesis, and then decreased rapidly after the initiation of oocyte maturation . In contrast, lhr expression is upregulated at the completion of vitellogenesis and maintained at a high level during the early phase of oocyte maturation. Thus, in individual ovarian follicles, elevated fshr expression during vitellogenesis is switched to rapid expression of lhr by the completion of vitellogenesis, preparing for the initiation of oocyte maturation. Although we did not analyze the oocyte size-dependent gene expression, the significant increases in fshr and lhr expression in the ovary at the completion of vitellogenesis may reflect the active egg production in preparation for spawning.
In salmonids, the ovarian gene expression of steroidogenic proteins (cyp11a1, cyp17a1, cyp19a1a, hsd3b, and star) is upregulated during vitellogenesis, in correlation with the production of E2 [61, 62]. In the present study, the expression of cyp11a1, cyp17a1, and cyp19a1a was lowest in immature ovary, gradually increased during ovarian development, and peaked at the completion of vitellogenesis. However, no significant changes in the expression of star and hsd3b occur, although star expression tended to be higher in the late-vitellogenic ovary than in immature ovary. In salmonids, Fsh upregulated cyp11a1, cyp17a1, and cyp19a1a expression in the previtellogenic ovarian follicles, but it was more potent in inducing star and hsd3b expression [63, 64]. Such differences in the Fsh response to gene expression may lead to differential expression profiles found in sardine steroidogenic genes.
The in vitro cultivation of ovarian follicles with radiolabeled steroid precursors has revealed that E2 is synthesized via T in ovarian follicles of chub mackerel (Scomber japonicus) and medaka, whereas E2 is synthesized via E1 in bambooleaf wrasse (Pseudolabrus japonicus) and red seabream (Pagrus major) [65,66,67,68]. Thus, there are species differences in the main substrate precursor of E2. In the present study, the serum levels of AD, which is the substrate precursor of T and E1, were elevated during ovarian development. Since aromatase catalyzes the conversion of T to E2 and AD to E1, the upregulation of cyp19a1a during ovarian development suggests that the catalytic activity of aromatase on E2 conversion from T or E1 conversion from AD may be activated. However, serum levels of T remained low throughout ovarian development, consistent with previous reports on sardine [55, 56]. We also found that E1 was not detected in the serum of female sardines. Hence, there is a need for further research on the E2 synthetic pathway of sardine ovary. Here, the physiological characterization of each type of Hsd17 will be important to consider since Hsd17b catalyzes the activation of AD to T and E1 to E2.
In the present work, we have identified reproduction-related genes from the BPG axis of Japanese sardine using NGS techniques. Genes encoding major regulators in the BPG axis were characterized in terms of their tissue- and sex-specific expression. Here, four of eight Hsd17b genes (hsd17b8, hsd17b10, hsd17b12a, and hsd17b12b) identified showed higher expression in ovary than in testis. This is of particular interest since the function of Hsd17bs in ovarian steroidogenesis is not well established. Despite the ecological and commercial importance worldwide, transcriptome information has rarely been available in the physiological study of clupeoid fish. The present transcriptomic sequences will provide a valuable resource for studying the molecular biology of clupeoid fish.
In the present study, a detailed molecular picture of the ovarian steroidogenic pathway was indicated from our transcript data of endocrine factors (Fig. 8). The information provided a good foundation for understanding the molecular mechanisms in the physiological process of sardine ovarian development. External factors such as temperature and feeding, and internal factors such as age and growth are known to influence gene expression in the BPG axis [14, 69]. Additional studies on the particular endocrine factors and the responses to changes in external and internal factors should provide further insight into the endocrine regulation of ovarian development in this species. Studies on wild sardine have demonstrated that characteristics of ovarian development such as fecundity and timing of maturity varied in association with changes in population size [23, 24]. Therefore, the physiological regulation of ovarian development at an individual level may lead to a new understanding of mechanisms associated with population fluctuation in the sardine.
Lastly, RNA-seq analysis of the reproduction-related transcriptome can serve as a foundation for studying the reproductive physiology of fisheries’ target species. We performed captive experiments since these often provide useful information for understanding the reproductive characteristics of fish with high reproducibility . Although specific traits in wild fish as compared to captive fish should also be considered, endocrine regulation of gonadal development, as well as molecular properties, might be conserved within each species. A combination of molecular analysis and field-research-dependent biological data would open up interesting avenues for the study of reproductive biology and stock assessment.
Fish and sample collection
The Japanese sardine used in this study were commercially caught in the wild as juveniles (0+ years old). About 400 fish were transferred to indoor concrete tanks at Hakatajima Field Station, Fisheries Technology Institute (formerly Hakatajima Laboratory, National Research Institute of Fisheries and Environment of Inland Sea), Japan Fisheries Research and Education Agency, Kinoura, Imabari, Ehime, Japan. The fish were then maintained for 1 to 2 years. Throughout the maintenance of the general farming system, the fish were regularly fed with commercial dry pellets. Of those reared, a total of 31 fish was used in the present study. On 18 February 2017, three females and three males were sampled at age 3. This period was chosen since the expression of many genes related to the BPG axis is activated during active gonadal development [14, 61]. For each fish, brain, pituitary, and a portion of gonad were frozen and stored at − 80 °C for RNA-seq analysis. Other fish were sampled on 17 October and 4 December, 2017, and 16 January and 18 September, 2018. Five to six females at age 2 were taken in October, December, January, and September, while four females at age 3 were taken in September. The BL, BW, and gonadal weight (GW) of each specimen were measured. Blood samples were centrifuged (3500×g, 15 min, 4 °C), and the serum was collected and stored at − 30 °C for steroid hormone measurements. A portion of each ovary was fixed in Bouin’s solution for histological observation. The rest of the ovary and pituitary were immersed in RNAlater (Qiagen, Hilden, Germany) for 24 h and then stored at − 80 °C for molecular cloning and quantitative real-time PCR analysis. The GSI was calculated as follows: GSI = (GW / BW) × 100. In total, 25 fish were used for the quantitative analysis of gene expression and steroid hormones (Table 5). The sample size of n = 4–6 per group was determined according to previous studies on fish [61, 71]. For each sampling, the fish were euthanized by rapid decapitation. The treatment of fish adhered to the guidelines for animal experiments set by the National Research Institute of Fisheries and Environment of Inland Sea.
A portion of each gonad fixed in Bouin’s solution was dehydrated through a graded series of ethanol. Dehydrated tissues were embedded in paraffin. Sections were cut (5 μm) and stained with hematoxylin–eosin.
RNA extraction and sequencing
Total RNA was extracted from brain, pituitary, and gonad (three fish for each sex) that had been stored at − 80 °C using a TRIzol Plus RNA Purification Kit (Thermo Fisher Scientific, Waltham, MA, USA) with PureLink DNase (Thermo Fisher Scientific), in accordance with the manufacturer’s instructions. The concentration and quality of purified total RNA were determined using a Qubit RNA HS Assay Kit (Thermo Fisher Scientific) and a Bioanalyzer 2100 with RNA 6000 Nano Kit (Agilent Technologies, Inc., Santa Clara, CA, USA).
The cDNA libraries were constructed using the total RNA (1 μg) and SuperScript III reverse transcriptase (Thermo Fisher Scientific), in accordance with the TruSeq RNA Sample Prep ver. 2 (LS) protocol (Illumina, Inc., San Diego, CA, USA). The cDNA libraries were sequenced into 150-bp paired-end reads by NextSeq500 (Illumina, Inc.) at the Bioinformatics and Biosciences Division of the Fisheries Resources Institute (formerly Research Center for Bioinformatics and Biosciences of the National Research Institute of Fisheries Science), Yokohama, Japan.
Transcriptome assembly and annotation
The assembly was performed as previously described . Briefly, the sequence trimming was performed using Trimmomatic , and the quality confirmation was performed using FastQC . The remaining paired-end reads were assembled using Trinity . From the assembled sequences, ORFs of more than 300 bp were extracted and translated into protein sequences using TransDecoder . ORFs of 95% similar protein sequences were clustered using the CD-HIT program  to remove redundant protein sequences.
The proteins were annotated based on their homology to protein sequences of the Atlantic herring, zebrafish, and medaka, using the BLASTP program, with a threshold e-value of <1e–05. A homology search was also performed using DIAMOND software  and the NCBI nr database with a threshold e-value of <1e–05. The cumulative relative frequency for amino acid identity against each database was plotted using R software ver. 3.3.1 . The GO and EC numbers, which are shared with the accession numbers used in zebrafish and medaka, were assigned from the best hits of the BLASTP results against these two species. The retrieved accession numbers of both zebrafish and medaka were converted to UniProt accession numbers at the UniProt website (https://www.uniprot.org/) and then GO terms were retrieved from the UniProt sequence data. The retrieved GO terms were converted to GO-slim which is cut-down version of GO  by OmicsBox software (; https://www.biobam.com/omicsbox). Based on the assigned EC numbers and annotated descriptions in Japanese sardine, the protein sequences related to reproduction were obtained referring to the Gnrh signaling pathway (Pathway ID: 04912) and steroid hormone biosynthesis (Pathway ID: 00140) in Kyoto Encyclopedia of Genes and Genomes (KEGG).
Differential expression analysis
The transcription level of each gene was determined from the number of mapped reads. The mapping of reads to each transcript and the counting were performed using Bowtie2  and RSEM , respectively. Transcription levels were normalized among the libraries using the trimmed mean of M-values method  and statistically compared with the edgeR package  in R software ver. 3.3.1. Fisher’s exact test was performed in OmicsBox software to estimate whether the differentially expressed genes were particularly associated with specific GO categories when compared with the background genes. P-values were adjusted by FDR , and the FDR of < 0.05 was chosen as the threshold for significance. The normalized fragments per kilobase of transcript per million mapped read (FPKM) of genes was used for generating Venn diagrams and a heatmap with R software ver. 3.3.1.
The complete nucleotide sequences for cDNAs encoding Japanese sardine Gth subunits (fshb, lhb, and cga), Gthrs (fshr and lhr), three cytochrome P450 enzymes (cyp11a1, cyp17a1, cyp19a1a), Hsd3b (hsd3b), and Star (star) were determined using basic molecular cloning techniques as previously described . fshb, lhb, and cga were cloned from the pituitary, while fshr, lhr, cyp11a1, cyp17a1, cyp19a1a, hsd3b, and star were cloned from the ovary. Total RNA was extracted from the pituitary using an RNeasy Lipid Tissue Mini Kit (Qiagen), while total RNA was extracted from the ovary using an RNeasy Mini Kit (Qiagen). RNA extracts were treated with DNase using the Turbo DNA-free Kit (Thermo Fisher Scientific). First-strand cDNA was synthesized using a Superscript III reverse transcriptase (Invitrogen, Carlsbad, CA, USA) with random hexamer primers. For fshb, lhb, and cga, the cDNAs from four pituitaries were subjected to PCR amplification. Since several transcripts encoding different amino acid sequences of Fshb were isolated, the cDNAs from the other six pituitaries were further used for cloning of fshb. For other genes, the cDNAs from two individuals were subjected to PCR amplification.
Based on the sequences of contigs, primers were designed in the 5′ and 3′ untranslated region (5′ and 3′ UTR) to obtain the clones containing full ORFs (Additional file 6: Table S3). For fshb, lhb, cga, cyp11a1, cyp17a1, cyp19a1a, hsd3b, and star, these contigs contain the complete ORF with the 5′ and 3′ UTR. For fshr and lhr, only contigs containing 5′ and 3′ ORFs were found. GenBank accession numbers of the contigs used for designing primers are listed in Additional file 7: Table S6.
The similarity of protein-coding sequences was calculated through pairwise alignment using BioEdit software.
Expression levels of pituitary and ovarian genes were determined by real-time PCR using a Thermal Cycler Dice Real Time System (Takara, Shiga, Japan). Gene-specific primers were designed for fshb, lhb, cga, fshr, lhr, cyp11a1, cyp17a1, cyp19a1a, hsd3b, and star (Additional file 6: Table S4 and S5). Although several isoforms of Fshb were found in the molecular cloning, primers for real-time PCR were designed to count all fshb isoforms because it is currently unclear whether these isoforms exhibit different functions.
Using the same methods as described in the section on cDNA cloning, total RNAs were extracted from the pituitary and from a piece of ovary, and template cDNAs were synthesized after DNase treatments. For each PCR reaction, 20 μl of PCR reaction mix was prepared, which contained 10 μl of TB Green Premix Ex Taq II (Takara), 0.6 μM of forward and reverse primers, and template cDNA corresponding to 0.3 ng (pituitary) or 50 ng (ovary) of total RNA. Real-time PCR runs were performed in 96-well plates in duplicate with no template controls. The thermal cycling conditions were 30 s at 95 °C, and 40 cycles of 5 s at 95 °C and 30 s at 60 °C. For each PCR, a standard curve from serial dilutions of purified PCR fragments containing a partial sequence of a target gene was constructed for quantification.
Steroid hormone measurements
Serum AD was measured using time-resolved fluoroimmunoassay (TR-FIA), as previously described . Serum E1, E2, and T were measured using the enzyme-linked immunosorbent assay (ELISA) kit (Abnova, CA, USA for E1; Cayman Chemical, MI, USA for E2 and T), in accordance with the manufacturer’s instructions. Before measurements of AD, E2, and T, the serum steroids were extracted into diethyl ether, which was then evaporated, and the remaining steroids were dissolved in suitable assay buffers for each assay. In the TR-FIA, AD-3-carboxymethyloxime-bovine serum albumin (AD-3-CMO-BSA) and anti-AD-3-CMO-BSA were used as coating antigen and primary antibody, respectively. All samples were assayed in duplicate for both TR-FIA and ELISA.
Data are presented as mean ± standard error of the mean (SEM). GraphPad Prism version 8.00 for Windows (GraphPad Software, La Jolla, CA, USA) was used for statistical evaluation. Multiple paired comparisons of FPKM values between ovary and testis were made using the Holm–Sidak method (α = 0.05). Results of gene expressions from real-time PCR were analyzed by one-way ANOVA followed by Tukey’s multiple comparisons test.
Availability of data and materials
The raw RNA sequencing reads from the Japanese sardine (Sardinops melanostictus) brain, pituitary, and gonad are available in DRA division of DDBJ under accession number DRA010129 (http://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA010129). The assembled contigs are also available in TSA division of DDBJ under accession numbers ICPT01000001–ICPT01290119 (http://getentry.ddbj.nig.ac.jp/top-e.html). The results of Protein Basic Local Alignment Search Tool are shown in Additional file 1. The GO ID linked to enriched GO terms are shown in Additional file 3. GenBank accession numbers for the sequences determined in this work are LC545596–LC545610. The reference genomes were obtained from the NCBI database (Atlantic herring [Clupea harengus], GCF_900700415.1; zebrafish [Danio rerio], GCF_000002035.6; medaka [Oryzias latipes], GCF_002234675.1). GenBank accession numbers used for sequence similarity calculation are included in Additional file 4. GenBank accession numbers of contigs used for clustering analysis are shown in Additional file 7.
Protein Basic Local Alignment Search Tool
- BPG axis:
Bovine serum albumin
Glycoprotein hormones, alpha polypeptide
Cytochrome P450, family 1, subfamily A, polypeptide 1
Cytochrome P450, family 1, subfamily B, polypeptide 1
Cytochrome P450, family 11, subfamily A, polypeptide 1
Cytochrome P450, family 17, subfamily A, polypeptide 1
Cytochrome P450, family 19, subfamily A, polypeptide 1a
Cytochrome P450, family 19, subfamily A, polypeptide 1b
DNA Data Bank of Japan
DDBJ Sequence Read Archive
Enzyme-linked immunosorbent assay
False discovery rate
Fragments per kilobase of transcript per million mapped read
Follicle-stimulating hormone subunit beta
Follicle-stimulating hormone receptor
Gonadotropin-releasing hormone receptor
Kyoto Encyclopedia of Genes and Genomes
Luteinizing hormone subunit beta
Luteinizing hormone receptor
National Center for Biotechnology Information
- NCBI nr database:
NCBI non-redundant protein database
Open reading frame
Standard error of the mean
Steroidogenic acute regulatory protein
Transcriptome shotgun assembly
Thyroid-stimulating hormone receptor
Anitha A, Gupta YR, Deepa S, Ningappa M, Rajanna KB, Senthilkumaran B. Gonadal transcriptome analysis of the common carp, Cyprinus carpio: identification of differentially expressed genes and SSRs. Gen Comp Endocrinol. 2019;279:67–77.
Bar I, Cummins S, Elizur A. Transcriptome analysis reveals differentially expressed genes associated with germ cell and gonad development in the southern bluefin tuna (Thunnus maccoyii). BMC Genomics. 2016;17:217.
Monson C, Forsgren K, Goetz G, Harding L, Swanson P, Young G. A teleost androgen promotes development of primary ovarian follicles in coho salmon and rapidly alters the ovarian transcriptome. Biol Reprod. 2017;97:731–45.
Rozenfeld C, García-Carpintero V, Pérez L, Gallego V, Herranz-Jusdado JG, Tveiten H, Johnsen HK, Fontaine R, Weltzien FA, Cañizares J, Asturiano JF, Peñaranda DS. Cold seawater induces early sexual developmental stages in the BPG axis of European eel males. BMC Genomics. 2019;20:597.
Wang ZP, Wang D, Wang CL, Xie WJ, Zhu YF, Chen XW. Transcriptome characterization of HPG axis from Chinese sea perch Lateolabrax maculatus. J Fish Biol. 2017;91:1407–18.
Barth JMI, Berg PR, Jonsson PR, Bonanomi S, Corell H, Hemmer-Hansen J, Jakobsen KS, Johannesson K, Jorde PE, Knutsen H, Moksnes PO, Star B, Stenseth NC, Svedäng H, Jentoft S, André C. Genome architecture enables local adaptation of Atlantic cod despite high connectivity. Mol Ecol. 2017;26:4452–66.
Suda A, Nishiki I, Iwasaki Y, Matsuura A, Akita T, Suzuki N, Fujiwara A. Improvement of the Pacific bluefin tuna (Thunnus orientalis) reference genome and development of male-specific DNA markers. Sci Rep. 2019;9:14450.
Gioacchini G, Marisaldi L, Basili D, Candelma M, Pignalosa P, Cigliano RA, Sanseverino W, Hardiman G, Carnevali O. A de novo transcriptome assembly approach elucidates the dynamics of ovarian maturation in the swordfish (Xiphias gladius). Sci Rep. 2019;9:7375.
Chavez FP, Ryan J, Lluch-Cota SE, Ñiquen CM. From anchovies to sardines and back: multidecadal change in the Pacific Ocean. Science. 2003;299:217–21. https://doi.org/10.1126/science.1075880.
Checkley DM, Alheit J, Oozaki Y, Roy C. Climate change and small pelagic fish. Cambridge: Cambridge University Press; 2009.
Takasuka A. Biological mechanisms underlying climate impacts on population dynamics of small pelagic fish. In: Aoki I, Yamakawa T, Takasuka A, editors. Fish population dynamics, monitoring, and management. Tokyo: Springer; 2018. p. 19–50.
Martinez Barrio A, Lamichhaney S, Fan G, Rafati N, Pettersson M, Zhang H, Dainat J, Ekman D, Höppner M, Jern P, et al. The genetic basis for ecological adaptation of the Atlantic herring revealed by genome sequencing. eLife. 2016;5:e12081.
Sí K, Mikalsen SO, Eí H, Jacobsen JA, Flicek P, Dahl HA. Using long and linked reads to improve an Atlantic herring (Clupea harengus) genome assembly. Sci Rep. 2019;9:17716.
Levavi-Sivan B, Bogerd J, Mañanós EL, Gómez A, Lareyre JJ. Perspectives on fish gonadotropins and their receptors. Gen Comp Endocrinol. 2010;165:412–37.
Swanson P, Dickey JT, Campbell B. Biochemistry and physiology of fish gonadotropins. Fish Physiol Biochem. 2003;28:53–9.
Zohar Y, Muñoz-Cueto JA, Elizur A, Kah O. Neuroendocrinology of reproduction in teleost fish. Gen Comp Endocrinol. 2010;165:438–55.
Somoza GM, Miranda LA, Strobl-Mazzulla P, Guilgur LG. Gonadotropin-releasing hormone (GnRH): from fish to mammalian brains. Cell Mol Neurobiol. 2002;22:589–609.
Pierce JG, Parsons TF. Glycoprotein hormones: structure and function. Annu Rev Biochem. 1981;50:465–95.
Lubzens E, Young G, Bobe J, Cerdà J. Oogenesis in teleosts: how fish eggs are formed. Gen Comp Endocrinol. 2010;165:367–89.
Yan L, Swanson P, Dickhoff WW. A two-receptor model for salmon gonadotropins (GTH I and GTH II). Biol Reprod. 1992;47:418–27.
Takasuka A, Yoneda M, Oozeki Y. Disentangling density-dependent effects on egg production and survival from egg to recruitment in fish. Fish Fish. 2019;20:870–87.
Furuichi S, Yasuda T, Kurota H, Yoda M, Suzuki K, Takahashi M, Fukuwaka M. Disentangling the effects of climate and density-dependent factors on spatiotemporal dynamics of Japanese sardine spawning. Mar Ecol Prog Ser. 2020;633:157–68.
Morimoto H. Relationship between batch fecundity and egg size in Japanese sardine Sardinops melanostictus in Tosa Bay and off Kii Channel, southwestern Japan from 1990 to 1993. Fish Sci. 1998;64:220–7.
Morimoto H. Age and growth of Japanese sardine Sardinops melanostictus in Tosa Bay, southwestern Japan during a period of declining stock size. Fish Sci. 2003;69:745–54.
Takasuka A, Yoneda M, Oozeki Y. Density dependence in total egg production per spawner for marine fish. Fish Fish. 2019;20:125–37.
Howe K, Clark MD, Torroja CF, Torrance J, Berthelot C, Muffato M, Collins JE, Humphray S, McLaren K, Matthews L, et al. The zebrafish reference genome sequence and its relationship to the human genome. Nature. 2013;496:498–503.
Kasahara M, Naruse K, Sasaki S, Nakatani Y, Qu W, Ahsan B, Yamada T, Nagayasu Y, Doi K, Kasai Y, et al. The medaka draft genome and insights into vertebrate genome evolution. Nature. 2007;447:714–9.
Pettersson ME, Rochus CM, Han F, Chen J, Hill J, Wallerman O, Fan G, Hong X, Xu Q, Zhang H, et al. A chromosome-level assembly of the Atlantic herring genome —detection of a supergene and other signals of selection. Genome Res. 2019;29:1919–28.
Liu H, Lamm MS, Rutherford K, Black MA, Godwin JR, Gemmell NJ. Large-scale transcriptome sequencing reveals novel expression patterns for key sex-related genes in a sex-changing fish. Biol Sex Differ. 2015;6:26.
Tao W, Yuan J, Zhou L, Sun L, Sun Y, Yang S, Li M, Zeng S, Huang B, Wang D. Characterization of gonadal transcriptomes from Nile tilapia (Oreochromis niloticus) reveals differentially expressed genes. PLoS One. 2013;8:e63604.
Manousaki T, Tsakogiannis A, Lagnel J, Sarropoulou E, Xiang JZ, Papandroulakis N, Mylonas CC, Tsigenopoulos CS. The sex-specific transcriptome of the hermaphrodite sparid sharpsnout seabream (Diplodus puntazzo). BMC Genomics. 2014;15:655.
MacKenzie DS, Jones RA, Miller TC. Thyrotropin in teleost fish. Gen Comp Endocrinol. 2009;161:83–9.
Goto KR, Kazeto Y, Trant JM. Molecular cloning, characterization and expression of thyroid-stimulating hormone receptor in channel catfish. Gen Comp Endocrinol. 2009;161:313–9.
Ponce M, Infante C, Manchado M. Molecular characterization and gene expression of thyrotropin receptor (TSHR) and a truncated TSHR-like in Senegalese sole. Gen Comp Endocrinol. 2010;168:431–9.
Miccoli1 A, Olivotto I, De Felice A, Leonori I, Carnevali O. Characterization and transcriptional profiles of Engraulis encrasicolus’ GnRH forms. Reproduction. 2016;152:727–739.
Sukhan ZP, Kitano H, Selvaraj S, Yoneda M, Yamaguchi A, Matsuyama M. Identification and distribution of three gonadotropin-releasing hormone (GnRH) isoforms in the brain of a clupeiform fish, Engraulis japonicus. Zool Sci. 2013;30:1081–91.
Lumayno SDP, Ohga H, Selvaraj S, Nyuji M, Yamaguchi A, Matsuyama M. Molecular characterization and functional analysis of pituitary GnRH receptor in a commercial scombroid fish, chub mackerel (Scomber japonicus). Gen Comp Endocrinol. 2017;247:143–51.
Okubo K, Ishii S, Ishida J, Mitani H, Naruse K, Kondo M, Shima A, Tanaka M, Asakawa S, Shimizu N, Aida K. A novel third gonadotropin-releasing hormone receptor in the medaka Oryzias latipes: evolutionary and functional implications. Gene. 2003;314:121–31.
Okubo K, Nagata S, Ko R, Kataoka H, Yoshiura Y, Mitani H, Kondo M, Naruse K, Shima K, Aida K. Identification and characterization of two distinct GnRH receptor subtypes in a teleost, the medaka Oryzias latipes. Endocrinology. 2001;142:4729–39.
Alvarado MV, Servili A, Molís G, Gueguen MM, Carrillo M, Kah O, Felip A. Actions of sex steroids on kisspeptin expression and other reproduction-related genes in the brain of the teleost fish European sea bass. J Exp Biol. 2016;219:3353–65.
Hodne K, Weltzien FA, Oka Y, Okubo K. Expression and putative function of kisspeptins and their receptors during early development in medaka. Endocrinology. 2013;154:3437–46.
Ohga H, Selvaraj S, Matsuyama M. The roles of kisspeptin system in the reproductive physiology of fish with special reference to chub mackerel studies as main axis. Front Endocrinol. 2018;9:147.
Diotel N, Page YL, Mouriec K, Tong SK, Pellegrini E, Vaillant C, Anglade I, Brion F, Pakdel F, Chung B, Kah O. Aromatase in the brain of teleost fish: expression, regulation and putative functions. Front Neuroendocrinol. 2010;31:172–92.
Mishra S, Chaube R. Distribution and localization of 3β-hydroxysteroid dehydrogenase (3β-HSD) in the brain and its regions of the catfish Heteropneustes fossilis. Gen Comp Endocrinol. 2017;241:80–8.
Scornaienchi ML, Thornton C, Willett KL, Wilson JY. Cytochrome P450-mediated 17β-estradiol metabolism in zebrafish (Danio rerio). J Endocrinol. 2010;206:317–25.
Labrie F, Luu-The V, Labrie C, Bélanger A, Simard J, Lin SX, Pelletier G. Endocrine and intracrine sources of androgens in women: inhibition of breast cancer and other roles of androgens and their precursor dehydroepiandrosterone. Endocr Rev. 2003;24:152–82.
Luu-The V, Labrie F. The intracrine sex steroid biosynthesis pathways. Prog Brain Res. 2010;181:177–92.
Mindnich R, Adamski J. Zebrafish 17beta-hydroxysteroid dehydrogenases: an evolutionary perspective. Mol Cell Endocrinol. 2009;301:20–6.
Suzuki H, Ozaki Y, Ijiri S, Gen K, Kazeto Y. 17β-Hydroxysteroid dehydrogenase type 12a responsible for testicular 11-ketotestosterone synthesis in the Japanese eel, Anguilla japonica. J Steroid Biochem Mol Biol. 2020;198:105550.
Zou C, Wang L, Zou Y, Wu Z, Wang W, Liang S, Wang L, You F. Characteristics and sex dimorphism of 17β-hydroxysteroid dehydrogenase family genes in the olive flounder Paralichthys olivaceus. J Steroid Biochem Mol Biol. 2020;199:105597.
Hilborn E, Stål O, Jansson A. Estrogen and androgen-converting enzymes 17β-hydroxysteroid dehydrogenase and their involvement in cancer: with a special focus on 17β-hydroxysteroid dehydrogenase type 1, 2, and breast cancer. Oncotarget. 2017;8:30552–62.
Young G, Kusakabe M, Nakamura I, Lokman PK, Goetz FW. Gonadal steroidogenesis in teleost fish. In: Melamed P, Sherwood N, editors. Hormones and their receptors in fish reproduction. USA: World Scientific; 2005. p. 155–223.
So WK, Kwok HF, Ge W. Zebrafish gonadotropins and their receptors: II. Cloning and characterization of zebrafish follicle-stimulating hormone and luteinizing hormone subunits—their spatial-temporal expression patterns and receptor specificity. Biol Reprod. 2005;72:1382–96.
Lien S, Koop BF, Sandve SR, Miller JR, Kent MP, Nome T, Hvidsten TR, Leong JS, Minkley DR, Zimin A, et al. The Atlantic salmon genome provides insights into rediploidization. Nature. 2016;533:200–5.
Matsubara T, Honda S, Wada T, Soyano K, Hara A. Seasonal changes in serum levels of vitellogenin and estradiol-17β related to sexual maturation in rearing female Japanese sardine Sardinops melanostictus. Bull Hokkaido Natl Fish Res Inst. 1995;59:19–29.
Matsuyama M, Adachi S, Nagahama Y, Kitajima C, Matsuura S. Annual reproductive cycle of the captive female Japanese sardine Sardinops melanostictus: relationship to ovarian development and serum levels of gonadal steroid hormones. Mar Biol. 1991;108:21–9.
Montserrat N, González A, Méndez E, Piferrer F, Planas JV. Effects of follicle stimulating hormone on estradiol-17β production and P-450 aromatase (CYP19) activity and mRNA expression in brown trout vitellogenic ovarian follicles in vitro. Gen Comp Endocrinol. 2004;137:123–31.
Tyler CR, Sumpter JP, Kawauchi H, Swanson P. Involvement of gonadotropin in the uptake of vitellogenin into vitellogenic oocytes of the rainbow trout, Oncorhynchus mykiss. Gen Comp Endocrinol. 1991;84:291–9.
Nagahama Y, Yamashita M. Regulation of oocyte maturation in fish. Develop Growth Differ. 2008;50:S195–219.
Nyuji M, Kitano H, Shimizu A, Lee JM, Kusakabe T, Yamaguchi A, Matsuyama M. Characterization, localization, and stage-dependent gene expression of gonadotropin receptors in chub mackerel (Scomber japonicus) ovarian follicles. Biol Reprod. 2013;88:148.
Guzmán JM, Luckenbach JA, Yamamoto Y, Swanson P. Expression profiles of Fsh-regulated ovarian genes during oogenesis in coho salmon. PLoS One. 2014;9:e114176.
Nakamura I, Evans JC, Kusakabe M, Nagahama Y, Young G. Changes in steroidogenic enzyme and steroidogenic acute regulatory protein messenger RNAs in ovarian follicles during ovarian development of rainbow trout (Oncorhynchus mykiss). Gen Comp Endocrinol. 2005;144:224–31.
Luckenbach JA, Dickey JT, Swanson P. Follicle-stimulating hormone regulation of ovarian transcripts for steroidogenesis-related proteins and cell survival, growth and differentiation factors in vitro during early secondary oocyte growth in coho salmon. Gen Comp Endocrinol. 2011;171:52–63.
Luckenbach JA, Iliev DB, Goetz FW, Swanson P. Identification of differentially expressed ovarian genes during primary and early secondary oocyte growth in coho salmon, Oncorhynchus kisutch. Reprod Biol Endocrinol. 2008;6:2.
Kobayashi M, Tanaka M, Fukada S, Nagahama Y. Steroidogenesis in the ovarian follicles of the medaka (Oryzias latipes) during vitellogenesis and oocyte maturation. Zool Sci. 1996;13:921–7.
Matsuyama M, Shiraishi T, Sundaray JK, Rahman M, Ohta K, Yamaguchi A. Steroidogenesis in ovarian follicles of chub mackerel, Scomber japonicus. Zool Sci. 2005;22:101–10.
Ohta K, Mine T, Yamagichi A, Matsuyama M. Steroidogenic pathway to estradiol-17β synthesis in the ovarian follicles of the protogynous wrasse, Pseudolabrus sieboldi. Zool Sci. 2001;18:937–45.
Ohta K, Yamaguchi S, Yamaguchi A, Okuzawa K, Gen K, Kagawa H, Matsuyama M. Biosynthesis of estradiol-17β in the ovarian follicles of the red seabream Pagrus major during vitellogenesis. Fish Sci. 2002;68:680–7.
Taranger GL, Carrillo M, Schulz RW, Fontaine P, Zanuy S, Felip A, Weltzien FA, Dufour S, Karlsen Ø, Norberg B, Andersson E, Hansen T. Control of puberty in farmed fish. Gen Comp Endocrinol. 2010;165:483–515.
Lambert Y, Thorsen A. Integration of captive and wild studies to estimate egg and larval production of fish stocks. J Northwest Atl Fish Sci. 2003;33:71–9.
Nyuji M, Kazeto Y, Izumida D, Tani K, Suzuki H, Hamada K, Mekuchi M, Gen K, Soyano K, Okuzawa K. Greater amberjack Fsh, Lh, and their receptors: plasma and mRNA profiles during ovarian development. Gen Comp Endocrinol. 2016;225:224–34.
Hongo Y, Yabuki A, Fujikura K, Nagai S. Genes functioned in kleptoplastids of Dinophysis are derived from haptophytes rather than from cryptophytes. Sci Rep. 2019;9:9009.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc. Accessed 22 Feb 2018.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12:59–60.
R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2016. https://www.r-project.org/. Accessed 22 Feb 2018.
Gene Ontology Consortium. The gene ontology (GO) database and informatics resource. Nucl Acids Res. 2004;32:D258–61.
Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357–9.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57:289–300.
Nyuji M, Hamada K, Kazeto Y, Mekuchi M, Gen K, Soyano K, Okuzawa K. Photoperiodic regulation of plasma gonadotropin levels in previtellogenic greater amberjack (Seriola dumerili). Gen Comp Endocrinol. 2018;269:149–55.
We thank Dr. Yuichi Ozaki (Fisheries Technology Institute, Japan Fisheries Research and Education Agency) for assistance with time-resolved fluoroimmunoassay. We thank Edanz (www.edanzediting.co.jp) for editing the English text of a draft of this manuscript.
The present work was financially supported in part by Grant-in-Aid for Young Scientists (B) (No. 15 K18738) and Challenging Research (Pioneering) (No. 19H05540) from the Japan Society for the Promotion of Science.
Ethics approval and consent to participate
All experimental procedures followed the guidelines for animal welfare of Fisheries Research and Education Agency, Japan (50322001) and were approved by the committee of animal welfare of National Research Institute of Fisheries and Environment of Inland Sea (Permit Number: 2016-2).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Results of homology search against Atlantic herring (Clupea harengus), zebrafish (Danio rerio), and medaka (Oryzias latipes).
GO classification of assembled contigs.
Complete list of enriched GO terms in the gonads. Table S1. Enriched GO terms in the ovary. Table S2. Enriched GO terms in the testis.
Sequence alignments of the amino acid sequences of Gth subunits (Figure S2), Gthrs (Figure S3 and S4), and steroidogenic proteins (Figure S5–S9).
Molecular characterization of Japanese sardine Fshb. Figure S10. Sequence alignment of the amino acid sequences of Japanese sardine Fshb isoforms.
Primers used for cDNA cloning and real-time PCR. Table S3. Primers used for cloning cDNAs. Table S4. Primers used for amplifying PCR products to construct a standard curve. Table S5. Primers used for real-time PCR.
GenBank accession numbers of contigs used for the clustering analysis.
About this article
Cite this article
Nyuji, M., Hongo, Y., Yoneda, M. et al. Transcriptome characterization of BPG axis and expression profiles of ovarian steroidogenesis-related genes in the Japanese sardine. BMC Genomics 21, 668 (2020). https://doi.org/10.1186/s12864-020-07080-1
- Clupeoid fish
- Steroidogenic enzyme