Abstract Background The Japanese eel (Anguilla japonica) holds significant economic value in East Asia, but limitations in understanding its reproductive biology have hindered advancements in artificial breeding techniques. Previous research has primarily focused on conserved sex differentiation genes, offering limited insights into the broader molecular mechanisms driving gonadal development and sexual dimorphism. To address these limitations, this study aims to investigate key genes and pathways involved in gonadal development through a comprehensive transcriptomic analysis of male and female eel gonads. Results PacBio Iso-Seq and Illumina RNA-Seq technologies were combined to conduct a full-length transcriptome analysis of male and female Japanese eel gonads at a post-differentiation, pre-maturation stage. A total of 24,661 unigenes were identified in ovaries and 15,023 in testes, along with genomic regulatory elements such as transcription factors, simple sequence repeats, and long non-coding RNAs. Additionally, 1,210 differentially expressed genes were detected. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses revealed significant pathways involved in cell cycle regulation, metabolic processes, apoptosis, and hormone activity. Notably, several reproductive-related genes, including bambi, ccnb1, cdc20, gdf9, prlh, ccdc39, chrebp, tspo, syce3, and ngb, demonstrated significant dimorphic expression in eel gonads. Conclusions This study provides valuable insights into the molecular mechanisms of gonadal differentiation and sexual dimorphism in Japanese eels. The findings expand the genetic resources available for the eel breeding industry and could facilitate the development of improved artificial breeding techniques focused on reproductive development. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-025-11279-5. Keywords: Japanese eel (Anguilla japonica), Transcriptome analysis, Testis, Ovary, Sexual dimorphism Background The Japanese eel (Anguilla japonica) is a catadromous species widely distributed across the Western Pacific, particularly in Vietnam, the Philippines, China, the Korean Peninsula, and Japan [[42]1]. Valued for its flavor, rich nutritional profile, and economic importance, it is a crucial aquaculture species in East Asia [[43]2]. However, limited knowledge of its life cycle and reproductive biology has restricted the success of artificial breeding, with larvae still mainly sourced from wild populations [[44]3]. Overfishing and habitat loss have further diminished wild stocks, threatening the sustainability of the eel industry [[45]4–[46]7]. Therefore, advancing the understanding of the reproductive biology of the Japanese eel is essential for both resource conservation and the long-term viability of the industry. The eel’s life cycle encompasses both marine and freshwater stages. Gonadal differentiation does not occur during the marine leptocephalus (larvae) and glass eel (juvenile) stages, but begins in the freshwater yellow eel stage (first metamorphosis), with differentiation influenced by factors such as growth rate and population density [[47]8]. Individuals exhibiting faster growth rates generally develop as males, while slower-growing individuals tend to become females. These factors play a pivotal role in sex determination, ultimately preparing the eel for reproductive migration during the silver eel stage (second metamorphosis), at which point gonads are fully differentiated, though still immature [[48]8, [49]9]. Sexual dimorphism in eels is intricately tied to their respective survival strategies. Males adopt a "time-minimizing strategy," reaching sexual maturity at smaller body sizes to maximize reproductive opportunities. For example, studies on Japanese eels during early sex differentiation (body weight < 200 g) showed that males had significantly higher growth rates, with daily growth rates of 2.06% vs. 1.55% in females [[50]10]. Annual growth rates were 124 mm/year for males and 116 mm/year for females until the yellow eel stage [[51]11]. Similar trends have been observed in European and American eels, where males exhibit faster early growth to accelerate reproductive readiness [[52]12]. In contrast, females follow a "size-maximizing strategy," prioritizing prolonged growth to achieve significantly larger asymptotic sizes [[53]8]. As growth progresses, females surpass males, ultimately reaching final body lengths (L∞) of 596–794 mm, significantly larger than the 428–557 mm seen in males [[54]9–[55]13]. This strategy enhances both fecundity and survival during migration. Environmental factors, particularly population density, further influence these growth patterns and sex ratios, with higher densities favoring male differentiation—likely due to elevated cortisol levels that suppress estrogen production—while lower densities promote female development and larger body sizes [[56]9]. These interactions underscore the complexity of sexual dimorphism in eels. At the molecular level, sexual dimorphism is regulated by sex-biased gene expression. Key genes such as dmrt1, gsdf, and amh are instrumental in driving testicular differentiation in males, while cyp19a1a and zp3 are involved in ovarian development in females[[57]14–[58]17]. These conserved genes play essential roles in sex differentiation across species, yet a comprehensive understanding of the broader molecular mechanisms underpinning reproductive development and sexual dimorphism in eels remains elusive. This study combines Illumina RNA-Seq and PacBio Iso-Seq technologies to provide a comprehensive transcriptomic analysis of gonadal tissues from male and female Japanese eels at a post-differentiation, early developmental stage, before sexual maturation. By comparing sex-specific gene expression profiles, this research aims to identify key genes and signaling pathways involved in sexual dimorphism and gonadal development, offering new insights into the molecular regulation of sex differentiation in Japanese eels. Method Sample collection Japanese eels were collected from a commercial eel farm in Guangdong Province, China. The eels selected for this study were farmed individuals that had not undergone artificial induction for maturation. The samples included three female eels (O1, O2, O3; 997.3 ± 143.1 g; 841.7 ± 18 mm) and three male eels (T1, T2, T3; 407.0 ± 18.5 g; 691.3 ± 2.1 mm). To minimize suffering, the eels were anesthetized using MS-222 at a concentration of 0.1% prior to sampling. Gonadal tissues (ovaries and testes) were dissected, with portions of each sample fixed in Bouin's solution for histological examination. The remaining tissues were immediately frozen in liquid nitrogen and stored at −80 °C until total RNA extraction. Tissue sectioning Gonadal tissues fixed in Bouin’s solution were kept at room temperature overnight. After fixation, the samples were embedded in paraffin and sectioned at a thickness of 4 μm. The tissue sections were deparaffinized, rehydrated, and stained with hematoxylin–eosin (HE) for histological examination. RNA extraction and quality assessment Total RNA was extracted from the six gonadal tissue samples using RNAiso Plus (Takara, Japan). The RNA quality and integrity were initially evaluated through 1% agarose gel electrophoresis. The purity of the RNA (OD260/OD280 ≈ 1.8–2.0) was measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA), and RNA integrity was further assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Only RNA samples that met the quality criteria were used for PacBio and Illumina library construction and sequencing, which were performed by Novogene Co., Ltd. (Beijing, China). Full-length transcriptome sequencing, library construction, and data analysis For PacBio sequencing, total RNA from the ovaries of three female Japanese eels was pooled into one sample, and total RNA from the testes of three male eels was pooled into another, resulting in one female and one male RNA pool. Full-length transcriptome sequencing was performed using the PacBio Sequel II platform. To construct the PacBio sequencing library, mRNA with polyA tails was enriched from the pooled RNA samples using Oligo(dT) magnetic beads and reverse-transcribed into cDNA using the SMARTer PCR cDNA Synthesis Kit. The cDNA was amplified via PCR, optimized for cycle conditions, and subjected to damage repair, end repair, and ligation with SMRTbell adapters to create the full-length transcriptome library. Exonuclease digestion was performed to remove unligated adapters, and the final SMRTbell library was completed for sequencing. The raw PacBio sequencing data were processed using the SMRT Link software V8.0.0 to correct subreads. Circular consensus sequences (CCS) were generated, and full-length non-chimeric (FLNC) reads were identified based on the presence of 5’ and 3’ primers and polyA tails. FLNC sequences were clustered using the n*log(n) clustering method provided by the official PacBio Iso-Seq software to generate consensus sequences, which were polished using the Arrow software (Pacific Biosciences) to produce high-quality sequences. LoRDEC V0.7 software was employed to correct errors in the PacBio long-read sequences using the high-accuracy short-read Illumina RNA-seq data[[59]18]. Finally, redundancy was removed from the high-quality transcripts using CD-HIT V4.6.8, resulting in the final full-length transcriptome reference for the Japanese eel [[60]19]. Illumina transcriptome library construction, sequencing, and data analysis For Illumina RNA-seq, total RNA was extracted from each of the six individual samples (three female ovaries, three male testes), representing three biological replicates per sex. mRNA with polyA tails was enriched from each sample using Oligo (dT) magnetic beads and fragmented into short segments. These fragmented mRNA segments were used as templates for cDNA synthesis. The double-stranded cDNA was purified, end-repaired, A-tailed, and ligated with sequencing adapters. Fragments of 250–300 bp were selected and amplified via PCR to construct the sequencing library. The libraries were sequenced on the Illumina NovaSeq 6000 platform. Raw sequencing data were processed to obtain clean reads by removing adapter contamination, low-quality reads, and reads containing undetermined bases. Clean reads were assessed for Q20, Q30, and GC content to ensure sequencing quality. Gene annotation The deduplicated sequences obtained after CD-HIT processing were aligned to the following nucleotide and protein databases for annotation: NR (Non-Redundant Protein Database), NT (NCBI Non-Redundant Nucleotide Sequences), Pfam (Protein Families Database), KOG (Cluster of Orthologous Groups of Proteins), Swiss-Prot (Swiss-Prot Protein Sequence Database), KEGG (Kyoto Encyclopedia of Genes and Genomes), and GO (Gene Ontology). Prediction of Transcription Factors (TFs), Simple Sequence Repeats (SSRs), and Long Non-Coding RNAs (LncRNAs) Transcription factors (TFs) were predicted by aligning the sequences to the animal TFDB V2.0 database. Simple sequence repeats (SSRs) were identified using MISA V1.0 software[[61]20, [62]21]. Long non-coding RNAs (LncRNAs) were predicted using a combination of CNCI (Coding-Non-Coding Index V2.0) [[63]22], PLEK (Predictor of Long Non-Coding RNAs and Messenger RNAs Based on an Improved K-mer Scheme V1.2), CPC (Coding Potential Calculator V2.0) [[64]23], and Pfam Scan (V1.6) [[65]24]. The transcripts with lengths greater than 200 nt but without coding function were defined as candidate lncRNAs. Gene expression quantification The PacBio Iso-Seq data from both male and female samples were merged to construct a unified reference transcriptome. Clean reads from each of the six individual Illumina RNA-seq samples (three male, three female) were first aligned to the reference genome using Bowtie2. The aligned reads were then processed with RSEM software to obtain read counts [[66]25]. To account for differences in sequencing depth and transcript length, all read counts were converted to FPKM (Fragments Per Kilobase of transcript per Million mapped reads) for gene expression analysis. Principal Component Analysis (PCA) To visualize the overall variation in gene expression between male and female samples, PCA was performed using FPKM values as input. PCA was conducted using the FactorMineR package in R (v3.1.2), and the major sources of variation between the samples were plotted. Differential expression analysis and functional annotation For differential expression analysis, read counts were used as input for DESeq2 (v1.10.1) [[67]26]. Normalization and statistical analysis were performed to identify differentially expressed genes between male and female samples. Genes were considered differentially expressed if |log2(fold change)|> 1 and adjusted p-value (P[adj]) < 0.05. Functional annotation was performed using GO terms, and pathway analysis was conducted using the KEGG database to interpret the biological significance of the differentially expressed genes. Real-Time Quantitative PCR (RT-qPCR) validation To validate the transcriptome sequencing results, differentially expressed genes were selected for RT-qPCR analysis. Quantitative primers (Table [68]1) were designed based on sequence information from NCBI. Total RNA from male and female Japanese eel gonadal tissues was extracted using RNAiso Plus (TaKaRa, Japan) and reverse-transcribed using the HiScript III 1st Strand cDNA Synthesis Kit (+ gDNA wiper) (Vazyme, China). RT-qPCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, China), and amplification products were analyzed through melting curve analysis. Relative gene expression levels were calculated using the 2^−ΔΔCt method with β-actin as the internal control[[69]27]. Statistical analysis was performed using t-tests, with p-values less than 0.05 considered significant. Bar graph generation was performed using GraphPad Prism 9 and R. In the transcriptome data, some genes showed zero expression, leading to undefined values during ratio calculation. To resolve this, a small constant (0.1) was added during log transformation, ensuring proper visualization and preventing errors. Table 1. Primer sequences Gene name Primer Primer sequence (5’ → 3’) Fragment size (bp) prlh F CTACAGGGCAGGGAACTGAAA 134 R AGTGACAGCACACCACTTCA ngb F GTTACAATGTCATATACTGTTGCCA 102 R TACTGTGAGAAGGCCAACGG syce3 F GTTGGCACTGGAAATTGCGG 163 R ACCCATTTGCTCAAACACAACA tspo F CTCGTTCCTTTGGGCCTGTA 197 R AGCCACATCAGGTAAGGGGT chrebp F CGGAATGGTGTGTCGGGATT 133 R GCTCCAGCAGTGTAAAATGGC ccdc39 F CCTGTGTTCCCGATAGCCATT 155 R GCGCTACCAGCAGAGCTTA gdf9 F CCAAGGGCCATGGGATTCAT 163 R TGTAGGCAATGGAGCCATCC bambi F TGCACACTCCTCTGTATGTTCA 127 R CTGAGCAGAAGTGCACCAGTT ccnb1 F GGACTGCAATTGGGAGCAAAC 140 R AGATTGGAACCCCCTGCAC cdc20 F CTGCCTGTCTGGGCTACTTTG 145 R ATAGTCCACACCTCTCCAGCC amh F TCTCCTAGCAACGGATTGGC 80 R CACAGGAAGTGGAAAGTCCGG gsdf F GGCAAGGTGAACAGATGCCA 151 R GGGCAGTGTCTCCTAAGACC figla F GAGCGCCTACGAATCAGGAA 203 R GCCATGCCCTCTGGTATCTC sox3 F GCGATGGGCTTGGGTTCTAT 239 R GTTAGGGGTAGCGTTCCGTT transcript26973/f2p0/1058 F ATGGTCTTGCATGTGTCGATG 191 R TTAGCGGCACACGCTTTTGA transcript28036/f3p0/754 F CGTGCGTGCGTTAATTGTTAGAA 132 R ATCCATCTTCCATTGTCTGGC transcript18948/f2p0/1878 F TCTTGCTAAGGGGGCTGTCCT 141 R ATACCTCAACAGCATTGACTCGT transcript24482/f2p0/1416 F CAGCACACACTGAAAACCACTCC 100 R CTCTTTGCAGAAACCTTGCCCC transcript25535/f2p0/1266 F TGATACCAAGTGATTAAATGCGGC 138 R CTGGATTGGGGCTCTCAAAG transcript6321/f2p0/3250 F CTTGGCTCCTCCTCCCCT 151 R AAGGCTCACACTTCCTGTATTGTG transcript46862/f2p0/2206 F ATTTGGTTGGCAGTGAAGTTGGAT 210 R TGGTGGAAGCCAGCTCTCC transcript50895/f5p0/1514 F AGGCAGTTGGGTCCTACTCAT 138 R CTCTTTGTCCAATCCATTGCCTTGA transcript52751/f2p0/1264 F TCTTGAACACTCCGATTCCGTC 158 R GCAGATTTTCAGCTTCCTCCAAA transcript53752/f7p0/1016 F GAGGATTCTAAAGCCCGGGAG 101 R TCCCAAAATAAACACAAAACAGCG transcript54430/f2p0/569 F GCAGCAGGTGGAAGTACCATA 179 R CCTGTTGGGGCCTGGTATTG transcript54502/f12p0/429 F CCATTGAGAGATAGCCTGGGAC 125 R TTGCACACTGATGAAACTGCCA β-actin F TCACCACCACAGCCGAAAGG 219 R CGCAGGATTCCATTCCCAGGA [70]Open in a new tab Results Sex identification of Japanese Eel To ensure accurate sex identification of the Japanese eel samples, we employed both histological examination and the analysis of sex differentiation gene expression. Gonadal development in Japanese eels involves spermatogenesis in males (spermatocytes, spermatids, and spermatozoa) and oogenesis in females (primary oocytes, meiotic oocytes, and previtellogenic oocytes). Ovarian development is divided into six stages based on oocyte diameter: pre-vitellogenic (PV) (< 200 μm), early vitellogenic (EV) (200–400 μm), mid-vitellogenic (MV) (400–600 μm), late-vitellogenic (LV) (600–800 μm), and migratory nucleus (MN) (800–1000 μm), reflecting the progression of oocyte maturation [[71]28]. As shown in Fig. [72]1, all selected eels were in the immature stage of gonadal development. In female eels, the ovaries contained oocytes with diameters of less than 200 μm (Fig. [73]1A-B), corresponding to the PV (pre-vitellogenic) stage. At this stage, the ovaries consisted primarily of oogonia and early-stage oocytes, with no visible yolk granules. The cytoplasm of the oocytes was filled with lipid vacuoles, and the nuclei contained clear nucleoli, indicative of early ovarian development [[74]29]. In male eels, the testes contained only spermatogonia, confirming their immature stage (Fig. [75]1C-D). Spermatogonia were the predominant cell type, with no progression to later stages of spermatogenesis. Fig. 1. [76]Fig. 1 [77]Open in a new tab Morphological Observation of Gonadal Tissue in Japanese Eel. A Anatomical diagram of female Japanese eel, with the dashed border indicating the ovarian tissue. B HE staining result of ovarian tissue section. Red arrow indicates the nucleolus (N), blue arrow indicates oil droplets (O). C Anatomical diagram of male Japanese eel, with the dashed border indicating the testicular tissue. D HE staining result of testicular tissue section. sg, spermatogonia. Scale bar: 50 μm. E Relative mRNA levels of sex differentiation genes in male and female gonads. The asterisks indicate statistical significance, where * denotes P < 0.05; *** denotes P < 0.001; and **** denotes P < 0.0001 Additionally, the expression of conserved sex differentiation genes, including sox3, figla, amh, and gsdf, was analyzed to further confirm the sex of the samples. As shown in Fig. [78]1E, sox3 and figla were highly expressed in females, while amh and gsdf were predominantly expressed in males. This gene expression profile further validated the histological findings, confirming the sex of the eels at the molecular level. Transcriptome data of male and female Japanese Eel Gonads The PacBio sequencing produced 20,802,769 subreads from the ovarian sample and 10,587,337 subreads from the testicular sample. After correction, the ovarian sample generated 840,775 CCS with an average length of 2,828 bp, while the testicular sample produced 351,284 CCS with an average length of 2,589 bp. FLNC sequences were identified, yielding 742,832 FLNC sequences from the ovaries and 310,823 from the testes. After correction using Illumina data and deduplication via CD-HIT, 24,661 unigenes were identified in the ovaries and 15,023 unigenes in the testes (Table [79]2). Table 2. Statistics on the full-length transcriptome sequence types of Japanese eel gonads Type Total number Ovary Testis Subreads 20,802,769 10,587,337 Subreads base (G) 53.46 23.14 CCS 840,775 351,284 Average read length of CCS 2828 2589 FLNC 742,832 310,823 Consensus isoforms 54,534 28,373 Unigene 24,661 15,023 [80]Open in a new tab CCS Circular consensus sequence FLNC Full-length non chimera Illumina sequencing generated a total of 263,909,230 raw reads across the six transcriptome libraries (three biological replicates per sex). After quality filtering, a total of 254,789,618 clean reads were retained, with clean bases exceeding 5.92 Gb per sample. The sequencing error rates for all libraries were ≤ 0.03%, with Q20 values ≥ 97.06% and Q30 values ≥ 92.77%, indicating high sequencing quality. Alignment of the Illumina clean reads to the PacBio-derived full-length transcriptome unigenes for differential gene analysis showed alignment rates ranging from 59.73% to 81.89% across the six samples (Table [81]3). Table 3. Quality monitoring of the transcriptome sequencing data of Japanese eel gonads using Illumina Sample Raw data Clean data Clean bases (Gb) Error rate (%) Q20 (%) Q30 (%) Mapped ratio (%) O1 47,064,916 45,247,062 6.79G 0.03 97.63 93.69 77.64% O2 41,101,038 39,493,324 5.92G 0.03 97.78 94.02 77.91% O3 47,655,174 45,765,970 6.86G 0.03 97.56 93.61 81.89% T1 40,896,232 39,703,236 5.96G 0.03 97.43 93.52 69.45% T2 45,463,992 44,200,344 6.63G 0.03 97.43 93.46 60.97% T3 41,727,878 40,379,682 6.06G 0.03 97.06 92.77 59.73% [82]Open in a new tab O Ovary T Testis Transcript annotation Functional annotation of the deduplicated unigenes was performed by aligning them to seven major public databases. The annotation results showed that 24,209 unigenes (98.17%) from the ovaries and 14,205 unigenes (94.56%) from the testes were successfully annotated using the public databases (Table [83]4). Table 4. Database annotation results Annotation database Ovary Testis GO 18,133 (73.53%) 10,308 (68.61%) KEGG 23,879 (96.83%) 13,645 (90.83%) KOG 18,596 (75.41%) 10,644 (70.85%) Pfam 18,133 (73.53%) 10,308 (68.61%) Swissprot 23,253 (94.29%) 13,225 (88.03%) NT 21,150 (85.76%) 12,080 (80.41%) NR 23,988 (97.27%) 13,767 (91.64%) All annotated 24,209 (98.17%) 14,205 (94.56%) [84]Open in a new tab GO Gene Ontology KEGG Kyoto Encyclopedia of Genes and Genomes KOG Cluster of Orthologous Groups of proteins Pfam Protein families database Swissprot SwissProt protein sequence database NT NCBI non-redundant nucleotide sequences NR Non-Redundant Protein Database Identification of TFs, SSRs, and LncRNAs TFs are proteins that bind to specific DNA sequences to regulate gene expression. We identified 1,499 TFs from 60 families in the ovaries (Supplementary data 1) and 830 TFs from 44 families in the testes (Supplementary data 2). The top 30 TF families are shown in Fig. [85]2A-B. (details in Supplementary data 1& 2). The most abundant TF families in the ovaries were zf-C2H2 (Zinc Finger-C2H2) with 470 members, ZBTB (Zinc Finger and BTB Domain Containing) with 150 members, and BHLH (Basic Helix-Loop-Helix) with 79 members. In the testes, the dominant families were zf-C2H2, ZBTB, and HMG (High Mobility Group) with 329, 80, and 52 members, respectively. Fig. 2. [86]Fig. 2 [87]Open in a new tab Prediction of Transcription Factors (TF), and Simple Sequence Repeats (SSR) in the Ovaries and Testes of Male and Female Japanese Eel. A and B Distribution of transcription factor (TF) families in the ovaries and testes of Japanese eel; C and D SSRs in the ovaries and testes of Japanese eel SSRs are short repetitive sequences composed of 1–6 nucleotide units that are uniformly distributed throughout the eukaryotic genome. Using MISA software, we detected SSRs in the unigenes of Japanese eel ovaries and testes (Table [88]5). A total of 26,456 potential SSRs were identified from 143,533 unigenes in the ovaries (Supplementary data 3), while 13,317 potential SSRs were identified from 7,907 unigenes in the testes (Supplementary data 4). Mononucleotide, dinucleotide, and trinucleotide repeats were the predominant types, accounting for 79.53% of all SSRs in the ovaries and 81.87% in the testes. In both tissues, the most common dinucleotide repeat motif was CA/TG, representing 38.05% in the ovaries and 34.93% in the testes. For trinucleotide repeats, the predominant motif was GAG/CTC, accounting for 19.21% in the ovaries and 19.01% in the testes. Another notable difference was observed in CAG/CTG, which was more abundant in the ovaries (9.1%) compared to the testes (6.08%). Additionally, GCA/TGC repeats had a slightly higher proportion in the ovaries (4.55%) compared to the testes (3.57%), while ATT/AAT repeats were more frequent in the testes (5.25%) than in the ovaries (2.53%). Among tetranucleotide repeats, TTTG/CAAA was the most abundant type in both tissues, representing 7.05% in the ovaries and 7.69% in the testes. Other high-frequency repeat motifs also displayed distinct differences between the ovaries and testes (Fig. [89]2C-D), allowing for a clearer comparative understanding of SSR distribution across the two tissues. Table 5. Prediction of SSR loci in Japanese eel testes and ovaries Type Ovary Testis number percentage (%) number percentage (%) Mono- 14,600 55.19% 7733 58.07% Di- 3669 13.87% 1855 13.93% Tri- 2770 10.47% 1315 9.87% Tetra- 312 1.18% 169 1.27% Penta- 67 0.25% 48 0.36% Hexa- 53 0.20% 23 0.17% C 4985 18.84% 2174 16.32% [90]Open in a new tab Mono- Mono-nucleotide repeat Di- Di-nucleotide repeat Tri- Trinucleotide repeat Tetra- Tetra nucleotide repeat Penta- Penta-nucleotide repeat Hexa- Hexa nucleotide repeat C Complex repeat type LncRNAs are RNA molecules longer than 200 nucleotides that do not encode proteins. In total, 347 high-confidence lncRNAs were identified in the ovaries and 694 in the testes (Fig. [91]3A-B, see Supplementary data 5 and 6). Comparative analysis revealed 54 lncRNAs with significant differences in expression between the two tissues (|log2(fold change)|> 1, P[adj] < 0.05). Of these, 22 lncRNAs showed higher expression in the ovaries, while 32 were more highly expressed in the testes (Fig. [92]3C, Supplementary data 7). RT-qPCR validation of 12 selected lncRNAs confirmed the consistency of the transcriptome sequencing results (Fig. [93]3D). Fig. 3. [94]Fig. 3 [95]Open in a new tab Analysis of Long-noncoding RNAs (LncRNA) in the Ovaries and Testes of Japanese Eel. A and B Venn diagrams of LncRNA identified in the ovaries and testes of Japanese eel using four different screening methods. C Volcano plot showing differentially expressed LncRNAs between ovaries and testes. Red dots represent genes with higher expression in ovaries, green dots represent genes with higher expression in testes, and blue dots indicate no significant difference. D Validation of LncRNAs using RT-qPCR PCA analysis PCA was performed to assess the overall variation in gene expression between male and female gonadal tissues. The PCA plot (Fig. [96]4A) showed a clear separation between ovary samples and testis sample. Although there was some variation within the male samples, as shown by their slightly greater dispersion, the separation between the sexes remained distinct. Fig. 4. [97]Fig. 4 [98]Open in a new tab Enrichment Analysis and Validation of Differential Expression Genes (DEGs) in Male and Female Japanese Eel. A Two-dimensional clustering map of principal component analysis (PCA) for six samples of male and female Japanese eel; B Volcano plot illustrating DEGs between ovaries and testes. Blue dots represent genes with higher expression in testis, green dots represent genes with higher expression in ovaries, and gray dots indicate genes with no statistically significant difference in expression between the two tissues; C GO analysis of DEGs in Japanese eel gonads; D KEGG enrichment of DEGs in Japanese eel gonads; E Validation of RNA-Seq data using RT-qPCR Identification and enrichment analysis of Differentially Expressed Genes (DEGs) A total of 1,210 DEGs were identified (|log2(fold change)|> 1, P[adj] < 0.05), with 688 showing higher expression in the ovaries and 522 exhibiting higher expression in the testes (Fig. [99]4B, see Supplementary data 8 and 9). GO analysis categorized DEGs into three main groups: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) Fig. [100]4C shows the top 15 pathways significantly enriched across the three major categories of the GO database. In the BP category, DEGs were enriched in terms covering fundamental biological processes. In the CC category, the most enriched term was "membrane," suggesting a strong association with membrane components. The only enriched MF term was "glutathione hydrolase activity," indicating involvement in oxidative stress response [[101]30]. KEGG pathway enrichment analysis identified 258 pathways, with the top 20 most significantly enriched pathways shown in Fig. [102]4D. These pathways are related to cell cycle regulation, apoptosis, phagocytosis, and cell adhesion, all essential for reproductive tissue development and maintenance. GO and KEGG pathway enrichment analysis also revealed several DEGs enriched in GO terms related to hormone activity, synaptonemal complex assembly, as well as KEGG pathways such as neuroactive ligand-receptor interaction, reproduction, TGF-beta signaling pathway, and oocyte meiosis. The genes enriched in these terms nd pathways are listed in Supplementary Table 10. Validation of RNA-seq Results by RT-qPCR To verify the DEGs between testes and ovaries, we selected 10 genes that are involved in reproductive-related pathways for RT-qPCR validation (Table [103]6). Table 6. Identification of differentially expressed gene candidates in the gonads of Japanese eel Name Function annotation Ovary-FPKM Testis-FPKM log2FoldChange P[adj] prlh prolactin releasing hormone 61.15 ± 10.49 0 Inf 1.66 × 10^–5 ngb neuroglobin-like 0 5.58 ± 5.41 -Inf 4.83 × 10^–2 syce3 synaptonemal complex central element protein 3 0.08 ± 0.13 20.21 ± 9.52 −7.52 1.39 × 10^–2 tspo translocator protein 0.35 ± 0.17 75.81 ± 14.50 −7.33 1.29 × 10^–3 chrebp carbohydrate-responsive element-binding protein 0.27 ± 0.07 69.65 ± 30.99 −7.65 8.85 × 10^–4 ccdc39 coiled-coil domain-containing protein 39-like 0 9.31 ± 1.70 -Inf 1.15 × 10^–3 gdf9 growth differentiation factor 9 37.34 ± 12.50 1.16 ± 0.79 5.37 1.47 × 10^–2 bambi BMP and activin membrane-bound inhibitor 36.51 ± 19.18 0.85 ± 0.73 5.83 6.38 × 10^–3 ccnb1 cyclin B1 106.20 ± 121.22 0.14 ± 0.25 9.94 4.69 × 10^–2 cdc20 cell division cycle 20 homolog 265.99 ± 223.58 5.41 ± 2.89 5.99 2.60 × 10–2 [104]Open in a new tab Inf Infinity Among the genes with significantly higher expression in ovaries compared to testes, we chose bambi (BMP and activin membrane-bound inhibitor) and gdf9 (growth differentiation factor 9), which were enriched in the KEGG pathways for Ovarian steroidogenesis and TGF-beta signaling pathway, respectively. Additionally, we included cdc20 (cell division cycle 20), enriched in the Oocyte meiosis pathway, and ccnb1 (cyclin B1), enriched in the FoxO signaling pathway, as well as prlh (prolactin-releasing hormone), which are involved in hormone synthesis and secretion regulation. For the genes showing higher expression in testes, we selected ccdc39 (coiled-coil domain-containing protein 39-like) and syce3 (synaptonemal complex central element protein 3), both of which are involved in meiosis and gametogenesis, as well as tspo (translocator protein) and chrebp (carbohydrate-responsive element-binding protein), related to steroid hormone synthesis. Additionally, we included ngb (neuroglobin-like), known for its role in developmental and reproductive regulation. The RT-qPCR results align closely with the transcriptome data (Fig. [105]4E), confirming the reliability of the transcriptome analysis for these genes. Discussion The artificial breeding of eels has long been hindered by a limited understanding of reproductive regulation. Despite extensive efforts, key aspects of gonadal differentiation and development remain unclear, posing challenges to breakthroughs in aquaculture. In recent years, several research groups have employed transcriptome analyses to explore reproductive development. For example, Burgerhoute et al. [[106]31]analyzed ovarian tissues during the artificial maturation of European eels, identifying key genes involved in steroidogenesis. Lai et al.[[107]28]examined dynamic transcriptome changes in ovarian follicles of Japanese eels, revealing genes that regulate vitellogenesis and oocyte maturation across various stages of follicle development. Babio et al.[[108]32] explored liver tissues in Anguilla australis, focusing on how cell junctions regulate vitellogenin uptake during vitellogenesis, while Hsu et al. [[109]33]studied sexually dimorphic gene expression in the pectoral fins of silver eels, uncovering regulatory differences between the sexes tied to secondary sexual characteristics. Although these studies have advanced our understanding of eel reproductive biology, they have largely focused on different tissues or developmental stages. Research specifically targeting early gonadal development after sex differentiation in both male and female Japanese eels has not been thoroughly addressed. To fill this gap, we combined PacBio Iso-Seq and Illumina RNA-Seq technologies, leveraging their complementary strengths. While Illumina RNA-Seq provides highly accurate quantification of gene expression, its short-read lengths limit the ability to capture full-length transcripts. In contrast, PacBio Iso-Seq offers ultra-long reads, enabling the identification of complete isoforms and alternative splicing events. However, PacBio’s higher error rates and lower throughput are limitations[[110]34–[111]36]. By combining these technologies, we generated a comprehensive dataset that offers deeper insights into the transcriptome dynamics of early gonadal development in male and female Japanese eels. In this study, we explored gene regulatory elements in the gonads of Japanese eel, including TFs, SSRs, and lncRNAs. We identified numerous genes belonging to TF families, many from the HMG family. Members of the HMG family, particularly the SOX gene subgroup, are known to play significant roles in sex determination and gonadal development [[112]37, [113]38]. For instance, Sox3 knockout in female zebrafish results in delayed follicle development and reduced fertility [[114]39]. In Japanese flounder (Paralichthys olivaceus), Sox9a and Sox9b are highly expressed in the testes, indicating their crucial role in testicular development [[115]40]. Our findings revealed a substantial number of the HMG family members in the gonads of Japanese eels, indicating their critical role in gonadal function. We also predicted and compared SSRs in the male and female gonads of Japanese eels. SSRs, which are widely present in eukaryotic genomes, are linked to gene expression regulation, genome maintenance, and mutations [[116]41, [117]42]. We observed significant differences in SSR types and proportions between male and female gonads, indicating that further studies could investigate sex-related SSRs for sex identification in Japanese eels. Additionally, we detected a large number of lncRNAs expressed in the testes and ovaries of Japanese eels. LncRNAs have been reported to regulate gonadal function in fish. For example, a gonad-specific lncRNA in medaka (Oryzias latipes), termed lnc4, acts as a primary transcript for miR-202, which regulates androgen synthesis and spermatogenesis. The deletion of lnc4 led to impaired spermatogenesis and reduced testis size, indicating its crucial role in gonadal development [[118]43]. However, the roles of lncRNAs in eel gonadal development require further investigation. By comparing the DEGs between male and female Japanese eel gonads, we identified several genes involved in fundamental biological processes such as cell cycle regulation, metabolism, and reproductive development. Among the key genes verified in this study, many have established roles in reproductive pathways and hormone regulation. In the ovaries, bambi, gdf9, ccnb1, cdc20, and prlh were highly expressed. Bambi and gdf9 are enriched in TGF-beta signaling pathways. Bambi acts as a negative regulator of TGF-beta family signaling pathways by inhibiting the formation of receptor complexes. Its inhibitory role in signaling pathways is essential for embryogenesis and ovarian function [[119]44, [120]45]. In fish, its specific role requires further investigation. Gdf9 is crucial for ovarian development and oocyte maturation. Studies in Japanese flounder (Paralichthys olivaceus) and short-finned eel (Anguilla australis) revealed that gdf9 promotes steroidogenesis and facilitates communication between oocytes and somatic cells, ensuring proper follicular development and ovulation [[121]46, [122]47]. Cdc20 and ccnb1 are key regulators of the cell cycle, mitosis, and meiosis. Mutations in the cdc20 gene in humans lead to abnormalities in oocyte maturation and embryonic development, resulting in female infertility [[123]48, [124]49]. Additionally, studies in blunt snout bream (Megalobrama amblycephala) have shown that cdc20 is a maternal gene primarily expressed in oocytes, suggesting its involvement in reproduction and embryonic development[[125]50]. Ccnb1, as a key cell cycle regulatory protein, is involved in cell division and is essential for development [[126]51]. Studies indicate that the deletion of ccnb1 leads to intrauterine death in mice [[127]52]. While there are no current reports on the role of this gene in fish ovarian development and function maintenance, comparative transcriptomic studies between zebrafish ovaries and testes have identified ccnb1 and cdc20 as sex-biased genes with higher expression in ovaries [[128]53]. Our results in Japanese eels also confirm these findings. Further studies are needed to validate the specific role of these genes. Prlh has been shown to participate in the neuroendocrine regulation of fish reproduction. The addition of estradiol to primary cultured grass carp brain tissue significantly alters the expression levels of prolactin and its receptor. This suggests that prlh may be located on the brain-pituitary-gonad axis, playing a role in female reproductive development [[129]54]. We also verified five highly expressed genes in the testes of Japanese eels: ccdc39, syce3, tspo, chrebp, and ngb. Ccdc39 is related to sperm flagellar morphological abnormalities and motility disorders. It mainly localizes to the ciliary axoneme, and its deficiency affects the assembly of the dynein regulatory complex and normal ciliary motility [[130]50, [131]55, [132]56]. Studies in humans have reported that CCDC39 mutations cause severe sperm motility disorders due to abnormal sperm flagellar morphology, leading to male infertility [[133]57, [134]58]. Although there are no reports on the role of ccdc39 in fish gonadal development, our analysis suggests that it may also have an important function in fish spermatogenesis. Syce3, a central element of the synaptonemal complex, is crucial for meiosis [[135]59, [136]60]. Studies have shown that the absence of this gene leads to meiotic disruption and infertility in mice [[137]59]. Furthermore, overexpression and knockdown of syce3 affect the expression of steroidogenesis-related genes in Leydig and Sertoli cells, thereby regulating testosterone levels [[138]61]. This indicates that syce3 is a key gene in meiosis and reproduction and is involved in regulating testosterone synthesis. Tspo and chrebp are involved in steroid synthesis. Tspo, a transport protein, plays an important role in steroid biosynthesis by promoting the transport of cholesterol to the mitochondrial inner membrane, facilitating testosterone synthesis [[139]62, [140]63]. In Chrebp knockout homozygous mutant mice, plasma cholesterol levels were found to decrease [[141]64]. Since cholesterol is crucial for steroid hormone synthesis, it is speculated that Chrebp may indirectly participate in adrenal steroid hormone synthesis and release through this regulation. Our findings suggest that they may be involved in regulating testicular function through hormone regulation. Moreover, this study identified some hormone-related neuroglobin-like proteins, such as ngb. Although there is currently no literature explaining their role in gonadal development, GO annotations suggest they are related to developmental regulation, cell development processes, and reproduction. It is speculated that this gene may be functional in the testes of Japanese eels, and its specific role requires further verification. As with any study, there are limitations that should be acknowledged. In our analysis, we did not utilize the available reference genome for the Japanese eel. Integrating the latest reference genome could address some of the constraints in transcriptome-based studies and enhance future analyses. Conclusion This study provides the first comprehensive transcriptomic analysis of immature male and female Japanese eel gonads, utilizing both PacBio Iso-Seq and Illumina RNA-Seq technologies to generate a robust dataset, offering new insights into early post-differentiation gonadal development. Our findings identified numerous sex-biased genes involved in key processes such as cell cycle regulation, apoptosis, hormone activity, and reproductive development, alongside a range of TFs, SSRs, and lncRNAs that may contribute to gonadal differentiation. These results enhance our understanding of the molecular mechanisms governing gonadal development in Japanese eels and lay a valuable foundation for future functional studies aimed at improving reproductive regulation and breeding success in aquaculture. Supplementary Information [142]Supplementary Material 1^ (79.6KB, xls) [143]Supplementary Material 2^ (44.7KB, xls) [144]Supplementary Material 3^ (2.1MB, xls) [145]Supplementary Material 4^ (1MB, xls) [146]Supplementary Material 5^ (276.1KB, xls) [147]Supplementary Material 6^ (241.3KB, xls) [148]Supplementary Material 7^ (77.4KB, xlsx) [149]Supplementary Material 8^ (124.9KB, xls) [150]Supplementary Material 9^ (19.1MB, xls) [151]Supplementary Material 10^ (12KB, xlsx) Acknowledgements