Gonadal development is driven by a complex genetic cascade in vertebrates. However, related information remains limited in molluscs owing to the long generation time and the difficulty in maintaining whole life cycle in the lab. The dwarf surfclam Mulinia lateralis is considered an ideal bivalve model due to the short generation time and ease to breed in the lab.
To gain a comprehensive understanding of gonadal development in M. lateralis, we conducted a combined morphological and molecular analysis on the gonads of 30 to 60 dpf. Morphological analysis showed that gonad formation and sex differentiation occur at 35 and 40–45 dpf, respectively; then the gonads go through gametogenic cycle. Gene co-expression network analysis on 40 transcriptomes of 35–60 dpf gonads identifies seven gonadal development-related modules, including two gonad-forming modules (M6, M7), three sex-specific modules (M14, M12, M11), and two sexually shared modules (M15, M13). The modules participate in different biological processes, such as cell communication, glycan biosynthesis, cell cycle, and ribosome biogenesis. Several hub transcription factors including SOX2, FOXZ, HSFY, FOXL2 and HES1 are identified. The expression of top hub genes from sex-specific modules suggests molecular sex differentiation (35 dpf) occurs earlier than morphological sex differentiation (40–45 dpf).
This study provides a deep insight into the molecular basis of gonad formation, sex differentiation and gametogenesis in M. lateralis, which will contribute to a comprehensive understanding of the reproductive regulation network in molluscs.
Mulinia lateralis reaches sexual maturation within 2 months, facilitating comprehensive studies of gonadal development.
A combined morphological and molecular analysis reveals the timing of gonad formation, sex differentiation and gametogenesis.
WGCNA on 40 transcriptomes identifies key modules and molecular drivers of clam gonadal development.
SOX2, FOXZ, HSFY, FOXL2 and HES1 are hub transcription factors for gonadal development of Mulinia lateralis.
Gonadal development is a complex process consisting of gonad formation, sex differentiation and gametogenesis. This process requires the coordinated expression of a specific set of genes in a strict spatiotemporal manner. In mammals, stem cell-like primordial germ cells migrate to the genital ridge, which gives rise to the bipotential gonad. In this period, GATA4, EMX2, LHX9, and NR5A1 (also known as SF1) are crucial for the formation of genital ridge [1,2,3,4], while SOX2, BMP2/4, OCT4, SOX17, and PRDM1/14 are involved in early germ cell specification [5,6,7,8,9,10,11]. Sex differentiation occurs under the control of Y-linked transcription factor SRY. It activates SOX9 to drive the testis differentiation [12, 13]. In the absence of SRY, female pathways involving RSPO1/WNT4/β-catenin and FOXL2 are induced, giving rise to the ovaries [14,15,16]. Finally, gonads undergo proliferation and meiosis to produce mature oocyte or spermatozoa . Many genes are involved in this process, and disruption of the complex mechanisms regulating oogenesis/spermatogenesis may cause infertility.
Mollusca, including bivalves, gastropods, cephalopods, composes the second-largest phylum of animal kingdom after arthropods, with around 200,000 extant species widespread in marine, fresh water and on land . Many molluscs, particularly bivalves such as oysters, scallops, clams and mussels, are an important food source for humans. As economically important aquaculture shellfish, molluscs (17.3 million tons) represented 56.3% of the production of marine and coastal aquaculture in 2018 . They have also gained attention as model species for climate change studies . Research on the gonadal development of molluscs not only helps to understand the molecular mechanism of sex differentiation and gametogenesis, but also assists in unraveling the impact of temperature changes on reproduction.
Molluscs exhibit a diversity of sexual systems, including gonochorism, simultaneous hermaphroditism, and sequential hermaphroditism [21, 22]. Although intensive efforts have been made to elucidate the key genes and pathways regulating gonadal development, most studies focus on germline specification and sexually differentiated gonads. According to previous reports, NANOS, VASA and PIWI are critically important for germ cell development [23,24,25], and FOXL2, DMRT1L, FEM-1-like, SOXH, VIT and TSSK1 play pivotal roles in sex maintenance or gametogenesis [26,27,28,29]. Till now, a comprehensive knowledge on the molecular mechanisms of gonadal development including gonad formation, sex determination/differentiation and the transition of gametogenic cycle is still lacking in molluscs. The main obstacle is that most molluscs are difficult to maintain in the lab throughout their lives mainly due to two reasons: (1) many of them have a biphasic life cycle and undergo metamorphosis, which requires special settlement cues that are difficult to obtain; (2) most species have long generation times (generally > 1 year), and cannot be raised under laboratory conditions throughout their lives.
The dwarf surfclam Mulinia lateralis is a small marine bivalve native to the Gulf of Mexico, West Indies, and the Atlantic coast . It is considered an ideal mollusc model for genetic research because it has a short generation time (2–3 months) and is easy to maintain and breed in the lab . Since being introduced to China in 2017, M. lateralis has been successfully cultured in our lab . Its gonad reaches maturation in 2 months, making it a good material for the study of gonadal development. In present study, we collected the dwarf surfclam from 30 days post-fertilization (dpf) to 60 dpf, and investigated their gonadal development from the morphological and molecular perspectives. It will contribute to a better understanding on the molecular basis of gonad formation, sex differentiation and gametogenesis in M. lateralis, and potentially useful for other molluscs.
Materials and methods
All experimental dwarf surfclams were cultured in the recirculating aquaculture system in our lab. At least 15 individuals were collected every 5 days from 30 dpf until 60 dpf. Gonads were carefully excised under a stereoscope. Part of them were immediately frozen in liquid nitrogen and then transferred to − 80 °C for RNA extraction, and the rest were fixed in Bouin’s solution overnight at room temperature for histological analysis. Since the gonads cannot be observed until 35 dpf, the molecular studies were initiated at this time point. The testicular and ovarian samples are indicated by the suffix “T” or “O”.
After the overnight fixation in Bouin’s solution, the specimens were dehydrated in graded ethanol, cleared in xylene, and embedded in paraffin. Tissue blocks were then cut into 5-µm serial sections on a rotary microtome (Leica, Wetzlar, Germany), and incubated overnight at 37 °C. Finally, the sections were stained with hematoxylin and eosin, and observed under a Nikon’s Eclipse E600 research microscope.
RNA library construction and sequencing
Total RNA was extracted using the conventional guanidinium isothiocyanate method, followed by DNase I (Takara Bio, Shiga, Japan) digestion to eliminate potential DNA contamination. RNA concentration and purity were determined by Nanovue Plus spectrophotometer (GE Healthcare, New Jersey, USA), and the quality was assessed by agarose gel electrophoresis.
High-quality RNA was used for RNA-Seq library construction following the instructions of the VAHTS Universal V6 RNA-seq Library Prep Kit for Illumina (Vazyme, Nanjing, China). Forty libraries (3–6 biological replicates for each sex at each time point) were subjected to 150 bp paired-end sequencing on the Illumina NovaSeq. All data have been submitted to the NCBI (BioProject ID: PRJNA862073).
Data preprocessing and normalization
To obtain high-quality (HQ) reads, adapters were removed with Trimmomatic-0.35, and homemade Perl scripts were used to remove reads containing undetermined bases (‘N’) or excessive numbers of low-quality positions (> 10 positions with quality scores < 20). The HQ reads were then mapped to the M. lateralis genome (unpublished) using STAR  with the parameters of ‘STAR -runThreadN 10 -genomeDir fasta -readFilesIn R1.fq R2.fq -sjdbGTFtagExonParentGene gff3 -outFilterMultimapNmax 100,000 -outFileNamePrefix -outReadsUnmapped fasta > STAR_sample.log’. Finally, the raw counts for each gene were obtained by Htseq-count with the parameters of ‘htseq-count -f bam -r name -s reverse -mode = intersection-nonempty -i Parent sample.sorted.bam gff3 > sample.count’, and converted into transcripts per million (TPM) using the RNA-Seq by Expectation Maximization (RSEM).
Gene co-expression network construction
A signed network was constructed using the R package weighted gene co-expression network analysis (WGCNA) following the procedures described . Briefly, a power of 15 was chosen for the network to exhibit approximate scale-free topology (model fitting index R2 = 0.75). Subsequently, an adjacency matrix was generated by calculating Spearman correlation coefficient between gene pairs. All genes were hierarchically clustered based on dissimilarity measure of topological overlap , and the resulting gene dendrogram was used for module detection using the cutreeDynamic method (minModuleSize = 100 and detectCutHeight = 0.99).
Identification of gonadal development-related modules
To identify gonadal development-related modules, we first identified gonadal development-related genes, which are expected to be differentially expressed among time points or between sex. Here, one-way analysis of variance (ANOVA) was conducted for each gene, and false discovery rate (FDR) was calculated using the q-value package to account for multiple tests. Only genes with q-value < 0.01 were considered as differentially expressed.
Gonadal development-related modules were determined by overrepresentation analysis of differentially expressed genes (DEGs) using a hypergeometric test, and modules with q-value < 0.1 were considered as gonadal development-related modules.
Functional annotation of genes and modules
Functional annotation of the genes was carried out using BLASTP (E-value cutoff 1e−5) against the Swiss-Prot database, and the results were imported into Blast2GO. Functional annotation of a module was conducted by GO and KEGG enrichment analyses using GOstats , and terms or pathways with adjusted q-value < 0.1 were considered as significantly more enriched than random chance in the particular module.
Hub gene selection and visualization
Hub genes are highly connected genes in a module. A gene’s connection strength to other genes can be indicated by intramodular connectivity Kin. Therefore, top 15% genes with high Kin were considered as the hub genes for a given module. Co-expression patterns of top 50 hub genes were visualized using the heatmap.
Histological analysis of the gonads
The anatomy of Mulinia lateralis is shown in Fig. 1A. Considering that M. lateralis reaches sexual maturation within 2 months, we examined the gonadal development procedures using juveniles aged 30 (shell height 2.27 ± 0.20 mm) to 60 dpf (shell height 7.70 ± 0.59 mm) (Fig. 1B). As shown in Fig. 1C, germ stem cells (GSCs) are surrounded by a monolayer of follicle cells near the crystalline style at 30 dpf, and the gonads can be observed between digestive gland and foot at 35 dpf; more follicles appear at 40 dpf, each containing sexually indistinguishable gonia. From 45 to 60 dpf, sex can be determined based on the histological analysis (Fig. 1D, E). In the testes, a small number of spermatocytes appear at 45 dpf (shell height 4.92 ± 0.46 mm). Subsequently, the follicles grow bigger, and various types of germ cells are visible at 55 dpf including spermatogonia, spermatocytes, spermatids, and spermatozoa; at 60 dpf, the follicles are predominantly filled with spermatids and spermatozoa (Fig. 1D). In the ovaries, oocytes show up in some follicles at 45 dpf, and more oocytes are observed with the growth of follicles from 50 to 60 dpf. Mature oocytes are visible at 55 and 60 dpf (Fig. 1E).
The histological analysis showed that in M. lateralis, gonad formation and sex differentiation occur at 35 and 40–45 dpf, respectively. Subsequently, the gonads go through gametogenic cycle, including a proliferating stage at 45–50 dpf, a growing stage at 50–55 dpf and a mature stage at 55–60 dpf.
Gonadal transcriptome sequencing and analysis
To better understand the molecular mechanisms of gonadal development in M. lateralis, we constructed 40 RNA-Seq libraries using 35–60 dpf gonads (Table 1). A total of 917,474,183 raw reads of 150 bp are obtained, of which 873,981,969 (95.26%) clean reads remained after adapter removal. Among them, 835,574,378 (95.61%) HQ reads are screened by homemade Perl scripts.
To obtain the relationship between the 40 transcriptomes, we performed a principal component analysis (PCA) on the 18,074 genes expressed in at least one group. According to Fig. 2A, the samples display a trend of differentiation from sexually indistinguishable gonads (35 and 40 dpf) to the mature testes and ovaries. The first principal component (PC1), which explains 55.3% of the variation, discriminates most of the testes from the ovaries and sexually indistinguishable gonads. The gonadal development stages are organized along the second principal component (PC2).
Network construction and identification of gonadal development-related modules
Gene co-expression network enables identification of modules that represent functional categories [37,38,39]. Here, 18,074 genes from the 40 gonadal transcriptomes were used for network construction, resulting in 15 modules (M1 ~ M15) with module size ranging from 103 to 2350 genes (Fig. 2B).
A total of 6986 DEGs (q-value < 0.01) are identified by differential expression analysis, of which 6768 (~ 97%) are included in the network construction. Overrepresentation analysis revealed that seven modules (q-value < 0.1) including M15 (q-value = 6.53e−202), M7 (q-value = 2.54e−185), M12 (q-value = 8.84e−57), M14 (q-value = 7.07e−40), M11 (q-value = 1.42e−27), M13 (q-value = 3.95e−3) and M6 (q-value = 0.07) are significantly enriched with DEGs (Fig. 2C). The proportion of DEGs is greater than 40% in each of these modules.
To determine to which gonadal development process the DEGs enriched modules are related, we investigated the gene expression profiles of top 50 hub genes in these modules (Additional file 1: Table S1). Results showed that the top 50 hub genes of M6 and M7 modules are highly expressed during early stage at 35 and 40 dpf, and for M7, the expression continues to the early differentiated testes and ovaries (Fig. 3A). Hub genes of M14, M12 and M11 display sex-specific expression patterns, with M14 and M12 being testis specific, and M11 being ovary specific (Fig. 4A). Hub genes of M15 and M13 are expressed at all stages of gonads, but the expression level is higher in the mature gonads than in the early ones (Fig. 5A). According to the expression patterns, these modules are considered as gonad-forming modules (M6, M7), sex-specific modules (M14, M12, M11), and sexually shared modules (M15, M13), respectively.
Functional annotation and hub genes of gonad-forming modules
To further understand the biological roles of these gonadal development-related modules, GO and KEGG enrichment analyses were performed. The gonad-forming module M6 is enriched with genes involved in cell communication (q-value = 5.03e−9), nuclear receptors (q-value = 1.11e−4), transcription factors (q-value = 4.22e−3) and signaling pathways regulating pluripotency of stem cells (q-value = 8.10e−3) (Fig. 3B). Transcription factor SOX2 has the highest intramodular connectivity, suggesting it may be a key driver for early gonadal development. The other genes such as transcription factors SOX9, GATA2 (GATA-binding factor 2), OVO (transcriptional regulator ovo), nuclear receptors NR2F1 (nuclear receptor subfamily 2 group F member 1), NR1D3 (nuclear receptor subfamily 1 group D member 3) and ER (estrogen receptor), cell communication-related genes WNT7B (protein Wnt-7b precursor) and HTR4 (5-hydroxytryptamine receptor 4) are among the top 15% hub genes.
The other gonad-forming module M7 is enriched for G protein-coupled receptors (q-value = 1.50e−13), cell adhesion molecules (q-value = 2.41e−13), GnRH signaling pathway (q-value = 1.28e−5) and ECM-receptor interaction (q-value = 4.02e−5) (Fig. 3B). The top 15% hub genes include G protein-coupled receptors 5-HTR1 (5-hydroxytryptamine receptor 1), FZD10 (Frizzled-10) and SSTR2 (Somatostatin receptor type 2), cell adhesion molecules NRXN3 (Neurexin-3) and NLGN4Y (Neuroligin-4, Y-linked), Rut (Ca(2+)/calmodulin-responsive adenylate cyclase) and GNRHR (gonadotropin-releasing hormone receptor) from GnRH signaling pathway, and several collagens (COL4A5, COL6A6 and COL1A2).
Figure 3C shows the dynamic expression profiles of some representative genes of M6 and M7 during gonadal development. Except for the above-mentioned hub genes, we found some previously reported early gonadal development-related genes, such as GATA4 (GATA-binding factor 4) and SIX4 (homeobox protein SIX4), are also in module M6 and M7, and display similar expression pattern with the hub genes.
Functional annotation and hub genes of sex-specific modules
The ovary-specific module M11 is enriched with genes participating in N-glycan biosynthesis (q-value = 1.25e−3), biosynthetic process (q-value = 1.04e−2), and fatty acid biosynthesis (q-value = 6.71e−2) (Fig. 4B). Top 15% hub genes including HES1 (transcription factor HES-1), GAL3ST3 (Galactose-3-O-sulfotransferase 3), FOXL2 (forkhead box protein L2) from biosynthetic process, FASN (Fatty acid synthase) from fatty acid biosynthesis, RPN1 (26S proteasome regulatory subunit RPN1) and DDOST (dolichyl-diphosphooligosaccharide–protein glycosyltransferase 48 kDa subunit) from N-glycan biosynthesis are specifically expressed in the ovaries. Similar expression patterns are also detected for VIT6 (Vitellogenin-6), a female-specific glycoprotein, as well as some cell cycle genes, such as MP62 (mitotic apparatus protein p62) and CDC20 (cell division cycle protein 20) (Fig. 4C).
The testis-specific modules are enriched with spermatogenesis-related categories. Specifically, M14 is overrepresented with genes involved in meiotic cell cycle (q-value = 7.63e−4) and reproductive process (q-value = 2.95e−3) (Fig. 4B). Representative genes such as transcription factors FOXZ (Forkhead box protein Z), HSFY (heat shock transcription factor, Y-linked) and SOXH (transcription factor SOXH), chromosome and associated proteins including tubulins, PRDM9 (histone-lysine N-methyltransferase PRDM9) and SYCP3 (Synaptonemal complex protein 3), DNA repair and recombination proteins including HFM1 (Probable ATP-dependent DNA helicase HFM1) and MSH4 (MutS protein homolog 4), display similar expression patterns, with specific expression in the testes and especially high in later stages. Spermiogenesis related genes, such as testis-specific serine/threonine-protein kinase family genes (TSSK1, TSSK4) are also present in module M14 and highly expressed in the testes (55 and 60 dpf) (Fig. 4C).
The other testis-specific module M12 is enriched for cilium and associated proteins (q-value = 1.42e−225), oxidative phosphorylation (q-value = 2.16e−7), and protein kinases (q-value = 2.46e−6) (Fig. 4B). Representative hub genes from these terms include cilia- and flagella-associated proteins, axonemal dyneins and radial spoke head proteins, which are highly expressed in 60 dpf testes and moderately in 55 dpf testes and 35 dpf gonads (Fig. 4C).
Functional annotation and hub genes of sexually shared modules
The two sexually shared modules M15 and M13 are enriched with biological processes commonly exist in spermatogenesis and oogenesis. M15 is enriched with genes participating in cell cycle (q-value = 3.60e−24), DNA replication (q-value = 1.05e−30) and repair (q-value = 8.44e−50) (Fig. 5B). Representative hub genes such as SMC5 (structural maintenance of chromosomes protein 5), MRE11 (double-strand break repair protein MRE11), CDK1 (cyclin-dependent kinase 1) and RFC1 (replication factor C subunit 1) are primarily expressed in the late stages of testes and ovaries, with the peak in 60 dpf testes (Fig. 5C).
Module M13 is enriched for ribosome biogenesis (q-value = 4.75e−219), translation (q-value = 8.32e−73), and spliceosome (q-value = 3.94e−31) (Fig. 5B). Hub genes such as NOP14 (nucleolar protein 14), DDX51 (ATP-dependent RNA helicase DDX51), DDX27 (Probable ATP-dependent RNA helicase DDX27), DDX54 (ATP-dependent RNA helicase DDX54), EIF3C (eukaryotic translation initiation factor 3 subunit C), PRPF3 (U4/U6 small nuclear ribonucleoprotein Prp3) and PABP1 (polyadenylate-binding protein 1) are highly expressed in the late stages of testes and ovaries, with the peak in 60 dpf ovaries (Fig. 5C).
Timing of molecular sex differentiation
To determine the timing of molecular sex differentiation in M. lateralis, we selected sex-specific genes for hierarchical clustering analysis and PCA on the 40 samples. Considering that M11 and M14 are sex-specific modules with genes consistently highly expressed from 45–60 dpf, we chose top 40 hub genes that show differential expression (p-value < 0.05) between genders at 45 dpf from each module. As shown in Fig. 6, similar results were obtained for hierarchical clustering and PCA analyses, showing that the 40 samples can be clustered into three groups (ovaries, testes and undifferentiated gonads). The genders of dwarf surfclams at 45–60 dpf are easily distinguished, consistent with the morphological results. In addition, the sex of most of the sexually indistinguishable gonads at 35 and 40 dpf can also be determined, with the ratio of 67% (4/6) and 83% (5/6), respectively.
Summary of the developmental process of M. lateralis gonads
Based on the above results, the gonadal developmental process of M. lateralis is summarized at histological and molecular levels (Fig. 7). Gonad formation occurs as early as 35 dpf, in which some cell communication-related genes, G protein-coupled receptors, and transcription factors are involved, such as SOX2, SOX9, GATA2/4, OVO, WNT7B and FZD10. At 35–40 dpf, molecular sex differentiation occurs, with the elevation of HSFY and FOXZ in the testis, and the induction of FOXL2 and HES1 in the ovary. At 45 dpf, sex can be distinguished based on histological observation, and the gonads start to go through gametogenic cycle. Cell cycle and ribosome biogenesis-related genes are highly expressed in the ovary and testis during the gametogenesis, such as SMC5, CDK1, DDX51 and DDX27. But chromosome and cilium-related genes such as PRDM9, SYCP3, SOXH and CFAP157 are specifically expressed in the testis, participating in the spermatogenesis. In the ovary, N-glycan and fatty acid biosynthesis related genes such as RPN1, DDOST, VIT6 and FASN are highly expressed, as well as some meiosis arrest-related genes like MP62, CDC20, CCNO-B and CCNB2, which could play vital roles in oocyte development.
Gonadal development involves various biological processes, such as gonad formation, sex differentiation, and gametogenesis. Comprehensive studies of gonadal development require continuous sampling until sexual maturation. For example, in Patinopecten yessoensis , Argopecten irradians  and Crassostrea sikamea , individuals aged 5–13, 2–8 and 1–15 months old were tracked to illustrate sex differentiation and gametogenesis, respectively. Long generation time (usually 1–2 years)  limits the studies of gonadal development for most molluscs. In contrast, Mulinia lateralis has been reported to reach sexual maturation within 2–3 months [31, 43], which facilitates comprehensive studies of gonadal development. Our study confirms that M. lateralis become mature at 60 dpf. Most importantly, based on the histological analysis of gonads from 30–60 dpf every 5 days, we determined the timing of gonad formation, sex differentiation and gametogenesis. It provides a valuable resource for studies of the molecular mechanisms of gonadal development, and will facilitate more precise regulation of M. lateralis reproduction in the lab.
Transcriptome analysis has been widely used to elucidate the molecular mechanisms underlying various traits. However, it often produces many DEGs, and the key drivers are difficult to identify. WGCNA represents a transition from gene-centric to network or module-centric approaches. It enables researchers to identify important modules and central players within modules, and has proven its power in many organisms including molluscs [39, 44,45,46,47]. In this study, a gene co-expression network was constructed based on 40 transcriptomes covering the entire gonadal development process of M. lateralis, allowing a comprehensive understanding of the molecular changes and key drivers for gonad formation, sex differentiation and gametogenesis.
Among the 7 gonadal development-related modules, the three sex-specific modules are enriched with GO terms or KEGG pathways that have been reported in previous studies [26, 28, 48,49,50,51]. Specifically, enrichment of N-glycan biosynthesis, fatty acid biosynthesis and GABAergic synapse in the ovary is also found in P. yessoensis [26, 48], Pinctada margaritifera  and Chlamys farreri , suggesting these biological processes may be critical for ovary development in various mollusc species. Cell cycle, sexual reproduction and oxidative phosphorylation have also been reported to be overrepresented in the testis-biased genes in P. yessoensis , Nodipecten subnodosus  and Crassostrea gigas . According to the functional annotation and key genes, the two testis-specific modules of M. lateralis correspond to the turquoise and green modules of P. yessoensis , suggesting the biological processes involved in the two modules are widely present in different molluscs.
In addition to sex-specific modules, we identified two sexually shared modules (M15 and M13). These modules participate in general biological processes such as DNA replication, repair, and translation, which occur in the testes and ovaries. Although similar modules or functional categories have been reported in other molluscs, they are generally considered to be sex biased [48, 51]. Specifically, DNA replication and repair are enriched in testis-biased genes, likely because there are more cells going through cell cycles in the testes than ovaries. Ribosome biogenesis and translation are enriched in ovary-biased genes, possibly because these processes occur more frequently in the oocytes than the male gametes. The sexual-biased expression patterns may cause these biological processes to be overlooked in the other sex and should receive more attention in future studies.
As is discussed above, previous studies of mollusc gonad development primarily focused on the gametogenic cycle and sexual differentiation [52,53,54,55], whereas studies on gonad formation are lacking. Since consecutive gonadal samples can be obtained for M. lateralis in the early stages, they can be used to construct gene network and identify gonad-forming modules. The two modules we identified consist of genes that show highest expression in 35 dpf gonads, which is quite different from the expression patterns of sex-specific or sexually shared modules. Functional annotation suggests involvement of cell communication, hormone-mediated signaling pathway, G protein-coupled receptors and ECM–receptor interaction in gonad formation. Similarly, several wnt ligands and their receptors Fzs from Wnt signaling pathway have been demonstrated to play a vital role in the stemness regulation of mammalian spermatogonial stem cells . In Xenopus laevis, signaling molecules (e.g., WNT7B), extracellular matrix molecules (e.g., collagens) and transcription factors (e.g., LHX8) also play important roles at the initial stage of gonadogenesis . These studies suggest a conserved role of wnt signaling pathway and ECM–receptor interaction in gonadogenesis in different animals.
Transcription factors are well known to play critical roles in diverse biological processes, including gonadogenesis [4, 58] and sexual development [16, 59]. Our study reveals that transcription factors SOX2, SOX9, GATA2 and OVO are important for gonad formation, FOXL2 and HES1 for ovary development, and HSFY, FOXZ and SOXH for testis development in M. lateralis. Among them, SOX2, FOXL2, FOXZ and SOXH have been reported to be vital for gonadal development in other molluscs [40, 41, 60,61,62,63,64]. Studies on FOXL2 have demonstrated its central role in sex differentiation and ovary maintenance in several mollusc species [40, 41, 60,61,62]. The other genes (SOX2, SOXH and FOXZ) are reported to be essential for spermatogenesis and testis development in scallops and oysters [60, 63, 64]. Except for the protostome-specific FOXZ, the remaining three genes seem to play conserved roles in vertebrates. Specifically, FOXL2 also exerts critical roles in sex determination , sex differentiation , and gametogenesis . SOX2 is essential for PGC (primordial germ cell) proliferation , and SOXH is associated with testis development and required for male fertility [67, 68]. The remaining transcription factors are rarely studied in molluscs, but some of them have been demonstrated to be involved in the sexual development of other animals. For example, the roles of SOX9 and HSFY in testis development have been well studied in mammals. SOX9 is the target of SRY, which mediates testis determination and differentiation [69, 70]. HSFY is a Y-linked gene in human whose deletion results in male infertility, and is also involved in the spermatogenesis in other animals [71,72,73,74]. OVO, HES1 and GATA2 are also reported to play essential roles in ovary development in various organisms. For example, OVO is important for female germline survival and oogenesis in insects [75, 76]. HES1 is the major target genes of NOTCH signaling, which regulates ovarian somatic cell development and is necessary for oocyte maturation in mammals [17, 77]. GATA2 is transiently expressed in the ovarian germ cell lineage during mouse gonadogenesis, and the signal location coincides with germ cell marker OCT4 . Some of the important transcription factors (e.g., FOXL2, SOX9 and SOXH) we found in M. lateralis are also identified as key regulators of sex development in mouse , tilapia  and sea bass  by WGCNA, suggesting the drivers of gonadal development may be conserved across animal phyla.
The onset of sex differentiation can be determined by the morphological or molecular method [80, 82]. In this study, we found molecular sex differentiation occurs as early as 35 dpf, while morphological sex differentiation occurs at 40–45 dpf in M. lateralis. It suggests sex differentiation occurs earlier at the molecular level than at the morphological level, similar to the findings in scallops P. yessoensis  and A. irradians . But unlike scallops in which FOXL2 and DMRT1L are the key sex differentiation genes, we found that only FOXL2 shows sexual dimorphic expression in M. lateralis. This confirms that the role of FOXL2 in the ovary development is likely to be deeply conserved in molluscs. DM domain genes play central roles in sex differentiation in many metazoans, including fish, birds, nematodes and arthropods . However, it may not be involved in sex differentiation of M. lateralis. We speculate that other transcription factors, such as FOXZ and HSFY, may replace DMRT1L to promote testis development.
Perspectives and significance
Mollusca composes the second-largest phylum of animal kingdom after arthropods, and is an important food source for humans. Our study on the gonadal development of dwarf surfclam M. lateralis, which is considered an ideal bivalve model for genetic research, provides detailed information on the gene network of gonad formation, sex differentiation and gametogenesis. It will facilitate reproductive control in molluscs. The hub transcription factors we identified could be key drivers in mollusc gonadogenesis, sex differentiation or gamete development. In the future, comprehensive studies on the expression, regulation and function of these transcription factors are needed to elucidate the molecular mechanism of mollusc sexual development and to reveal its consistency/differences with other animal phyla.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its additional information files.
Kranz AM, Tollenaere A, Norris BJ, et al. Identifying the germline in an equally cleaving mollusc: Vasa and Nanos expression during embryonic and larval development of the vetigastropod Haliotis asinina. J Exp Zool Part B Mol Dev Evol. 2010;314(4):267–79.
Li Y, Zhang L, Sun Y, et al. Transcriptome sequencing and comparative analysis of ovary and testis identifies potential key sex-related genes and pathways in scallop Patinopecten yessoensis. Mar Biotechnol. 2016;18(4):453–65.
Li J, Zhou Y, Zhou Z, et al. Comparative transcriptome analysis of three gonadal development stages reveals potential genes involved in gametogenesis of the fluted giant clam (Tridacna squamosa). BMC Genomics. 2020;21(1):872.
Teaniniuraitemoana V, Huvet A, Levy P, et al. Gonad transcriptome analysis of pearl oyster Pinctada margaritifera: identification of potential sex differentiation and sex determining genes. BMC Genomics. 2014;15(1):491.
Fu X, Sun Y, Wang J, et al. Sequencing-based gene network analysis provides a core set of gene resource for understanding thermal adaptation in Zhikong scallop Chlamys farreri. Mol Ecol Resour. 2014;14:184–98.
Xu W, Cui Z, Wang N, et al. Transcriptomic analysis revealed gene expression profiles during the sex differentiation of Chinese tongue sole (Cynoglossus semilaevis). Comp Biochem Physiol D: Genomics Proteomics. 2021;40:100919.
Xu R, Pan L, Yang Y, et al. Characterizing transcriptome in female scallop Chlamys farreri provides new insights into the molecular mechanisms of reproductive regulation during ovarian development and spawn. Gene. 2020;758:144967.
Llera-Herrera R, García-Gasca A, Abreu-Goodger C, et al. Identification of male gametogenesis expressed genes from the scallop Nodipecten subnodosus by suppressive subtraction hybridization and pyrosequencng. PLoS ONE. 2013;8:e73176.
Yue C, Li Q, Yu H. Gonad transcriptome analysis of the Pacific Oyster Crassostrea gigas identifies potential genes regulating the sex determination and differentiation process. Mar Biotechnol. 2018;20:206–19.
Broquard C, Saowaros S, Lepoittevin M, et al. Gonadal transcriptomes associated with sex phenotypes provide potential male and female candidate genes of sex determination or early differentiation in Crassostrea gigas, a sequential hermaphrodite mollusc. BMC Genomics. 2021;22:609.
Piprek RP, Damulewicz M, Tassan JP, et al. Transcriptome profiling reveals male- and female-specific gene expression pattern and novel gene candidates for the control of sex determination and gonad development in Xenopus laevis. Berlin: Springer; 2019. https://doi.org/10.1007/s00427-019-00630-y.
Zhang X, Zhou J, Li L, et al. Full-length transcriptome sequencing and comparative transcriptomic analysis to uncover genes involved in early gametogenesis in the gonads of Amur sturgeon (Acipenser schrenckii). Front Zool. 2020;17(1):11.
Sato Y, Yoshida K, Shinka T, et al. Altered expression pattern of heat shock transcription factor, Y chromosome (HSFY) may be related to altered differentiation of spermatogenic cells in testes with deteriorated spermatogenesis. Fertil Steril. 2006;86(3):612–8.
We would like to thank all members of the laboratory for valuable discussions.
This work was supported by the Key R&D Project of Shandong Province (2021ZLGX03), National Natural Science Foundation of China (32172967, 32130107), China Agriculture Research System of MOF and MARA and Taishan Scholar Project Fund of Shandong Province of China. The funding bodies had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Authors and Affiliations
MOE Key Laboratory of Marine Genetics and Breeding & Sars-Fang Centre, Ocean University of China, 5 Yushan Road, Qingdao, China
Yajuan Li, Liangjie Liu, Lijing Zhang, Huilan Wei, Shaoxuan Wu, Tian Liu, Ya Shu, Yaxin Yang, Zujing Yang, Shi Wang, Zhenmin Bao & Lingling Zhang
Laboratory for Marine Biology and Biotechnology & Laboratory for Marine Fisheries Science and Food Production Processes, Pilot National Laboratory for Marine Science and Technology (Qingdao), Qingdao, China
Shi Wang, Zhenmin Bao & Lingling Zhang
Key Laboratory of Tropical Aquatic Germplasm of Hainan Province, Sanya Oceanographic Institution, Ocean University of China, Sanya, China
LLZ and ZB developed the research questions and designed the experiment. YL, LL, and LJZ performed research. YL completed the analyses and wrote the first draft of the manuscript. HW, SXW, TL, YS, YY and ZY aided in sample collection. LLZ and SW edited the manuscript. All authors contributed critically to the drafts. All authors read and approved the final manuscript.
Animal materials of M. lateralis used in this study were obtained from our lab. No field permissions were necessary to collect the animal samples for this study. The authors declared that the experimental research on the animals described in this paper was in compliance with institutional, national and international guidelines.
Consent to 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.
Detailed information on WGCNA, differential expression analysis and gene annotation for the seven gonadal development-related modules.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Li, Y., Liu, L., Zhang, L. et al. Dynamic transcriptome analysis reveals the gene network of gonadal development from the early history life stages in dwarf surfclam Mulinia lateralis.
Biol Sex Differ13, 69 (2022). https://doi.org/10.1186/s13293-022-00479-3