Sexual dimorphism in growth of C. maculata

To notarize sexual dimorphism in growth, the body weight and length of males and females across multiple developmental stages were measured. Between 30 and 60 dpf, no significant differences in body weight and length were observed between males and females (P > 0.05). However, discernible differences in growth emerged from 90 dpf onwards (P < 0.05). At 90 dpf, male snakeheads (173.1 ± 29.3 g, 21.5 ± 1.3 cm) were 13.07% heavier and 5.64% longer than females (153.6 ± 25.3 g, 20.4 ± 1.2 cm), respectively. This sexual dimorphism in growth intensified over time, with males ultimately weighing 29.91% heavier and measuring 12.69% longer than females at 360 dpf, respectively (Fig. 1A, B). Notably, a significant positive correlation between body length and weight of male and female snakeheads was evident (Fig. 1C).

Fig. 1

Characterization of sexual dimorphism in growth in blotched snakehead. A Body weight of males and females from 30 to 360 dpf, B body length of males and females from 30 to 360 dpf, C correlation analysis between body length and weight in males and females

Histology and morphology of gonadal development in C. maculata

To elucidate the interplay between sexual dimorphism and gonadal development, a thorough histological and morphological examination of C. maculata gonads were conducted at different developmental stages. Initially, at 10 dpf, primordial gonads in females were observed in pairs beneath the mesonephric duct, containing 1–2 primordial germ cells (PGCs) distinguishable from somatic cells by their large diameter (Fig. 2A). Subsequently, at 15 dpf, somatic cells rapidly proliferated, leading to an increase in the volume of primordial gonads, accompanied by the appearance of blood vessels in females (Fig. 2B). As gonadal development progressed, sectional cell growth and organic heave manifestation occurred at 20 dpf (Fig. 2C), followed by the formation of an incipient ovarian cavity at 25 dpf, alongside further extension of the organic heave, indicating the onset of morphological sex differentiation in females (Fig. 2D). Oogonia was first observed in the germinal epithelium at 25 dpf (Fig. 2D). By 30 dpf, a fully developed ovarian cavity was evident in female gonads, with numerous oogonia undergoing mitosis (Fig. 2E). The chromatin-nucleolus oocytes and continued enlargement of the ovarian cavity were observed at 35 dpf (Fig. 2F). At 40 dpf, primary oocytes emerged in the ovary (Fig. 2G), which exhibited significantly larger volumes compared to oogonia, indicating the onset of cytological ovarian differentiation. Their numbers continued to increase by 45 dpf (Fig. 2H). Between 60 and 90 dpf, the ovary progressed to Stage II, due to the proliferation of primary oocytes and the appearance of growing oocytes (Fig. 2I, J). At 120 dpf, the ovary transitioned to Stage III, assuming an elliptical and cylindrical shape with visible circular eggs. Concurrently, a mass of growing oocytes were observed in the ovary (Fig. 2K). At 150 dpf, a plethora of growing oocytes filled the ovary, exhibiting significantly augmented volumes due to continuous synthesis and accumulation of nutrients, surrounded by two layers of follicular membranes (Fig. 2L). At this stage, oil droplets in growing oocytes increased markedly compared to primary oocytes, while yolk granules began appearing within the oil droplets, scattered in the cytoplasm. By 180 dpf, the ovary advanced to Stage IV, rapidly increasing in size and occupying most of the abdominal cavity. The surface of the ovarian membrane was decorated with blood vessels, and golden-yellow eggs became visible to the naked eye. Mature oocytes predominated, characterized by increased yolk granules and gradual disintegration of the nucleus (Fig. 2M, N).

Fig. 2

Histological analyses of ovaries and testes across different developmental stages in C. maculata. A–M Represent the ovary at 10, 15, 20, 25, 30, 35, 40, 45, 60, 90, 120, 150 and 180 dpf, respectively; a–m represent the testis at 10, 15, 20, 25, 30, 35, 40, 45, 60, 90, 120, 150 and 180 dpf, respectively. N and n show localized magnifications of the ovary and testis at 180 dpf, corresponding to M and m, respectively. PGC primordial germ cell, MD mesonephric duct, NO notochord, BV blood vessel, GC germ cell, OH organic heave, OC ovarian cavity, OG oogonia, SG spermatogonia, CNO chromatin-nucleolus oocyte, EDA efferent duct anlage, POC primary oocyte, PSC primary spermatocyte, GOC growing oocyte, SSC secondary spermatocyte, MOC mature oocyte, ST spermatids, SZ spermatozoa

In contrast, the onset of testicular differentiation in males occurred later than that in ovaries. During the early post-fertilization phase, male juveniles exhibited few PGCs in their primordial gonads, with most remaining inactive in mitosis (Fig. 2a). Between 15 and 30 dpf, both gonad size and germ cell numbers increased dramatically, although no evident morphological signs for testicular differentiation were observed (Fig. 2b–e). A distinct histological transition occurred at 35 dpf, marked by the appearance of the efferent duct anlage and limited round spermatogonia (Fig. 2f), indicating the initiation of morphological sex differentiation in the testis. Subsequently, the number of spermatogonia increased at 40 dpf (Fig. 2g), with primary spermatocytes becoming evident by 45 dpf, signifying the onset of cytological sex differentiation in the testis (Fig. 2h). Compared to spermatogonia, primary spermatocytes were smaller and had more intensely basophilic nuclei. The testis entered Stage II at 60 dpf, with primary spermatocytes arranged in bundles, forming seminiferous lobules (Fig. 2i). By 90 dpf, plentiful secondary spermatocytes were observed within seminiferous lobules, arranged in a regular radial pattern, with a central cavity appearing (Fig. 2j). The testis transitioned to Stage III at 120 dpf, characterized by being light-red, slightly enlarged in volume, and flat in shape, with histological sections displaying the simultaneous presence of spermatogonia, primary spermatocytes, and secondary spermatocytes (Fig. 2k). With ongoing testicular development, the testis volume significantly enlarged, with a darkening color and prominent large blood vessels on the surface, indicative of Stage IV. At 150 dpf, in addition to a few spermatogonia, primary spermatocytes, and secondary spermatocytes, a large number of spermatids appeared with spermatogenetic cysts, nearly filling the entire testis (Fig. 2l). Shortly thereafter, spermatozoa were observed in testis at 180 dpf (Fig. 2m, n).

Global transcriptome profiles of gonads in C. maculata

Comparative transcriptome analyses of gonads in female and male C. maculata were conducted across ten distinct developmental stages, spanning from 10 to 180 dpf. After removing low-quality reads, 273.66 Gb high-quality data were obtained from 40 samples (BioProject Accession No.PRJNA1087676), averaging 5.75 Gb per sample, with Q30 values ranging from 90.74% to 97.97%. The clean data were successfully mapped to the high-quality reference genome of C. maculata (SRA Accession No. PRJNA730430) [41]. However, during the correlation assessment and principal component analysis (PCA) of samples, one biological replicate (XX-90d-3) showed significant divergence from the other two replicates (XX-90d-1 and XX-90d-2) (Fig. S1A). Specifically, the correlation coefficients between sample XX-90d-3 and samples XX-90d-1 and XX-90d-2 were only 0.08 and 0.06, respectively (Fig. S1B). To ensure the accuracy of the results, sample XX-90d-3 was excluded from the subsequent analyses. In total, 33,266 transcripts were acquired, of which 24,115 corresponded to annotated protein-coding genes, while the remaining 9151 transcripts represented newly identified genes. Inter-sample correlation analysis and PCA revealed notable aggregation in the testes and ovaries transcriptomes before 30 dpf, followed by progressive divergence after 30 dpf (Fig. S1A, C).

Sex-related DEGs in male and female gonads: identification, spatiotemporal dynamics, functional annotation, and PPI analysis

A total of 72,330 DEGs were identified from pairwise comparisons between male and female samples at the same developmental stages (10, 15, 20, 25, 30, 60, 90, 120, 150, and 180 dpf), with 30,524 up-regulated and 41,806 down-regulated DEGs (Figs. S2, S3). Significantly, the number of up-regulated DEGs rapidly increased from 10 to 60 dpf, sharply decreased at 90 dpf, and then markedly rose again between 120 and 180 dpf. Likewise, down-regulated DEGs increased from 10 to 60 dpf, experienced a slight decline at 90 dpf, and then elevated again between 120 and 180 dpf (Fig. S3). Comparative analyses revealed notable up-regulation of genes such as Dmrt1, Amh, Amhr2, Star, Gsdf, Sox9b, and Ar in males compared to females, primarily participating in testicular differentiation and spermatogenesis. Conversely, DEGs associated with ovarian differentiation and oogenesis, including Foxl2, Cyp19a1a, Figla, and Bmp15, were significantly up-regulated in females relative to males. To comprehensively investigate the functional roles of DEGs between females and males at each time point, GO and KEGG enrichment analyses were conducted (Figs. S4, S5, Tab. S1–S4).

Subsequently, the expression profiles of 165 sex-related DEGs (Tab. S5, S6) were analyzed across ten developmental stages in females (Fig. 3A) and males (Fig. 3B), respectively. To elucidate the temporal dynamics of the transcriptomic datasets across these stages, the Mfuzz clustering method was used to categorize these 165 sex-biased DEGs into 12 clusters (Fig. 3C, D). Notably, the expression patterns observed in Cluster 5 and Cluster 8 of the female samples were particularly noteworthy, wherein pivotal transcription factors including Sox8, Sox9b, Sox11a, Sox17, Foxo3, Gata3, and Lhx8 decreased sharply from 10 dpf, surged again between 20 and 25 dpf, and then sharply declined and manifested low expression levels between 60 and 180 dpf, suggesting a potential role in early ovarian differentiation (Fig. 3C). Additionally, significant up-regulation in the transcription levels of meiosis marker genes, such as Sycp3 and Spo11, was observed in Clusters 3 and 4 during 30–60 dpf. This up-regulation is consistent with the observation of primary oocytes at 40 dpf, reinforcing the pivotal role of these genes in the onset of meiotic process. Key genes involved in follicular development in Cluster 3 (e.g., Figla, Wnt4, Zp4, Hsd17b12, Er) were not expressed from 10 to 30 dpf, but showed continuous expression from 60 to 180 dpf (Fig. 3C), revealing the crucial roles of these genes in female gonadal development. In contrast, major male-differentiation-related genes in Cluster 6 (Dmrt1, Dmrtb1, Amhr2, Star, Ar) were prominently expressed between 30 and 60 dpf (Fig. 3D), aligning with the commencement of morphological testicular differentiation at 35 dpf (Fig. 2f) and cytological testicular differentiation at 45 dpf (Fig. 2h).

Fig. 3

Heatmaps illustrate the expression levels of 165 sex-related DEGs in ovaries (A) and testes (B), respectively. The cluster analyses of these DEGs across ten key developmental stages in females (C) and males (D). UpSet network plots reveal the distribution of sex-related DEGs during crucial phases of male and female gonadal development (E). The upper bar chart shows the gene count in each group, and the lower left bar chart displays the number of sex-related DEGs at each gonadal development stage

To better understand the interactions among these 165 sex-related DEGs over ten developmental stages, an UpSet plot was constructed (Fig. 3E). Notably, Star exhibited pronounced male-biased expression spanning six developmental stages post 30–180 dpf, suggesting its potential significance in the regulation of testicular development. Additionally, Figla, Sox11b, Cyp19a1a, and Ctnd2 showed significant female-biased expression between 60 and 180 dpf, likely related to the development and maintenance of female gonads. Meanwhile, Sox3 and Sox11a exhibited significant male-biased expression prior to 30 dpf, suggesting their crucial role in early testicular differentiation.

The GO and KEGG functional enrichment analyses were performed on these 165 sex-related DEGs. Several sex-related GO terms were identified, including development of primary sexual characteristics (GO:0045137), sex differentiation (GO:0007548), female gonad development (GO:0008585), and female sex differentiation (GO:0046660) (Fig. 4A, Tab. S7). Additionally, significantly enriched KEGG pathways were observed, such as ECM-receptor interaction (ko04512), Cell cycle (ko04110), Oocyte meiosis (ko04114), and GnRH signaling pathway (ko04912) (Fig. 4B, Tab. S8). To further explore the interrelationships among selected DEGs, a PPI network was developed, encompassing pivotal sex-related genes like Dmrt1, Amh, Foxl2, Cyp19a1a, Wnt4, Er, and Figla. Within this network, 10 core genes, including Dmrt1, Amh, Cyp19a1a, Gsdf, Er, Figla, Ar, Fshr, Cyp11a1, and Cyp17a1, emerged as being highly interconnected (Fig. 4C).

Fig. 4

The GO (A), KEGG enrichment (B) and PPI network (C) analyses of 165 sex-related DEGs during gonadal development of the ovaries and testes

DEGs in female gonads: identification, spatiotemporal dynamics, functional annotation, and PPI analysis

Comparative analyses of ovarian tissues across nine female sample groups (10, 15, 20, 25, 60, 90, 120, 150, and 180 dpf) relative to 30 dpf ovary benchmark revealed a total of 96,716 DEGs, including 50,854 up-regulated and 45,862 down-regulated DEGs (Fig. 5A, S6A). These DEGs were classified into 12 distinct clusters, each showing discernible expression patterns, as illustrated in Fig. 5B.

Fig. 5

Heatmap illustrates the expression levels of all DEGs across ten key stages of ovarian development (A), cluster analysis of DEGs at ten critical periods in ovarian development using Mfuzz (B), Venn (C) and UpSet plots (D) show the distribution of DEGs at ten critical periods in ovarian development. The upper bar chart shows the gene count in each group, and the lower left bar chart displays the number of sex-related DEGs at each gonadal development stage

In Cluster 9, Sox17, Bmp2, Gnrh3, and Tfap2c exhibited high expression levels at 10 dpf, followed by a decline, reaching another peak at 25 dpf before rapidly decreasing (Fig. 5B). These DEGs are associated with the development of PGCs and female germ cells (FGCs). The findings suggest that the initial time of early ovarian differentiation occurs at 20–25 dpf or earlier, consistent with histological observations where the ovarian cavity and oogonia were first observed at 25 dpf (Fig. 2D). Well-known genes involved in FGCs development and ovarian differentiation, such as Dazl, DDx4, Sycp3, and Spo11 in Cluster 1, exhibited a sharp increase from 30 to 90 dpf, followed by a slight decrease in expression levels, but they remained relatively high. Similarly, other well-known genes involved in FGCs development and female gonadal differentiation, also observed in Cluster 12, like Nanog, Sall4, Sycp1, Sycp2, Figla, and Zar1, exhibited a sharp increase from 30 to 60 dpf, maintained high expression levels at 90 dpf, and subsequently declined gradually. Combining the first observation of primary oocytes at 40 dpf (Fig. 2G), it is inferred that molecular ovarian differentiation occurs between 40 and 60 dpf. Importantly, DEGs involved in estrogen synthesis were also identified in Cluster 1, such as Foxl2, Cyp19a1a, Zp3 and Zp4. In Cluster 5, the expression levels of Gdf9, Sox3, and Bmp15 continuously increased from 30 to 150 dpf and slightly decreased at 180 dpf, mainly involved in the follicle growth/vitellogenesis process. Additionally, DEGs involved in the response to steroid hormone were also identified, including Er, Pgrmc1, and Pgrmc2.

Based on UpSet and Venn diagram analyses (Fig. 5C, D), a total of 900 overlapping DEGs were identified across ten female developmental stages, encompassing critical genes involved in female gonadal differentiation and steroid hormone biosynthesis, notably including Cyp1a1, Hsd11b1, Cyp7b1, Akr1d1, Ugt1a1, Ugt2a3, Sult1a4 and Sts. Furthermore, DEGs at various developmental stages in females underwent analysis for GO and KEGG pathways. The pivotal GO terms associated with sex included sexual reproduction (GO:0019953), female gonad development (GO:0008585), and female gamete generation (GO:0007292), as shown in Fig. 6A and Tab. S9. The majority of identified KEGG pathways were linked to steroid hormone biosynthesis (ko00140), oocyte meiosis (ko04114), and ECM-receptor interaction (ko04512), as indicated in Fig. 6B and Table S10. PPI analysis highlighted central genes, such as Cyp19a1a, Er, Ccnb1, and Cdk1 (Fig. 6C).

Fig. 6

The GO (A), KEGG enrichment (B) and PPI network (C) analyses of all DEGs at ten key stages of ovarian development

DEGs in male gonads: identification, spatiotemporal dynamics, functional annotation, and PPI analysis

In male samples, a total of 51,681 DEGs were revealed through comparisons between 30 dpf and other developmental stages (10, 15, 20, 25, 60, 90, 120, 150, and 180 dpf), comprising 26,283 up-regulated and 25,398 down-regulated DEGs. Remarkably, there was a discernible upward trend in the number of DEGs as testicular development progressed, indicating an increasing involvement of genes in the regulatory processes of testicular differentiation and development (Fig. S6B).

Temporal trend analysis revealed 12 unique expression profiles for these DEGs in males. Marker genes present in PGCs and the undifferentiated spermatogonia, including Tfap2c, Tet1, Bmp2, Id4, Gnrh3 and Gfar1, showed significant up-regulation in Cluster 7 from 10 to 30 dpf, despite minor fluctuations in expression levels at certain intervals (Fig. 7A, B). Subsequently, their expression levels declined, remaining relatively low. Notably, the expression levels of Dmrt1 and Hormad1 in Cluster 9 were significantly increased within 30–60 dpf (Fig. 7B). Combined with the initial detection of the efferent duct anlage and spermatogonia at 35 dpf (Fig. 2f), this suggests a potential occurrence of early testicular differentiation between 30 and 35 dpf or earlier. The expression of numerous meiosis-related genes (e.g., Spo11, Sycp3, Spag6, Dmc1) in Cluster 6 sharply increased from 60 dpf, peaked at 90 days, then gradually declined, but remained at high levels until decreasing again at 150 dpf. Given the first observation of secondary spermatocytes at 90 dpf, it is indicated that 60–90 dpf may be the key time window for molecular testicular differentiation. Furthermore, DEGs involved in testicular development, such as Amh and Star, were also identified in Cluster 6. Notably, Ar and Cyp17a2 in Cluster 8 continued to rise (Fig. 7B), which are associated with responses to steroid hormone stimulation.

Fig. 7

Heatmap illustrates the expression levels of all DEGs across ten key stages of testicular development (A), cluster analysis of DEGs at ten critical periods in testicular development using Mfuzz (B), Venn (C) and UpSet plots (D) show the distribution of DEGs at ten critical periods in testicular development. The upper bar chart shows the gene count in each group, and the lower left bar chart displays the number of sex-related DEGs at each gonadal development stage

Utilizing UpSet and Venn diagram analyses (Fig. 7C, D), a consensus of 173 DEGs shared across ten male developmental stages were identified. This cohort of DEGs was implicated in steroid metabolism processes, including Star, Nrob2, Npc1l1, Ugt2a3, Apob, Malrd1 and Igf1. GO functional analysis revealed that terms related to male reproductive development were primarily enriched in categories such as male meiosis I (GO:0007141), male gamete generation (GO:0048232), and male sex differentiation (GO:0046661) (Fig. 8A, Table S11). Likewise, KEGG enrichment analysis identified key pathways such as ECM-receptor interaction (ko04512), MAPK signaling (ko04010), and TGF-β signaling (ko04350), as indicated in Fig. 8B and Table S12. To further explore the interactions among these DEGs, the PPI network was constructed, identifying six central genes: Dmrt1, Amh, Ccnd1, Nrob1, Bcl2a, and Wt1 (Fig. 8C).

Fig. 8

The GO (A), KEGG enrichment (B) and PPI network (C) analyses of all DEGs at ten key stages of testicular development

Subsequently, we identified a subset of 81 overlapping genes (Fig. S7A) shared between 900 overlapping DEGs in females (Fig. 5C) and 173 in males (Fig. 7C). These genes exhibited significantly biased expression in the 30 dpf ovaries, despite their lower expression levels at other developmental stages (Fig. S7B). GO enrichment analysis revealed that these 81 overlapping genes were predominantly enriched in the brush border (GO:0005903), clusters of actin-based cell projections (GO:0098862), response to lipopolysaccharide (GO:0032496), and transaminase activity (GO:0008483) (Fig. S7C). Furthermore, KEGG enrichment analysis identified several key metabolic pathways, including fat digestion and absorption (ko4975), cholesterol metabolism (ko4979), and phenylalanine, tyrosine, and tryptophan biosynthesis (ko00400) (Fig. S7D).

Validation and structural analyses of sex-related DEGs

To validate the accuracy of our transcriptome data, 18 sex-biased DEGs were randomly selected, and their expression patterns in male and female gonads across different developmental stages (10, 15, 20, 25, 30, 60, 90, 120, 150, and 180 dpf) were analyzed. The expression trajectories of these DEGs were visually depicted over time and space, along with their chromosomal locations. Figures 9 and 10 demonstrate that the expression trajectories from qPCR analyses closely aligned with those inferred from RNA-seq data. Our analysis revealed that Amhr2, Amh, Dmrt1, Cyp17a2, Ar, Gsdf, and Star were situated on chromosomes LG05, LG09, LG13, LG15, LG16, LG17, and LG20, respectively (Fig. S8), exhibited significantly higher expression in males (Fig. 9), indicating their crucial roles in male sex differentiation and spermatogenesis. Notably, Sox11a, located on sex chromosome LG02 [30], exhibited high expression in males before 30 dpf, which is the key period of early sex differentiation. Then, its expression decreased sharply at 60 dpf in both sexes. Subsequently, its expression in males was extremely low, while its expression in ovaries gradually increased after 90 dpf. Conversely, Er, Cyp19a1a, Ctnd1, Foxl2, Sox3, Figla, Ctnd2, Bmp15, Sox11b, and Pax4, located on chromosomes LG02, LG04, LG06, LG08, LG12, LG13, LG14, LG15, LG18, and LG21, respectively (Fig. S8), showed female-biased expression patterns (Fig. 10), suggesting their significance in female sex differentiation and oogenesis.

Fig. 9

qPCR analyses of critical male-biased DEGs during gonadal development. A Amh, B Amhr2, C Ar, D Cyp17a2, E Dmrt1, F Gsdf, G Star and H Sox11a. The results of RNA-seq and qPCR in each panel were shown in the left and right, respectively

Fig. 10

qPCR analysis of critical female-biased DEGs during gonadal development. A Bmp15, B Ctnd1, C Ctnd2, D Cyp19a1a, E Er, F Figla, G Foxl2, H Pax4, I Sox11b, J Sox3. The results of RNA-seq and qPCR in each panel were shown in the left and right, respectively

The levels of sex steroid hormone in serum

The fish ELISA kit was utilized to ascertain the levels of E2 and T during gonadal differentiation and development in male and female snakeheads. E2 levels exhibited pronounced sexual dimorphism, being significantly more abundant in females, suggesting that E2 may play a pivotal role in female gonadal differentiation and maintenance (Fig. 11A). Likewise, T levels in males were significantly higher than those in females, highlighting the potential critical role of T in male sex differentiation and maintenance (Fig. 11B).

Fig. 11

Comparison of E2 (A) and T (B) levels across six developmental stages in female and male snakeheads. Asterisks represent significant differences between males and females. *P < 0.05, and **P < 0.01