Whole transcriptome analysis of brain, gonad, liver, and muscle were conducted to explore the different growth mechanisms involved in the female, male, and pseudomale C. semilaevis. A total of 3.43 × 10⁹ raw reads were obtained from 36 libraries, including FB1-3, MB1-3, PMB1-3, FG1-3, MG1-3, PMG1-3, FL1-3, ML1-3, PML1-3, FM1-3, MM1-3, and PMM1-3 (Supplementary Table S1). After data filtering by removing adapter and low-quality reads, 3.42 × 10⁹ cleaned reads were retrieved (Supplementary Figure S1). The subsequent Pearson correlation coefficients and PCA analysis (Supplementary Figure S2) revealed high correlation coefficients (>0.9222) in the three biological replicates of the gonad, liver, and muscle. The PMB1 sample had a lower correlation (0.6210–0.7007) with other PMB samples and was removed from subsequent analysis.

After transcripts reconstruction by genome mapping and FPKM calculations, differentially expressed genes (DEGs) were identified (|log2FC| > 1 and q < 0.05) and screened from the M-vs-F, M-vs-PM, and PM-vs-F groups (Figure 1A). The DEGs included: 1) 430, 296, and 116 genes from the brain; 2) 11,701, 418, and 11,720 genes from the gonad; 3) 426, 274, and 108 genes from the liver and 4) 2,369, 289, and 1,749 genes from the muscle. When compared with the female group, the overlapped genes in the male and pseudomale groups from the brain, gonad, liver, and muscle were 23, 10,766, 43, and 1,439, respectively (Figure 1B).

The KEGG enrichment analysis (Figure 1C) for DEGs demonstrated that the cell cycle was significantly enriched in the gonad and muscle (q < 0.05). In particular, 104 (71.72%) and 107 (73.79%) DEGs of the cell cycle (145 genes) were classified in the male-vs-female and pseudomale-vs-female gonads, respectively (Supplementary Table S2). Similarly, in the muscle, 53 (36.55%) and 49 (33.79%) DEGs of the cell cycle were screened in the male-female and pseudomale-female groups, respectively. Other cell growth and death-related pathways, including oocyte meiosis, p53 signaling pathway, and cellular senescence, were also significantly enriched in the different sexual muscle tissues. Several replication and repair pathways, DNA replication, mismatch repair, and homologous recombination were also significantly enriched. In the brain and liver, cell growth and death-related pathways, including oocyte meiosis and cell cycle, were identified from the male-vs-female group (p < 0.05). However, only fewer than ten DEGs were classified. The GO enrichment analysis of DEGs was shown in Supplementary Figures S3–S6.

A total of 13,239 DEGs were obtained by overlapping the DEGs across the four tissues between different sexes to evaluate their correlation with the female-biased SSD in C. semilaevis. These DEGs were subsequently submitted for WGCNA analysis, and the gene cluster dendrogram was constructed based on the correlation coefficients of each gene (Figure 2A). Thirteen modules were obtained with module sizes of 95–4,748. The subsequent calculation of module correlation coefficient and sample growth trait (shown in Supplementary Table S3) identified two modules as the most positive and most negative growth trait-related modules (turquoise and brown, respectively in Figure 2B). The high expression levels of genes in the turquoise and brown modules were observed in the female gonad, male and pseudomale gonad, respectively (Supplementary Figure S7).

Genes in the turquoise module were significantly expressed in the following GO terms: catalytic activity, nucleic acid binding, transferase activity, metabolic processes, chromosome organization, RNA processing, and cell cycle (p < 0.05) (Supplementary Figure S8). In the brown module, GO terms for cytoplasmic dynein complex, cytoskeleton, dynein complex, and spindle microtubule were significantly enriched (q < 0.05) (Supplementary Figure S9).

KEGG enrichment was subsequently employed to reveal the possible functions of these two modules. In the turquoise module, 31 pathways, including cell cycle, spliceosome, DNA replication, etc., were significantly enriched (q < 0.001, Figure 3A). Additionally, 76 DEGs in the cell cycle (Figure 4A) and 65 DEGs in the spliceosome exhibited a female-biased mRNA expression in the gonad. In the brown module, cellular senescence, signal transduction related pathways-Hippo signaling pathway, TGF-beta signaling pathway, and Jak-STAT signaling pathway were significantly enriched (p < 0.05) (Figure 3B). A total of 99 DEGs from cellular senescence and hippo signaling pathway displayed male- and pseudomale-biased mRNA expression in the gonad (Figure 4B). The relationship between these pathways revealed that the cell cycle and hippo signaling pathways exhibited the closest interaction with other pathways. Meanwhile, crucial genes-gadd45, cyclin D (ccnd), proliferating cell nuclear antigen (pcna), yorkie homologue (yap1, also known as yap), and tafazzin (taz) simultaneously participated in multiple pathways were also noted (Figure 3C).

Moreover, the WGCNA network generated using 60 hub genes in the turquoise module (GS. length and GS. weight > 0.50, MM > 0.90, and Weight > 0.60) revealed that yap1 shared the closest relationship with other genes (Figure 2C). Likewise, DNA methyltransferase 3b (dnmt3b) displayed the closest relationship with other genes in the brown network constructed using 71 hub genes (GS. length and GS. weight < −0.35, MM > 0.90, and Weight > 0.48) (Figure 2D).

The genes involved in the cell cycle and hippo signaling included gadd45g, cyclin-dependent kinase 7 (cdc7), pcna, minichromosome maintenance proteins 7 (mcm7), yap1, bmp2, origin recognition complex subunit 4 (orc5), and cyclin B1 (ccnb1). The other genes contained ubiquitin-conjugating enzyme E2L3 (ube2l3), cell death activator CIDE-3 (cidec), and ubiquitin C (ubc), which were all selected for expression pattern validation by qPCR (Primer information listed in Supplementary Table S4). Similar expression patterns of DEGs in transcriptomes were observed in the four tissues from three comparison groups (Figure 4C).

Transcriptome analysis showed many DEGs (10,766 (∼92%) and 1,439 (61–82%) in the gonad and muscle, respectively) overlapped when male and pseudomale were compared with the female group which indicated that males and pseudomales had similar transcriptome patterns, especially in the gonad. Besides, the methylome analysis showed that the DNA methylation levels were higher in the male and pseudomale than in the female groups. Previous studies also reported a similar comparison in the gonad (Shao et al., 2014).

This study found that there were negative correlations in the 2 kb upstream and 2 kb downstream. However, the negative relationship gradually changed to positive in the whole gene body. Studies in mammals and plants have shown that methylation in the promoter and downstream region negatively correlates with the transcript expression (Jones, 1999; Liang et al., 2014). Inversely, the methylation at the TSS inhibits transcript initiation, while the methylation in the gene body can stimulate transcript elongation or affect splicing (Jones, 1999; Jones, 2012). The active transcription is positively associated with gene body methylation in the active X chromosome (Hellman and Chess, 2007).

Herein, transcriptome analysis revealed the significant female-biased expression of most cell cycle-related genes in the gonad and muscle of C. semilaevis, implying their essential roles in the female-biased SSD phenomenon. Undoubtedly, cell cycle is a ubiquitous and complex process that promotes cell genome duplication, growth, and division (Poon, 2016). Cell growth in yeast, bacteria, and plants are usually accompanied with the promotion of cell cycle progression (Goranov et al., 2009; Sablowski and Carnier Dornelas, 2014; Levin and Taheri-Araghi, 2019). Members of cyclin-cdk complexes-involved in regulating cell cycle transition through G1, S, G2, M and G0 phases (Schafer, 1998), commonly displayed upregulated transcripts and hypomethylated upstream or downstream region in the female gonad and muscle. Cyclin D is an effective mitogenic sensor that increases expression to convey extracellular signals to the cell cycle via growth factors. The anaphase-promoting complex (APC 4, 5, 6, 10, 13, CDCs) is involved in the proteolysis of cell cycle regulatory proteins (Teixeira and Reed, 2013). E2F-DP complexes (E2F1-2) promote transcription of downstream genes, including cyclin A and E (Poon, 2016). Similarly, other crucial genes-pcna, gadd45, orc1-5, mcm4, 6, 7, and smad2 also exhibited female-biased transcripts and hypomethylated levels at the 2 kb upstream and downstream regions of the gonad and muscle tissues. However, detailed mechanisms of these cell cycle-related genes involved in growth regulation need further exploration.

The hippo pathway-related genes in the negative growth-related brown module suggested their possible roles in the sexual growth difference of C. semilaevis. The hippo pathway (conserved in Drosophila melanogaster and mammals) has been shown to regulate cell growth, fate decision, organ size, and regeneration (Halder and Johnson, 2011; Ma et al., 2019; Cho and Jiang, 2021). Besides, tissue overgrowth phenotypes in D. melanogaster are promoted by the following factors: 1) The mutation of its core kinase components, such as warts (wts) (Justice et al., 1995; Xu et al., 1995), hippo (hpo) (Harvey et al., 2003), and salvador (sav) (Tapon et al., 2002). 2) the loss-of-function of its upstream regulatory factors, including merlin (mer), expanded (ex) (Hamaratoglu et al., 2006), fat (ft) (Cho et al., 2006). 3) Overexpression of its downstream effectors, such as yorkie (yki) and scalloped (sd) (Huang et al., 2005; Wu et al., 2008). The overexpression of yorkie homologue-yap also enhances liver size in mouse by increasing cell proliferation (Dong et al., 2007). Herein, upstream factor, mer (nf2a and nf2b), were downregulated in the female gonad of C. semilaevis, and downstream effector-sd (tead3/tef5) were upregulated. Another key downstream factor, yap1 had two homologues in chromosomes 4 and 19, which separately displayed male- or pseudomale-biased transcripts in the gonad, and female-biased expression levels in the gonad and muscle, indicating the complexity of the hippo pathway involved in fish growth regulation. The integration of transcriptome and methylome data suggested that the inhibition of hippo signaling via the canonical and non-canonical pathways was possibly regulated by hypomethylated patterns in the gene body region, especially the gonad.

The tissue overgrowth phenotype induced by the hippo pathway in Drosophila are commonly accompanied and characterized by the activation of cell cycle regulator cyclin e (Harvey et al., 2003). Recent studies in various mammal cells have revealed that cell cycle are associated with the hippo pathway. For instance, yap knockdown causes cell cycle arrest at the G0/G1 phase by decreasing the cycle-related protein-Cyclin D1 (Chen et al., 2017). In contrast, yap activation promotes cell cycle progression by increasing Cyclin D1 (Mizuno et al., 2012; Xie et al., 2019). Besides cyclin e and cyclin d, other cell cycle genes, including cdc6 and e2f1, are also downstream targets of yap (Kapoor et al., 2014; Kim et al., 2019). Moreover, the cell cycle controlling factor, APC/Cyclosome (APC/C)^cdh1 E3 ubiquitin ligase complex, can increase YAP/TAZ activities by promoting LATS1/2 kinase degradation (Kim et al., 2019). Herein, the possible downstream targets, including cyclin d1, cyclin e1, cyclin e2, cdc6, and e2f1, the possible upstream, apc/c and cdh1 genes exhibited female-biased mRNA levels in the gonad of C. semilaevis.

Additionally, the GH-IGF system, highly conserved in somatotropic axis across various species, including mammals and fish (Baker et al., 1993; Reinecke et al., 2000; Ahmed and Farquharson, 2010), participate in growth regulation via JAK-STAT, PI3K-Akt signaling etc. (Brooks et al., 2014; Ipsa et al., 2019; Rotwein, 2020). IGF1, a crucial factor of the GH-IGF system, can inhibit the hippo pathway via PI3K-Akt signaling, thus activating YAP/TAZ effector to promote cell proliferation (Fan et al., 2013; Fruman et al., 2017; Chen et al., 2018; Thompson, 2020). Herein, igf1 had upregulated transcripts and hypomethylated patterns in the female gonad of C. semilaevis. Thus, this phenomenon seems to provide a possible linkage to hippo pathway and classical somatotropic axis regulators.

The de novo and maintenance of DNA methylation in mammals require DNA methyltransferase enzymes, including DNMT1, DNMT3A, and DNMT3B (Reik et al., 2001; Jones and Liang, 2009). Herein, the pattern of transcripts levels of dnmt3a and dnmt3b were FG < MG < PMG, indicating a negative correlation with the DNA methylation levels. Besides, the histone variant H2A.Z, strongly antagonistic to DNMTs (Zilberman et al., 2008; Conerly et al., 2010), had two orthologue genes in C. semilaevis chromosome w and z with 100% protein homology. Their total expression levels were FG > MG > PMG, while their DNA methylation levels were FG < MG < PMG. Recent studies have revealed that the recruitment of DNMT1 to replication sites occurs via an interaction with PCNA and UHRF1 (Jurkowska and Jeltsch, 2016; Jimenji et al., 2019), which both had the highest expression levels in the female gonad and muscle of C. semilaevis. Interestingly, uhrf1 overexpression can cause DNA hypomethylation, similar to uhrf1 mutants (Kent et al., 2016). Therefore, its female-biased transcripts could be responsible for the hypomethylation of the gonad and muscle tissues in the female group of C. semilaevis. Surprisingly, the highest expression levels of dnmt1 were detected in the female groups. Therefore, further investigation is required to confirm whether a negative interaction exists between the dnmt1 and DNA methylation levels.

In conclusion, this study identified many DEGs implicated in various cell growth and death-related pathways through transcriptome analysis of the four organ tissues from three sexual groups. Methylome data revealed that the DNA methylation patterns in the male and pseudomale tissues were higher than in the females. Moreover, the integration of transcriptome and methylome data revealed that cell cycle-related genes were upregulated and hypomethylated in the upstream or downstream regions of the female gonad and muscle. Conversely, the male and pseudomale-biased expression of hippo signaling pathway genes were positively correlated with their hypermethylation levels in the gene body. The activation of the cell cycle and the inhibition of the hippo signaling pathway could be implicated in the C. semilaevis female-biased SSD. Additional functional experiments involving essential genes of these pathways and the epigenetic regulatory factors are needed to assess the fish sexual size dimorphism.

In this study, individual fish were anesthetized with MS-222 before the sampling to relieve pain. Four somatotropic and reproductive tissues from the brain, liver, gonad, and muscle were isolated from nine females, nine males, nine pseudomales in triplicate and immediately stored in liquid nitrogen. The fish, each 1.5-years-old C. semilaevis species were cultivated in Haiyang Yellow Sea Fisheries Limited Company. The genetic sex identification of the male and the pseudomale was performed using the primers sex-F (CCT​AAA​TGA​TGG​ATG​TAG​ATT​CTG​TC) and sex-R (GAT​CCA​GAG​AAA​ATA​AAC​CCA​GG), described in a previous study (Liu Y. et al., 2014).

Each triplicate sample was pooled per treatment into one sample, generating 36 samples hereafter named FB1-3, MB1-3, PMB1-3, FG1-3, MG1-3, PMG1-3, FL1-3, ML1-3, PML1-3, FM1-3, MM1-3, and PMM1-3. Total RNA was extracted using the Trizol reagent kit (Invitrogen, United States) according to the manufacturer’s protocol.

The RNA quality was assessed on an Agilent 2,100 Bioanalyzer (Agilent Technologies, United States) and by agarose gel electrophoresis, then the total RNAs with RIN > 7.0 were used for library construction. Briefly, mRNA and ncRNA were retained by removing rRNAs using Ribo-ZeroTM Magnetic Kit (Epicentre, United States). The cDNA fragments were reverse-transcribed using DNA polymerase I and purified with a QiaQuick PCR extraction kit (Qiagen, Venlo, Netherlands). After end repair, poly(A) was added and ligated onto Illumina sequencing adapters. The cDNA libraries were sequenced on Illumina HiSeqTM 4,000 (Guangzhou Genedenovo Biotechnology Co. Ltd., China).

Clean reads were obtained by removing low-quality bases using fastp 0.18.0 software, while rRNA was removed using Bowtie 2.2.8. The remaining paired-end cleaned reads were mapped to the C. semilaevis reference genome (Cse_v1.0) using HISAT2 2.1.0 for transcript reconstruction and novel transcripts identification. FPKM (fragment per kilobase of transcript per million mapped reads) values were calculated using StringTie 1.3.1 to quantify expression abundance for the transcription region.

The differentially expressed genes (DEGs) were identified using DESeq2 software. All DEGs were mapped subjected to Gene Ontology (GO) analysis on the GO database (http://www.geneontology.org/) and KEGG database (http://www.kegg.jp), respectively, to infer functions of the DEGs. Gene numbers were calculated for every term and significantly enriched GO terms/KEGG pathways in DEGs compared to the genome background, as defined by hypergeometric tests. The GO terms or KEGG pathways with corrected FDR < 0.05 were considered significantly enriched.

Twelve DEGs (Supplementary Table S4) were selected for the qRT-PCR validation as previously described (Wang et al., 2018) using the C. semilaevis β-actin as the internal reference. Total RNA was reverse-transcribed into cDNA using the PrimeScript RT reagent Kit with gDNA Eraser (Takara Bio, Japan). The qPCR reaction was conducted using SYBR Premix Ex Taq (Takara Bio, Japan) in a 20 μL reaction on the ABI 7500 Fast Real-time PCR System (Applied Biosystems, United States). Fold change values for relative expression levels of the 12 genes were calculated using the 2^−ΔΔCt method (Livak and Schmittgen, 2001). The log2FC values retrieved by qPCR and RNA-seq counts were used for graphical presentation by GraphPad Prism 8.

Co-expression networks were constructed using WGCNA (v1.47) to describe the correlation patterns among DEGs across multiple samples and find modules of highly correlated to female-biased SSD in the R package following a previous procedure (Langfelder and Horvath, 2008). Briefly, DEGs were submitted for co-expression modules construction with the power = 8 and minModuleSize = 50. GO and KEGG pathway enrichment analysis were conducted for genes in each module to analyze the biological functions. Subsequently, module eigengenes were analyzed to determine the correlation relationship between modules and the growth performance of samples. Modules with the highest/lowest correlation values and p < 0.05 were considered positive-/negative-related, while the genes manifesting the closest relationships with other genes were recognized as hub genes. The network for hub genes characterization was generated by Gephi 0.9.2 (Bastian et al., 2009).

For interpreting the possible epigenetic mechanism involved in sexual size dimorphism, gonad and muscle tissues with more DEGs between different sexual individuals were submitted for whole-genome DNA methylation sequence analysis. The integrity of the extracted DNA was determined using NanoPhotometer^®, spectrophotometer (IMPLEN, United States), and agarose gel electrophoresis. Thereafter, the DNA was fragmented into 100–300 bp by Sonication (Covaris, United States) and purified using the MiniElute PCR Purification Kit (QIAGEN, United States). The purified DNA fragments were ligated to methylate sequencing adapters after end repair with the “A” nucleotide. Subsequently, bisulfite conversion was conducted using a Methylation-Gold kit (ZYMO, United States) to convert unmethylated cytosine to uracil. Finally, the converted DNA fragments were PCR amplified and sequenced using Illumina HiSeqTM 2,500 at Gene Denovo Biotechnology Co. (Guangzhou, China).

Raw reads were filtered by removing reads containing >10% unknown nucleotides and low-quality reads containing >40% low-quality (Q-value ≤ 20) bases. The cleaned, filtered reads were mapped to the C. semilaevis reference genome using BSMAP 2.90. The methylated cytosines were called using a Perl script and tested using a previous algorithm (Lister et al., 2009). The methylation level was calculated based on the percentage of methylated cytosine in the whole genome, on each chromosome, and at different genomic regions for each sequence context (CG, CHG, and CHH). The methylation profiles at flanking 2 kb regions and gene body were plotted based on the average methylation levels for each window to assess distinct methylation patterns in different genomic regions.

DMRs between two samples were identified using methylkit V1.4.1 with calling window size of 200 bp and the minimum read coverage of 4 bp. DMRs for each sequence context (CG, CHG, and CHH) based on this criteria: 1) For CG and CHG, numbers in each window ≥5, absolute values of the difference in methylation ratio ≥ 0.25, and q ≤ 0.05; 2) For CHH, numbers in a window ≥15, absolute values of the difference in methylation ratio ≥0.15, and q ≤ 0.05; 3) For all C, numbers in a window ≥20, absolute values of the difference in methylation ratio ≥0.2, and q ≤ 0.05. Subsequently, GO and KEGG pathway enrichment analysis were conducted for DMGs.

Genes were categorized into four classes based on their expression levels to determine whether gene expression influences DNA methylation. These classes included non-expressed (RPKM ≤ 1), low-expressed (1 < RPKM ≤ 10), middle-expressed (10 < RPKM ≤ 100), and high-expressed groups (RPKM > 100).

Spearman’s correlation analysis was performed to discern the statistical relationships between DNA methylation and gene expression within ± 2 kb flanking regions and the gene body. A positive correlation was indicated by rho >0, while rho <0 indicated a negative correlation.

Common genes between DMGs and DEGs were extracted, and their GO/KEGG enrichment analysis was conducted to explore the potential functions of DNA methylation responsible for DEGs.

The genes from the positive or negative growth correlated modules were separately selected for the integration by combining their transcripts and DNA methylation levels within the R software. The value for DEG-selection was set at |log2FC| > 1, and the value for DMG-selection was set at |meth.diff| > 25.