Dynamic transcriptome analysis reveals the gene network of gonadal development from the early history life stages in dwarf surfclam Mulinia lateralis
Biology of Sex Differences volume 13, Article number: 69 (2022)
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.
Transcripts per million
False discovery rate
Differentially expressed genes
Germ stem cells
Principal component analysis
- GATA2 :
GATA-binding factor 2
- OVO :
Transcriptional regulator ovo
- NR2F1 :
Nuclear receptor subfamily 2 group F member 1
- NR1D3 :
Nuclear receptor subfamily 1 group D member 3
- ER :
- WNT7B :
Protein Wnt-7b precursor
- HTR4 :
5-Hydroxytryptamine receptor 4
- 5-HTR1 :
5-Hydroxytryptamine receptor 1
- FZD10 :
- SSTR2 :
Somatostatin receptor type 2
- NRXN3 :
- NLGN4Y :
- Rut :
Ca(2+)/calmodulin-responsive adenylate cyclase
- GNRHR :
Gonadotropin-releasing hormone receptor
- GATA4 :
GATA-binding factor 4
- SIX4 :
Homeobox protein SIX4
- HES1 :
Transcription factor HES-1
- GAL3ST3 :
- FOXL2 :
Forkhead box protein L2
- FASN :
Fatty acid synthase
- RPN1 :
26S proteasome regulatory subunit RPN1
- DDOST :
Dolichyl-diphosphooligosaccharide–protein glycosyltransferase 48 kDa subunit
- VIT6 :
- MP62 :
Mitotic apparatus protein p62
- CDC20 :
Cell division cycle protein 20
- FOXZ :
Forkhead box protein Z
- HSFY :
Heat shock transcription factor, Y-linked
- SOXH :
Transcription factor SOXH
- PRDM9 :
Histone-lysine N-methyltransferase PRDM9
- SYCP3 :
Synaptonemal complex protein 3
- HFM1 :
Probable ATP-dependent DNA helicase HFM1
- MSH4 :
MutS protein homolog 4
Testis-specific serine/threonine-protein kinase family genes ¼
- SMC5 :
Structural maintenance of chromosomes protein 5
- MRE11 :
Double-strand break repair protein MRE11
- CDK1 :
Cyclin-dependent kinase 1
- RFC1 :
Replication factor C subunit 1
- 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
- PABP1 :
Polyadenylate-binding protein 1
Hu Y, Okumura L, Page D, et al. Gata4 is required for formation of the genital ridge in mice. PLoS Genet. 2013;9(7):e1003629.
Miyamoto N, Yoshida M, Kuratani S, et al. Defects of urogenital development in mice lacking Emx2. Development. 1997;124(9):1653–64.
Birk O, Casiano D, Wassif C, et al. The LIM homeobox gene Lhx9 is essential for mouse gonad formation. Nature. 2000;403(6772):909.
Eggers S, Ohnesorg T, Sinclair A. Genetic regulation of mammalian gonad development. Nat Rev Endocrinol. 2014;10(11):673–83.
Campolo F, Gori M, Favaro R, et al. Essential role of Sox2 for the establishment and maintenance of the germ cell line. Stem Cells. 2013;31:1408–21.
Boyer L, Mathur D, Jaenisch R. Molecular control of pluripotency. Curr Opin Genet Dev. 2006;16(5):455–62.
Ying Y, Zhao G. Cooperation of endoderm-derived BMP2 and extraembryonic ectoderm-derived BMP4 in primordial germ cell generation in the mouse. Dev Biol. 2001;232(2):484–92.
Pesce M, Scholer H. Oct-4: control of totipotency and germline determination. Mol Reprod Dev. 2000;55:452–7.
Irie N, Weinberger L, Tang W, et al. SOX17 is a critical specifier of human primordial germ cell fate. Cell. 2015;160:253–68.
Lin I, Chiu F, Yeang C, et al. Suppression of the SOX2 neural effector gene by PRDM1 promotes human germ cell fate in embryonic stem cells. Stem Cell Rep. 2014;2(2):189–204.
Okuzaki Y, Kaneoka H, Suzuki T, et al. PRDM14 and BLIMP1 control the development of chicken primordial germ cells. Dev Biol. 2019;455(1):32–41.
Sinclair A, Berta P, Palmer M, et al. A gene from the human sex-determining region encodes a protein with homology to a conserved DNA-binding motif. Nature. 1990;346:240–4.
Polanco J, Koopman P. Sry and the hesitant beginnings of male development. Dev Biol. 2007;302(1):13–24.
Boulanger L, Pannetier M, Gall L, et al. FOXL2 is a female sex-determining gene in the goat. Curr Biol. 2014;24(4):404–8.
Kazuma T, Kaori H, Rina Kitada Y, et al. R-spondin1 plays an essential role in ovarian development through positively regulating Wnt-4 signaling. Hum Mol Genet. 2008;17(9):1278–91.
Uhlenhaut NH, Jakob S, Anlag K, et al. Somatic sex reprogramming of adult ovaries to testes by FOXL2 ablation. Cell. 2009;139(6):1130–42.
Li L, Dong J, Yan L, et al. Single-cell RNA-seq analysis maps development of human germline cells and gonadal niche interactions. Cell Stem Cell. 2017;20:858–73.
Wanninger A, Wollesen T. Mollusca (chapter 7). In: Wanninger A, editor. Evolutionary developmental biology of invertebrates 2: Lophotrochozoa (Spiralia). 1st ed. Vienna: Springer; 2015.
FAO. The state of world fisheries and aquaculture 2020. In: Sustainability in action. Rome. 2020. https://doi.org/10.4060/ca9229en.
Matoo B, Neiman M. Bringing disciplines and people together to characterize the plastic and genetic responses of molluscs to environmental change. Integr Comp Biol. 2021;61(5):1689–98.
Breton S, Capt C, Guerra D, et al. Sex-determining mechanisms in bivalves. In: Leonard J, editor., et al., Transitions between sexual systems. Cham: Springer; 2018.
Rachel C. Phylogenetic patterns and phenotypic plasticity of molluscan sexual systems. Integr Comp Biol. 2013;4:723–35.
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.
Xu R, Li Q, Yu H. Expression pattern of Piwi-like gene implies the potential role in germline development in the Pacific oyster Crossosrea gigas. Aquac Rep. 2020;18:100486.
Liu L, Liu T, Wu S, et al. Discovery of Nanos1 and Nanos2/3 as germ cell markers during scallop gonadal development. Mar Biotechnol. 2022;24:408–16.
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.
Yu R, Li Q. Seed production technology for marine molluscs in China. Qingdao: Ocean University of China Press; 2016.
Calabrese A. Mulinia lateralis: Molluscan fruit fly?. In: Proceedings of the national shellfisheries association. 1969.
Lu J, Chen T, Allen S, et al. Production of transgenic dwarf surfclams, Mulinia lateralis, with pantropic retroviral vectors. Proc Natl Acad Sci. 1996;93:3482–6.
Yang Z, Huang X, Wang H, et al. Effects of microalgae diets and stocking density on larval growth, survival and metamorphosis of dwarf surfclam, Mulinia lateralis. Aquaculture. 2021;536:736440.
Dobin A, Davis C, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9(1):559.
Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4(1):17.
Beißbarth T, Speed TP. GOstat: find statistically overrepresented gene ontologies within a group of genes. Bioinformatics. 2004;20(9):1464–5.
Oldham M, Konopka G, Iwamoto K, et al. Functional organization of the transcriptome in human brain. Nat Neurosci. 2008;11:1271–82.
Ponomarev I, Wang S, Zhang L, et al. Gene coexpression networks in human brain identify epigenetic modifications in alcohol dependence. J Neurosci. 2012;32(5):1884–97.
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.
Li R, Zhang L, Li W, et al. FOXL2 and DMRT1L are yin and yang genes for determining timing of sex differentiation in the Bivalve Mollusk Patinopecten yessoensis. Front Physiol. 2018;9:1166.
Wei H, Li W, Liu T, et al. Sexual development of the hermaphroditic scallop Argopecten irradians revealed by morphological, endocrine and molecular analysis. Front Cell Dev Biol. 2021;9:646754.
Zhang Y, Qin Y, Ma L, et al. Gametogenesis from the early history life stages of the Kumamoto Oyster Crassostrea sikamea and their breeding potential evaluation. Front Physiol. 2019;10:524.
Guo X, Allen S. Sex determination and polyploid gigantism in the Dwarf Surfclam (Mulinia Lateralis Say). Genetics. 1994;138(4):1199–206.
Sabik OL, Calabrese GM, Taleghani E, et al. Identification of a core module for bone mineral density through the integration of a co-expression network and GWAS data. Cell Rep. 2020;32(11):108145.
Gareth MR, Stephanie S, Moises FA, et al. Sex differences in developmental patterns of neocortical astroglia: a mouse translatome database. Cell Rep. 2022;38(5):110310.
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.
Zhang L, Rui H, Su H, et al. Network analysis of oyster transcriptome revealed a cascade of cellular responses during recovery after heat shock. PLoS ONE. 2012;7(4):e35484.
Zhou L, Liu Z, Dong Y, et al. Transcriptomics analysis revealing candidate genes and networks for sex differentiation of Yesso scallop (Patinopecten yessoensis). BMC Genomics. 2019;20:671.
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.
Zhao W, Li Y, Wang Y, et al. Gonadal development and the reproductive cycle of the geoduck clam (Panopea japonica) in the Sea of Japan. Aquaculture. 2022;548:737606.
De Sousa JT, Milan M, Bargelloni L, et al. A microarray-based analysis of gametogenesis in two Portuguese populations of the European clam Ruditapes decussatus. PLoS ONE. 2014;9(3):e92202.
Galindo-Torres P, Garcia-Gasca A, Llera-Herrera R, et al. Sex determination and differentiation genes in a functional hermaphrodite scallop, Nodipecten subnodosus. Mar Genomics. 2018;37:161–75.
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.
Golestaneh N, Beauchamp E, Fallen S, et al. Wnt signaling promotes proliferation and stemness regulation of spermatogonial stem/progenitor cells. Reproduction. 2009;138(1):151–62.
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.
Qin M, Zhang Z, Song W, et al. Roles of Figla/figla in juvenile ovary development and follicle formation during zebrafish gonadogenesis. Endocrinology. 2018;159(11):3699–722.
Kopp A. Dmrt genes in the development and evolution of sexual dimorphism. Trends Genet. 2012;28(4):175–84.
Zhang N, Xu F, Guo X. Genomic analysis of the Pacific Oyster (Crassostrea gigas) reveals possible conservation of vertebrate sex determination in a mollusc. G3 Genes Genomes Genet. 2014;4(11):2207–17.
Teaniniuraitemoana V, Huvet A, Levy P, et al. Molecular signatures discriminating the male and the female sexual pathways in the Pearl Oyster Pinctada margaritifera. PLoS ONE. 2015;10(3):e0122819.
Evensen KG, Robinson WE, Krick K, et al. Comparative phylotranscriptomics reveals putative sex differentiating genes across eight diverse bivalve species. Comp Biochem Physiol D: Genomics Proteomics. 2022;41:100952.
Liang S, Liu D, Li X, et al. SOX2 participates in spermatogenesis of Zhikong scallop Chlamys farreri. Sci Rep. 2019;9(1):76.
Wu S, Zhang Y, Li Y, et al. Identification and expression profiles of Fox transcription factors in the Yesso scallop (Patinopecten yessoensis). Gene. 2020;733: 144387.
Shu C, Wang L, Zou C, et al. Function of Foxl2 and Dmrt1 proteins during gonadal differentiation in the olive flounder Paralichthys olivaceus. Int J Biol Macromol. 2022;215:141–54.
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.
Feng C, Spiller C, Merriner D, et al. SOX30 is required for male fertility in mice. Sci Rep. 2017;7:17619.
Han F, Dong Y, Liu W, et al. Epigenetic regulation of Sox30 is associated with testis development in mice. PLoS ONE. 2014;9(5):e97203.
Stévant I, Nef S. Genetic control of gonadal sex determination and development. Trends Genet. 2019;35(5):346–58.
She ZY, Yang WX. Sry and SoxE genes: How they participate in mammalian sex determination and gonadal development? Semin Cell Dev Biol. 2017;63:13–22.
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.
Hamilton C, Revay T, Domander R, et al. A large expansion of the HSFY gene family in cattle shows dispersion across Yq and testis-specific expression. PLoS ONE. 2011;6:e17790.
Skinner B, Lachani K, Sargent C, et al. Expansion of the HSFY gene family in pig lineages. BMC Genomics. 2015;16(1):442.
Kinoshita K, Shinka T, Sato Y, et al. Expression analysis of a mouse orthologue of HSFY, a candidate for the azoospermic factor on the human Y chromosome. J Med Invest. 2006;53:117–22.
Andrews J, Garcia-Estefania D, Delon I, et al. OVO transcription factors function antagonistically in the Drosophila female germline. Development. 2000;127(4):881–92.
Bi H, Xu X, Li X, et al. CRISPR disruption of BmOvo resulted in the failure of emergence and affected the wing and gonad development in the silkworm Bombyx mori. Insects. 2019;10(8):254.
Manosalva I, González A, Kageyama R. Hes1 in the somatic cells of the murine ovary is necessary for oocyte survival and maturation. Dev Biol. 2013;375(2):140–51.
Siggers P, Smith L, Greenfield A. Sexually dimorphic expression of Gata-2 during mouse gonad development. Mech Dev. 2002;111(1–2):159–62.
Wang J, Tian G, Zheng Z, et al. Comprehensive transcriptomic analysis of mouse gonadal development involving sexual differentiation, meiosis and gametogenesis. Biol Proced Online. 2019;21(1):20.
Tao W, Chen J, Tan D, et al. Transcriptome display during tilapia sex determination and differentiation as revealed by RNA-Seq analysis. BMC Genomics. 2018;19:363.
Nuria S, Laia R, Francesc P. Improved biomarker discovery through a plot twist in transcriptomic data analysis. BMC Biol. 2022;20(1):208.
Robledo D, Ribas L, Cal R, et al. Gene expression analysis at the onset of sex differentiation in turbot (Scophthalmus maximus). BMC Genomics. 2015;16(1):973.
Matson CK, Zarkower D. Sex and the singular DM domain: insights into sexual regulation, evolution and plasticity. Nat Rev Genet. 2012;13:163.
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.
Ethics approval and consent to participate
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.
About this article
Cite this article
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 Differ 13, 69 (2022). https://doi.org/10.1186/s13293-022-00479-3