Deep-sea mussel G. platifrons specimens were obtained from Site F (22°06′N, 119°17′E), a cold seep which located at a depth of 1100m in the South China Sea ( Figure 1A ). The environmental temperatures of Site F are about 4°C and other physiochemical parameters were documented by (Cao et al., 2021). The mussel materials were captured using the TITAN 4 manipulator of the ROV Faxian, operated from the R/V Kexue during September 2017, and the mussels were transported in an isothermal bio-box ( Figure 1B ) in the underframe of the ROV, which kept the water temperature around 8–9°C. The gonads of adult G. platifrons were dissected immediately after ROV retrieval on deck and then split into two pieces. Each gonad sample was fixed in 4% paraformaldehyde (PFA) overnight, washed twice with 0.01 M phosphate-buffered saline (1 × PBS) and then stored in 70% ethanol for sex identification. The remainder were flash-frozen in liquid nitrogen and stored at -80°C for RNA and DNA extraction. RNA and DNA were extracted from the testes and ovaries using Trizol (Thermo, Waltham, USA) and E.Z.N.A.^® Mollusc DNA kit (Omega, Norcross, USA), according to the manufacturer’s instructions, respectively.

The gonad samples were dehydrated with serial ethanol (80, 90, 100, and then 100%) and then transferred to xylene and embedded in paraffin wax. The 7 μm sections of gonad samples were cut on a Leica microtome and tiled on glass slides. Gonad sections were deparaffinized with xylene, hydrated with serial ethanol (100, 100, 90, 80, and then 70%), and stained with hematoxylin-eosin (HE) to determine the sex and development stages of the gonads ( Figures 1E, F , other gonad sections of males and females are shown in Supplementary Figure 1 ).

After sex identification, 12 female gonads and 12 male gonads were chosen for RNA extraction. RNA degradation and contamination were monitored on 1% agarose gels. RNA quality and integrity were checked using an Agilent 2100 Bioanalyzer, and the concentration of total RNA was quantified using a Nanodrop 1000 spectrophotometer. Qualified RNA was subsequently used for RNA-seq analysis. The RNA of each four female gonads was mixed in equal amounts to obtain three RNA pools of 12 female gonads and similar operation was performed for 12 male gonads. Herein, three female and three male gonadal RNA pools were prepared for RNA-seq. Six RNA samples were sent to Novogene Corporation (Beijing, China) for cDNA library construction. Then, library preparations were sequenced on an Illumina Hiseq 2500 platform and 125 bp paired-end reads were generated.

The WGBS sample scheme was same as RNA seq, and three female and three male gonadal DNA pools were obtained. For each DNA pool, 5.2 μg of genomic DNA spiked with 26 ng lambda DNA were fragmented by sonication to 200—300 bp with Covaris S220, followed by end repair and adenylation. These DNA fragments were treated twice with bisulfite using EZ DNA Methylation-GoldTM Kit (Zymo Research), before the resulting single-strand DNA fragments were PCR amplified using KAPA HiFi HotStart Uracil + ReadyMix (2X). Library concentration was quantified using Qubit^® 2.0 Flurometer (Life Technologies, CA, USA) and quantitative PCR, and the insert size was assayed on an Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and 125 bp paired-end reads were generated. Image analysis and base calling were performed with Illumina CASAVA pipeline, and generating 125bp paired-end reads.

Raw sequencing data were processed using fastp (version 0.20.0) with default parameters (Chen et al., 2018) to remove the adaptor sequences and low-quality reads of the RNA-seq and WGBS raw reads. Then the clean reads of RNA-seq and WGBS were prepared for downstream analysis.

The clean reads were mapped with an updated reference genome for G. platifrons (unpublished genome version, updated by Dr. Minxiao Wang) using HISAT2 (version 2.2.1). The mapped genes for G. platifrons were annotated in Nr, Swissprot, InterProscan, KEGG and the COG database. The transcripts of all the identified genes were quantified with featureCounts software (Version 2.0.1). Person correlation ( Figure 1C ) and principal component analysis (PCA) ( Figure 1D ) were then performed using custom-written R scripts to assess the relationships between gonad samples ( Figure 1B ). To identify sex-biased genes, the DESeq2 R package (version 1.10.1) was used. Only genes with | log₂ (fold change) | > 1 and a corrected p value (q value) < 0.05 were regarded as differentially expressed genes (DEGs). Gene ontology (GO) and KEGG enrichment analysis of DEGs were conducted using Blast2GO and KOBAS software, GO terms and KEGG pathways with q values less than 0.05 were considered as significantly enriched. The results of enrichment analysis were visualized using ggplot2 R package (version 3.2.1) and OmicShare tools (www.omicshare.com/tools). Genes with GO annotation of 0007530 (sex determination) and 0007548 (sex differentiation) terms were selected as sex-related gene sets.

To find the similarities and differences in the gonadal transcriptomes and sex-related genes of deep-sea and shallow water mussels, the published transcriptome data of shallow water mussel gonads were selected for comparative transcriptome analysis. Transcriptome data for female and male gonads of only one shallow water mussel L. fortunei are published. The transcriptome and genome data for L. fortunei were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/bioproject/587212) and GIGADB (https://ftp.cngb.org/pub/gigadb/pub/10.5524/100001_101000/100386/), respectively. All data were reanalyzed in line with the method mentioned previously. In short, first, the DEGs among the two sexes of L. fortunei were determined by DESeq2 and the GO and KEGG enrichment analysis were conducted using Blast2GO and KOBAS software. We then used OrthoFinder (version 2.5.2) to find the homologous protein sequences of G. platifrons and L. fortunei. Next, we used the G. platifrons sex-related genes and screened out single copy orthologous genes to obtain the orthologue sex-related genes in both G. platifrons and L. fortunei. Finally, we compared the expression pattern of shared sex-related genes in female and male gonads between G. platifrons and L. fortunei.

Bismark software (version 0.16.1) was used to perform alignments of bisulfite-treated reads to a reference genome with the default parameters. The G. platifrons reference genome was first transformed to a bisulfite-converted version (C-to-T and G-to-A converted) and then indexed using bowtie2. Sequence reads were also transformed into fully bisulfite-converted versions (C-to-T and G-to-A converted) before they are aligned to similarly converted versions of the genome in a directional manner. Sequence reads that produce a unique best alignment from the two alignment processes (original top and bottom strand) were then compared to the G. platifrons genomic sequence and the methylation state of all cytosine positions in the read was inferred. Read pairs that shared the same coordinates in the genome were regarded as duplicated, and removed before calling methylation state, thus avoiding potential methylation level calculation bias. The bisulfite non-coversion rate was calculated as the percentage of cytosines sequenced at cytosine reference positions in the lambda genome. To identify the methylation site, we modeled the sum of 5mC of methylated counts as a binomial (Bin) random variable with methylation rate (r), as 5mC ∼ Bin (5mC + umC*r). To calculate the methylation level of the sequence, we divided the sequence into multiple bins, with each bin size at 10 kb. The sum of methylated and unmethylated read counts in each window were calculated. The 5mC level (ML) for each window or C site shows the fraction of methylated Cs, and is defined as: ML(5mC) =reads (5mC)/(reads (5mC) + reads (umC)). The relative proportion of 5mCs in the three regions including 2kb upstream of the genes’ transcription start sites (TSS2K), the gene body, and 2kb downstream of transcription termination sites (TTS2K) were calculated as the proportion of mCG, mCHG, and mCHH of the total 5mC sites, respectively. Differentially methylated regions (DMRs) of the two sexes were identified using the DSS software (version 2.12.0). The regions with a p value less than 1e-05, length more than 50, average methylation level difference between groups greater than or equal to 20%, and at least three CpG sites were retained as a final DMR. We defined the genes related to DMRs as genes whose gene body region (from TSS to TTS) or TSS2K have an overlap (> 1bp) with the DMRs. Moreover, genes that were DEGs and had differential methylated promoter (DMP) regions were selected for the further analysis.

The average methylation level of all DEGs and sex-related genes in the three sequence regions (TSS2K, Gene body and TTS2K) were calculated and then compared between females and males. Conventionally, the major region of focus is promoter, so analyzed the 5mC level in the promoter of these genes between females and males.

Statistical analyses of DEGs and sex-related genes methylation level between females and males were undertaken in GraphPad (version 7.04). The p values were calculated using non-parametric tests. P < 0.05 was considered statistically significant.

To provide a global view of transcriptional differentiation and determine the molecular mechanisms underlying sex maintenance between male and female G. platifrons, six cDNA libraries were constructed and utilized for Illumina paired-end sequencing to generate representative transcripts of a wide range of biological processes. As shown in Table 1 , after eliminating the primer and adapter sequences as well as the low-quality reads by fastp, the number of clean reads ranged between 38391192–50489378 with a relatively high base quality (Q30 > 92.63% for all samples). The unique mapping rate of all samples ranged from 63.35 to 69.20% when mapped against the updated G. platifrons reference genome (unpublished data). We obtained a total of 32,918 expressed genes in female and male gonads according to featureCounts ( Supplementary Table 1 , genes with expression).

DEG analysis revealed 5923 sex-biased genes which were differentially expressed between the ovaries and testes ( Figure 1B ), 2711 of which were up-regulated in the ovaries (female-biased genes) and 3212 were up-regulated in the testes (male-biased genes) ( Supplementary Table 1 , all DEGs). To investigate the functions of DEGs, GO and KEGG analysis were implemented. For GO enrichment analysis, the results showed that 434 GO terms were enriched to 257 biological processes (BP), 96 to cellular components (CC) and 81 to molecular functions (MF) (q value < 0.05, Supplementary Table 2 ; GO_BP, GO_MF, and GO_CC). Among them, 3318 DEGs were enriched to BP, and most of the BP were enriched within single-organism processes and multicellular organismal processes. Altogether, 2995 DEGs were enriched to CC. In terms of CC, cell related terms contributed the greatest proportion. A total of 3738 DEGs were enriched within MF. For MF, protein binding and catalytic activity represented the most universal terms. 3318 DEGs assigned to BP and the top 20 enriched BP are showed in Figure 2A . The processes related to cell growth and development were extensively enriched, especially those related to cytoskeleton and cell motility. Additionally, 600, 116, and 11 DEGs were assigned to “reproduction”, “sex differentiation”, and “sex determination”, respectively, which account for 18.08, 3.5, and 0.33% of all BP terms, respectively. Furthermore, 2711 female-biased genes and 3212 male-biased genes were selected for separate GO enrichment analysis. We chose the top 10 enriched GO terms in BP, CC and MF for analyzing females and males ( Supplementary Table 2 , Top10 GO terms). Among them, nutrient metabolism related processes were the most enriched GO terms for female-biased genes while cytoskeletal correlation processes accounted for the dominant GO terms enriched for male-biased genes.

For KEGG enrichment analysis, a total of 1120 candidate genes in DEGs were noted from KEGG. Of these, 470 DEGs were significantly mapped onto 15 pathways (q value < 0.05, Supplementary Table 2 , KEGG). The KEGG enrichment results were comparable with GO results, nine were related to metabolism, including carbon, glutathione, cytochrome P450, amino acids, and glycometabolism, which indicated that metabolic process related pathways play a vital role in sex determination and differentiation ( Figure 2B ). We found that pathways of microbial metabolism in diverse environments were the most enriched pathways, the other 20 related microbial metabolism pathways are shown in Figure 2C . We also found that 17 DEGs mapped onto the WNT signaling pathway ( Supplementary Table 2 , WNT), including Wnt6 (CSbapl.23381) and Wnt7 (CSbapl.46743) which function as signaling molecules in sexually dimorphic development. Wnt6 and Wnt7 were up-regulated in female G. platifrons gonads and the canonical WNT signaling mediator, β-catenin, was also detected in transcriptome and up-regulated in female gonads, which indicated the conservatism of sex determination in the deep-sea mussel G. platifrons.

Of the total candidate genes, 161 were identified as sex-related genes generated using the “sex determination” and “sex differentiation” terms ( Supplementary Table 2 , sex-related genes). These genes could be categorized into transcriptional regulatory genes including the transcription factor, steroid hormone biosynthesis related genes, ovary and testis specific genes, gametogenesis related genes, some kinase genes related to cell differentiation, and development and other genes involved in sex determination or differentiation were mentioned in the literature. The expression patterns of 161 sex-related genes are shown in Figure 2D and Supplementary Table 2 . Among these candidate genes, six vitelline envelope zona pellucida genes (VEZP type 6/9/14/17/24), one vitelline membrane outer layer protein 1 (VMO1), and one vitellogenin (VTG) were classified as ovary or oocyte-specific genes which were expressed more in females than males. Twenty-one genes had testis-specific expressions and 17 were sperm-specific genes such as testis-specific serine kinase (TSSK), testis-expressed protein (TEX), spermatogenesis-associated protein (SPATA), and sperm surface protein (SP); all of these genes were more highly expressed in males. Additionally, four genes annotated as the SOX gene family were found, including SOX2 in the SOXB1 group, SOX4 in the SOX C group, SOX5 in the SOX D group, and SOX9 in the SOX E group. The main transcription factors such as FOXL2 and DMRT2 involved in female and male sex determination, respectively were also identified, and the FOXL2 were expressed more in females, while DMRT2 were high expression in males. Moreover, sex steroidogenesis-related genes such as estradiol 17β-dehydrogenase (E2DH), 17β-hydroxysteroid dehydrogenase (17β-HSD type 6/10/12/14) and CYP17a were all significantly high expressed in females. The roles of these genes need further study in G. platifrons, but most of them have been identified as playing important roles in sex-determination/differentiation or gametogenesis in other animals and molluscan species.

To find the similarities and differences between deep-sea mussels and shallow water relatives, we first used OrthoFinder (version 2.5.2) to screen the homologous genes of deep-sea mussel G. platifrons and shallow water mussel L. fortunei. A total of 2543 orthogroups were documented ( Supplementary Table 3 , othogroups), uncovering 35 shared sex-related genes by intersection with 161 sex-related genes in G. platifrons. Then, we compared the shared expression patterns of 35 sex-related genes between a deep-sea mussel and shallow water mussel. Thirty-five shared sex-related genes contained most of the genes that play an important role in sex determination and differentiation, such as SOX gene family, FOXL2, β-catenin and some testis-specific genes. For female high expressed genes, such as SOX2, FOXL2 and β-catenin, SOX2 high expressed in female than male about 10 times in both two species, FOXL2 were about 100 times more expressed in females than males of G. platifrons, and about 900 times more expressed in females than males of L. fortunei. The same expression patterns also found in β-catenin which were all more than 2 times expressed in females than males in both species. As similar expression pattern was observed in male-upregulated genes in both shallow water and deep-sea mussels, such as TEX36 and TEX47, which about 2-10 times more expressed in males than females ( Figure 3 and Supplementary Table 3 , shared sex-related genes).

To examine the roles of DNA methylation in sex maintenance and gonad development in adult G. platifrons, the methylation patterns and levels of gonads for female and male G platifrons were detected with WGBS analysis ( Supplementary File 1 ). We found that the 5mC levels of all genomic C sites for males (3.22%) were slightly higher than for females (3.09%) and the 5mC levels for the CG site (>96%) were higher than for the CHG (<1%) and CHH (~2%) sites. For the GC site, the average 5mC levels in the gene body (30-40%) were higher than that in other sequence regions such as promoter (~20%) and repeat regions (~10%) in the genome. We also found that the 5mC level showed a weak positive correlation with gene expression level (Spearman correlation r range from 0.313 to 0.543) at the three sequence regions (TSS2K, Gene body, TTS2K) and in both sexes ( Supplementary Table 4 , F_TTS2K_5mC_vs_gene_express, M_TTS2K_5mC_vs_gene_express, F_genebody_5mC_vs_gene_express, M_genebody_5mC_vs_gene_express, F_TSSK2k_5mC_vs_gene_express, and M_TSSK2k_5mC_vs_gene_express).

The differential 5mC level of DEGs has the same trend with all genes being detected in the three sequence regions, and the methylation levels were all significantly higher in males than females ( Figure 4A and Supplementary Table 4 , 5mC_DEGs). To identify the role of methylation more clearly in sex determination or sex differentiation, we screened the 1903 shared genes of DMRs and DEGs for further analysis ( Figure 4B , Supplementary Table 4 , DMRs_vs_DEGs). GO functional enrichment analysis showed that 1180 shared genes were enriched in the BP category, and 256 BP GO terms were significantly enriched (q value < 0.05). The top 20 enriched BP terms are shown in Figure 4C . Additionally, BP terms related to “sex differentiation” and “male sex differentiation” were all significantly enriched ( Supplementary Table 4 , BP_shared_genes_ DMRs_vs_DEGs). For KEGG analysis, 419 shared genes between DMRs and DEGs were associated with 335 pathways. The top 20 enriched pathways are shown in Figure 4D . A total of 128 shared genes were significantly enriched in metabolism and hormone biosynthesis related pathways (q value < 0.05). Some of the shared genes were enriched in the “WNT signaling pathway”, “oocyte meiosis”, and “FoxO signaling pathway”, which play important roles in sex determination or differentiation and germ cell development.

The 5mC levels of each sex-related gene and DEG were documented in Supplementary Table 4 . In total, 59 of 161 (36.65%) sex-related genes in G. platifrons displayed significantly differential methylation patterns between females and males ( Supplementary Table 4 , DMR_vs_sex-related genes), which contrasts to ~ 20.97% (6904/32,918) over whole genome (Fisher’s exact test, p < 0.0001). To find out if the 5mC levels of sex-related genes were different from other DEGs, the 5mC levels of the TSS2K, the gene body and the TTS2K regions of 161 sex-related genes were compared with the other DEGs (5762 genes) using Kolmogorov-Smirnov tests. We defined the OF and OM as the other DEGs in females and other DEGs in males. The result shows that for the TSS2K region, the 5mC levels of sex-related genes for females (F) were significantly higher than for the other DEGs in females (OF) ( Figure 5A , F vs OF, p=0.0169) while no significantly differences were found in males ( Figure 5A , M vs OM, p=0.0977). For the gene body region, the 5mC levels of sex-related genes were not significantly different from OF ( Figure 5B , F vs OF, p=0.0582), while 5mC levels of sex-related genes were significantly higher than OM ( Figure 5B , M vs OM, p=0.0273). For the TTS2K region, 5mC levels of sex-related genes were not significantly different from OF ( Figure 5C , F vs OF, p=0.2180) and OM ( Figure 5C , M vs OMp=0.1567), respectively. Moreover, the 5mC levels of 161 sex-related genes for males were all significantly higher than for females in the gene body and TTS2K regions ( Figure 5D , Gene body and TTS2K, p < 0.0001). However, for the TSS2K region, there were no significant differences between females and males ( Figure 5D , TSS2K, p = 0.2113).

We then divided 161 sex-related genes into 73 female-biased sex-related genes (FBSGs) and 88 male-biased sex-related genes (MBSGs) according to expression levels. We found that 5mC levels were all significantly higher in MBSGs than FBSGs for the two sexes and the three regions ( Figure 6 , p < 0.0001). To find if the 5mC levels of FBSGs and MBSGs were significantly different in females and males, we conducted Wilcoxon matched-pairs signed rank tests for the three regions ( Figure 6 ). The results showed that, except for the levels of FBSG 5mC levels in the gene body regions, there were no significant differences between the two sexes ( Figure 6 , p = 0.0554). Other 5mC levels of FBSGs and MBSGs were all higher in males than females in all three regions (p < 0.0001).

We subsequently investigated whether the sex-related genes overlapped with the DMPs ( Supplementary Table 4 , DMP_vs_sex-related genes). Only 12 sex-related genes were found by transcriptome and methylation associated analysis. These genes include male-biased genes such as testis-specific serine kinase (CSbapl.702), a WNT pathway gene β-catenin (CSbapl.3788) which is believed to be essential for female sex determination, and the transcription factor SOX family gene SOX5. Except for gene Prx, all 5mC levels of females were lower than those of males, which is in line with the trend of 5mC levels of all sex-related genes. For five female-biased genes, the hypomethylation of promoter regions corresponded with high methylation levels in promoter regions, which would suppress gene expression. However, for six male-biased genes except Prx, the hypermethylation of promoter regions resulted in high gene expression.

Many genes related to sex determination are transcription factors, especially in terms of the DMRT, FOX and SOX, gene families (Smith et al., 2009; Matson et al., 2011; Huang et al., 2017; Li et al., 2018). By analyzing the female and male gonad transcriptome of G. platifrons, we identified at least 161 sex-related genes which were all belonged to DEGs and annotated to sex determination or sex differentiation GO term. Sixteen transcription factors were detected, including the common genes DMRT2, FOXL2, and SOX9, which were shown to be deeply conserved in sex determination or differentiation in animals (Liu et al., 2012; Li et al., 2018; Afonso et al., 2019).

The DMRT family of genes contains a highly conserved DNA binding domain called the DM domain. There are seven DMRT genes that have been shown to play critical roles in sexual regulation, including sex differentiation, sexual dimorphism, and spermatogenesis in a broad variety of metazoans, including nematodes, insects, and vertebrates (Zarkower, 2013). Another famous sex determination transcription factor is FOXL2. FOXL2 is one of the FOX gene family members and preferentially expressed in the ovary of vertebrates, controlling ovarian differentiation and maintenance by repression of testis-specific genes from the early embryonic gonad throughout adult life (Uhlenhaut et al., 2009). Therefore, FOXL2 and DMRT are the key genes for female and male sex determination and maintenance, respectively. In addition, they are mutually antagonistic in determining the sex of mammals (Matson et al., 2011). In G. platifrons, the DMRT2 gene was examined and found to be up-regulated in male gonads (q < 0.05), and FOXL2 was at higher levels in females than males (q < 0.05). As in mammals, the balance between FOXL2 and DMRT2 is suggested to be important for sex determination in G. platifrons. This phenomenon was also reported in other mollusca species such as the blacklip pearl oyster Pinctada margaritifera, Pacific oyster Crassostrea gigas, Yesso scallop Patinopecten yessoensis and razor clam Sinonovacula constricta, although the genes in the DMRT family and FOX family are not always the same among different species, for example DMRT2 and FOXL2 in P. margaritifera, Dsx and FOXL2 in C. gigas, DMRT1L and FOXL2 in P. yessoensis, FOXA1, FOXD2, and DMRTA2 in S. constricta (Teaniniuraitemoana et al., 2014; Zhang et al., 2014; Li et al., 2018; Yao et al., 2021).

The male-determining gene in mammals is the SRY gene. SRY is a member of the SOX gene family of DNA-binding proteins. SRY complexed with the SF1 protein would upregulate other transcription factors including the important SOX9 which is the up-stream regulator of male sex determination (Sinclair et al., 1990; Kashimada and Koopman, 2010). SOX9 is one gene in the SOXE group that has an expression in the Sertoli cells necessary for testis determination, and is reinforced by FGF9 and DMRT1 which are thus essential to this network (Matson and Zarkower, 2012). In mammals, once SRY expression begins, other genes start to express only in one sex, so the SOX9 would not be active in female gonads (Gross, 2006). However, in our study, although no SRY was detected, the downstream gene SOX9 was found in the gonad transcriptome, and SOX9 was not only active in male testes but also in female ovaries. Moreover, SOX9 was not only highly expressed in males but also female. High expression in females may explain the phenomenon of hermaphroditism previously discovered in G. platifrons (Zhong et al., 2020). For other gonadal transcriptomes of shallow water bivalves, SOX9 expression levels were also determined, such as in P. yessoensis, C. gigas and L. fortunei, but for these bivalves the SOX9 expression levels were not significantly different between males and females (Zhang et al., 2014; Li et al., 2016; Afonso et al., 2019). Overall, we suggest that SOX9 could be a key factor in sex determination or differentiation in G. platifrons, but it may not function in the way that it does in mammals or indeed other shallow bivalves.

We also identified the β-catenin gene in G. platifrons. In mammals, β-catenin, the center in the canonical WNT signaling pathway, is a key transcriptional regulator and is essential for ovarian development and maintenance (Liu et al., 2009). RSPO1-WNT4-β-catenin in the canonical WNT signaling pathway is independent and complementary to the FOXL2-leading pathway (Kocer et al., 2008). β-catenin presents in the gonads of both sexes of G. platifrons but was expressed significantly higher in ovaries, while RSPO1 was not found, and Wnt4 was not a DEG in the transcriptome of G. platifrons. Therefore, the existence and specific function of β-catenin for female sex determination needs further investigation.

In non-mammalian vertebrates, sterol hormones play an important role in sex determination and sex differentiation (Navarro-Martin et al., 2011). Aromatase CYP19a (Cytochrome P450 family 19) can irreversibly converts androgens into estrogens, which influences the ratio of androgen-to-estrogen, and the different ratio then determines undifferentiated gonad sexually differentiates into a testis or ovary. Sixteen Cytochrome P450 family genes were found in the deep-sea mussel gonadal transcriptome and eight were found in females and males’ DEGs, however, CYP19a were not found in these data, but another important Cytochrome P450 family 17, one of the sex-related genes CYP17a was found both in two sexes of G. platifrons which were high expression in female. CYP17a has been shown to play important roles in female sex determination and male fertility by regulating sex steroid biosynthesis in fish (Yang et al., 2021). For G. platifrons, 9.88% were methylated in the testes and 8.21% were methylated in the ovaries for the TSS2K region of CYP17a, however, CYP17a was not defined as a DMR. So, the effect of methylation on CYP17a is not as pronounced as it is on CYP19a which hypermethylation occurs in testis and gene high expression in ovaries of fish (Wen et al., 2014), and it also possible that the mechanism of sterol hormones in mollusks is obviously different from that in vertebrates.

Specifically, we compared the sex-related genes between the deep-sea mussel G. platifrons and the shallow water mussel L. fortunei. The results showed that the common expression pattern of the male-biased sex-related genes and the expression level of those genes were all coincident among male samples, while for female-biased sex-related genes, the expression pattern and gene expression levels were inconsistent in L. fortunei. We suggest that the common expression pattern for female G. platifrons may be caused by the stable deep sea environment or the same development stages of gonads. Therefore, there are two factors that may lead to differences in the expression patterns of female-biased sex-related genes between G. platifrons and L. fortunei. First, the development cycle is different at the time of sampling for female gonads of L. fortunei. Second, the shallow water environment is not stable, and the maternal immunity may be influential, given that the SOX family genes are important genes that cause differential expression patterns in female individuals of L. fortunei, and the SOX family genes are reported to be related to immunity (Furukawa et al., 2021). The ovary is an important source of maternal immunity, therefore the differential expression between individuals may be caused by their different responses to the changing environment (Wang et al., 2015). In addition, the other important sex determination genes like FOXL2, β-catenin, and Wnt7 all have same the expression patterns between G. platifrons and L. fortunei ( Figure 3 ), suggesting some conserved sex determination mechanisms between these two mussels.Based on previous studies and our transcriptome data, all the aforementioned genes which are shared among G. platifrons and other animals are thought to be involved in sex determination, and thus we speculate that the sex determination or differentiation mechanism may be conserved in some degree among deep-sea mussels and shallow water species or other vertebrates, but there are also significant differences in the male and female expression patterns of specific sex-related genes. As for G. platifrons, based on sequence homology and functions inferred from transcriptome data, we suggest that although traditionally identified sex-related genes are not strictly expressed in male and female gonads, the expression of these genes or homologues, especially SOX, FOXL2, and DMRT genes, which are commonly referred to in other mollusks, are significantly different between males and females. We suggest that SOX may have high expression in the early developmental stage before G. platifrons gonadal maturation and in the up-stream regulator leading to male sex determination. However, in the gonadal maturation stage, the SOX gene would not highly express in testes. Thus, SOX would directly or indirectly activate DMRT2 for testis development in early stages and DMRT2 may keep testis functional through a life. FOXL2 may be highly expressed in ovaries and complemented by β-catenin leading to female sex-determination pathways in G. platifrons.

In the present study, we have shown that G. platifrons has some sex determination or differentiation mechanisms common with the shallow water mussel L. fortunei and other bivalves. Additionally, G. platifrons is shown to be a gonochoric species (Zhong et al., 2020), while the sex ratio is not coincident among groups of different shell lengths. Additionally, while hermaphroditism was observed in our samples, this was only the case in two large individuals, so we suspect that this species may not exhibit large-scale sex-reversal and sex maintenance. We recorded higher levels of global DNA methylation level in the testis (3.22%) than ovary (3.09%), suggesting that male germ cells may be more heavily methylated than female germ cells, which concurs with findings reported for P. yessoensis (Li et al., 2019). Additionally, we also found that 5mC is predominantly in the gene body (30-40%) in both sexes, which is in agreement with previous studies of C. gigas (Gavery and Roberts, 2013; Olson and Roberts, 2014; Wang et al., 2014). Through enrichment analysis of the intersection of DMRs and DEGs in two sexes of G. platifrons, we found that methylation mainly affects metabolic processes and morphological differentiation of organisms ( Figure 4C, D ). Therefore, we hypothesize that methylation plays an important role in the differential metabolism of males and females in deep-sea mussels, in particular high 5mC levels may have a vital effect on male sex determination.

To further explore the effect of methylation on sex determination or differentiation, we chose to study mainly the sex-related genes. The methylation levels presented negative correlations with expression levels of sex-related genes in the ovaries of G. platifrons, which corroborates the concept that increased methylation corresponds to decreased gene expression (Riviere et al., 2013). In males, our results indicated the positive correlation between 5mC and sex-related gene expression, where genes with high expression levels had high 5mC levels ( Supplementary Figure 2 ). This is similar to the observations reported for oyster male gamete (Gavery and Roberts, 2013; Olson and Roberts, 2014). We also found that the tendency for 5mC levels of all sex-related genes and MBSGs or FBSGs were all higher in males than females in the three sequence regions, although some comparisons were not statistically significantly different ( Figures 5D and 6 ). Moreover, similar patterns were found at the level of the single sex-related gene ( Supplementary Table 4 , 5mC_sex-related genes) indicating that females present more hypomethylation than males, but in males hypermethylation correlated with high levels of gene expression in all sequence regions ( Supplementary Figure 2 ). We hypothesize that DNA methylation may act as an enhancer rather than a repressor in testis development. As show in Figure 7 , we integrated the differential methylation of several key sex-related genes to speculate the sex-determining mechanism of deep-sea mussel G. platifrons.

For the result stated above, we suggest that the different methylation levels of sex-related genes in the two sexes and the inconsistent correlation between methylation level and sex-related genes expression level may play an important role for sex determination or differentiation of G. platifrons in the deep-sea.