Integration of Transcriptome and Methylome Highlights the Roles of Cell Cycle and Hippo Signaling Pathway in Flatfish Sexual Size Dimorphism

Sexual size dimorphism (SSD) is the difference in segments or body size between sexes prevalent in various species. Understanding the genetic architecture of SSD has remained a significant challenge owing to the complexity of growth mechanisms and the sexual influences among species. The Chinese tongue sole (Cynoglossus semilaevis), which exhibits a female-biased SSD and sex reversal from female to pseudomale, is an ideal model for exploring SSD mechanism at the molecular level. The present study aimed to integrate transcriptome and methylome analysis to unravel the genetic and epigenetic changes in female, male, and pseudomale C. semilaevis. The somatotropic and reproductive tissues (brain, liver, gonad, and muscle) transcriptomes were characterized by RNA-seq technology. Transcriptomic analysis unravelled numerous differentially expressed genes (DEGs) involved in cell growth and death-related pathways. The gonad and muscle methylomes were further employed for screening differentially methylated genes (DMGs). Relatively higher DNA methylation levels were observed in the male and pseudomale individuals. In detail, hypermethylation of the chromosome W was pronounced in the pseudomale group than in the female group. Furthermore, weighted gene co-expression network analysis showed that turquoise and brown modules positively and negatively correlated with the female-biased SSD, respectively. A combined analysis of the module genes and DMGs revealed the female-biased mRNA transcripts and hypomethylated levels in the upstream and downstream regions across the cell cycle-related genes. Moreover, the male and pseudomale-biased gene expression in the hippo signaling pathway were positively correlated with their hypermethylation levels in the gene body. These findings implied that the activation of the cell cycle and the inhibition of the hippo signaling pathway were implicated in C. semilaevis female-biased SSD. In addition, the dynamic expression pattern of the epigenetic regulatory factors, including dnmt1, dnmt3a, dnmt3b, and uhrf1, among the different sexes correspond with their distinct DNA methylation levels. Herein, we provide valuable clues for understanding female-biased SSD in C. semilaevis.

Introduction

Sexual size dimorphism (SSD) refers to the differences in segment or body size between sexes prevalent in various species including drosophilid flies, mammals, birds, reptiles, arachnida, and fish (Cox et al., 2007; Foellmer and Moya-Larano, 2007; Lindenfors et al., 2007; Székely et al., 2007; Mei and Gui, 2015; Halvorsen et al., 2016; McLean et al., 2018). Understanding the genetic architecture of SSD is a major challenge considering the complexity of the growth mechanisms and the sexual influences. In Drosophila, the diminutive (myc) gene located on the X-chromosome is crucial in the determination of SSD by contributing to the activation of a female-specific gene transformer (tra) (Mathews et al., 2017). In primates, the sexual brain size dimorphism is regulated by estrogen, suppressing the promoter activities of microcephaly genes (Shi et al., 2015). Similarly, the growth rate in reptiles is influenced by testosterone levels and environmental factors such as temperature and precipitation (Stillwell and Fox, 2007; Cox et al., 2009; Michael et al., 2014; Agha et al., 2018).

In fish species, a female-biased or a male-biased SSD causes a growth disadvantage within a single-sex, subsequently confining the sustainable development of aquaculture. For example, the Chinese tongue sole (Cynoglossus semilaevis), a typical female heterogamete (ZW/ZZ) flatfish, exhibits a female-biased SSD and sex reversal from female (ZW) to pseudomale (ZW), which significantly increases the proportion of phenotypic males in the aquaculture (Zhou et al., 2005; Chen et al., 2009; Chen et al., 2014). Previous studies of SSD in C. semilaevis have investigated transcriptomic patterns in both female and male individuals (Wang et al., 2014; Wang et al., 2016; Wang et al., 2018) and various growth-related genes including growth hormone (gh), growth hormone receptor (ghr), pituitary adenylate cyclase activating polypeptide (pacap), growth hormone releasing hormone (ghrh), leptin, bone morphogenetic factor 6 (bmp6), bmp7, and growth arrest and DNA damage inducible gamma (gadd45g) have also been identified by expression levels analysis (Ji et al., 2011; Ma et al., 2012a; Ma et al., 2012b; Liu W.-J. et al., 2014; Ma et al., 2017; Xu et al., 2018). However, the mechanism of SSD in pseudomale is still unknown.

Despite the female and pseudomale having the same genetic background, they exhibit substantial growth differences, which their epigenetic mechanism could explain. Epigenetic modifications such as DNA methylation play crucial roles in gene expression, histone modification, and maintenance of genome stability (Cedar and Bergman, 2009; Cedar and Bergman, 2012). DNA methylation regulates plant species’ growth and development processes (Bartels et al., 2018) and the body weight of male and female mammals (Gonzalez-Nahm et al., 2018). In Nile tilapia (Oreochromis niloticus), the epigenetic status of gh promoter is negatively correlated with male growth superiority (Zhong et al., 2014). However, in C. semilaevis, the relationship between DNA methylation and growth trait was studied only in a few genes, including gh, igf1, ghr1, and pacap (Zhao et al., 2015; Si et al., 2016; Li et al., 2017). Much of the research on the whole-genome DNA methylation in C. semilaevis has focused on its sexual reversal mechanism and regulation of disease resistance pathways (Shao et al., 2014; Xiu et al., 2019). Moreover, whole-genome DNA methylation and its effects on the SSD transcription patterns in C. semilaevis remains unknown.

Thus, the present study aimed to explore the female-biased SSD in C. semilaevis through whole-genome transcriptomic and weighted gene co-expression network analysis in four organs (brain, liver, gonad, and muscle), DNA methylation analysis in the gonad and muscle. Integrating transcriptome and methylome data in female, male and pseudomale individuals, will provide valuable insights into the crucial pathways and genes involved in female-biased SSD of C. semilaevis.

Results

Screening of Differentially Expressed Genes From Female, Male, and Pseudomale C. semilaevis Transcriptomes

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).

FIGURE 1

FIGURE 1. The identification and functional enrichment of DEGs in the whole transcriptome. (A,B) The identification of DEGs (A) and venn diagram demonstrating DEGs (B) from the brain, gonad, liver and muscle of M-vs-F, M-vs-PM, PM-vs-F groups. (C) Description of functional KEGG enrichment for DEGs. The color bar means q-value from low (red) to high (green).

DEGs Were Affiliated to Many Cell Growth and Death-Related Pathways

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.

Turquoise and Brown Modules Exhibiting Positive- or Negative-Correlation to Female-Biased SSD Were Identified by Weighted Gene Co-Expression Network Analysis

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).

FIGURE 2

FIGURE 2. Screening of modules and hub genes by WGCNA. (A) The hierarchical clustering of WGCNA. (B) The relationship between modules and growth traits. (C,D) The gene network constructed by the hub genes in the turquoise module (C) and brown module (D).

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).

FIGURE 3

FIGURE 3. The KEGG enrichment for turquoise and brown modules and the relationship among pathways and genes. (A,B) The top 20 enriched KEGG pathways in the turquoise module (A) and brown module (B). (C) The relationship among the enriched pathways and key genes.

FIGURE 4

FIGURE 4. The heatmap of DEGs from the cell cycle, cellular senescence and hippo signaling pathways. (A,B) The heatmap of gene expression patterns for cell cycle pathway (A), cellular senescence and hippo signaling pathways (B). (C) The qPCR validation of transcriptome gene expression. The color bar means expression levels from low (blue) to high (red).

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).

Expression Pattern Validation of DEGs From Cell Cycle and Hippo Signaling Pathways

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).

Whole-Genome Bisulfite Sequencing of Gonad and Muscle From Different Sexual Groups

Methylation Analysis of Gonad and Muscle From Different Sexual Groups

The gonad and muscle tissues with the highest number of DEGs between individuals of different sex were subjected to WGBS analysis to reveal the possible epigenetic mechanism. A total of 2.08 × 10⁹ reads with the sequencing depth of 27.36× were obtained, and 70.26% mapped to the reference genome (Supplementary Table S5). The sequencing depth and the cumulative distribution of effective sequencing depth by C base proved the uniformity and high quality of the WGBS sequencing (Supplementary Figure S10). The efficiency of bisulfite treatment was >99% after the detection with lambda DNA. The effective C base coverage rates by chromosome, different genomic region, and repeat region were 81.67–97.19%, 91.88–98.54%, and 67.72–95.33%, respectively. The whole genomic methylation analysis revealed that the DNA methylation levels in the male and pseudomale tissues were higher than in the females (Figure 5A). Moreover, the 3′UTR and 5′UTR exhibited the highest and lowest DNA methylation levels, respectively (Figure 5A). In the female gonad and muscle, the highest CG methylation levels was both observed on the W chromosome (Figure 5B). Interestingly, chromosome 17 demonstrated the lowest DNA methylation levels in all examined six groups. A detailed comparison (Figures 5B,C) showed that the total methylation levels of chromosome Z in the gonad and muscle were PM > M > F, while the total methylation levels of chromosome W in these two tissues were F < PM. Furthermore, the dynamic patterns of CG methylation in different genomic regions revealed that the highest and lowest CG methylation levels were observed at the start and end of the gene body, respectively (Figure 5D).

FIGURE 5

FIGURE 5. The DNA methylation patterns in three sexual groups by different genomic regions and chromosomes. (A) The DNA methylation patterns in different genomic regions of female, male and pseudomale gonads and muscles. (B) Illustration of DNA methylation levels in different chromosomes of the six groups. (C) Manifestation of methylation status for w and z chromosomes in the different groups. The left Y-axis illustrated mCG density in per 10 kb (mCG number/10 kb), and the right Y-axis presented mCG density per 100 kb by calculating each CG site. (D) The mean methylation levels in the gonad and muscle from three sexual groups. The color bar means methylation levels from low (blue) to high (red).

The expression patterns of dnmt1, dnmt3a, dnmt3b, and ubiquitin like with phd and ring finger domains 1 (uhrf1) were illustrated using IGV software (Supplementary Figure S11) to reveal potential roles of these epigenetic regulatory factors in dynamic methylation levels of different sexual groups. Female-biased dnmt1 and uhrf1 transcripts were observed in the gonad and muscle tissues. Additionally, the expression levels of dnmt3b and dnmt3a were F < M < PM in the gonad and muscle tissues.

The identification of differentially methylated regions (DMRs) between two samples resulted in: 1) 414,401 gonad and 2,727 muscle DMRs from the male-vs-female group, 2) 398,065 gonad and 3,949 muscle DMRs from the pseudomale-vs-female group, and 3) 2,872 gonad and 1,293 muscle DMRs from the male-vs-pseudomale group (Supplementary Figure S12A). Further annotation for DMRs revealed 19,968 and 3,318 differentially methylated genes (DMGs) from the gonad and muscle tissues, respectively (Figure 6A). A total of 1,214 gonad and 78 muscle DMGs were separately overlapped in the three comparison groups. The subsequent KEGG enrichment analysis discovered many pathways related to cellular community, signal transduction, and cancers (Supplementary Figure S12B). More than 95% of the genes in the KEGG enriched pathways exhibited differential methylation levels in the gonad comparisons (Supplementary Figures S13A,B). Meanwhile, less than 20% of the genes in the annotation pathway were differentially methylated in the muscle tissues (Supplementary Figures S13C,D).

FIGURE 6

FIGURE 6. The description of DMGs in the methylome and their linkage analysis within the samples. (A) The Venn diagram of DMGs in the six comparison groups. (B) The relationship between gene expression and DNA methylation of different regions in the samples. Four groups of DEGs with none (fpkm ≤ 1), low (1 < fpkm ≤ 10), middle (10 < fpkm ≤ 100) and high (fpkm > 100) expression were classified and their mean DNA methylation expression level were subsequently calculated in the 2 kb upstream, gene body, and 2 kb downstream regions.

Linkage Analysis of Gene Expression and DNA Methylation Within the Samples and Among the Groups

Four groups of DEGs with no expression (fpkm ≤ 1), low expression (1 < fpkm ≤ 10), moderate expression (10 < fpkm ≤ 100), and high expression (fpkm > 100) were classified. Their mean DNA methylation levels were subsequently calculated in the gene body, upstream and downstream regions. The results revealed negative correlations within 2 kb upstream and 2 kb downstream regions. While, the negative relationship gradually turned positive in the whole gene body region (Figure 6B).

The Spearman’s correlation coefficients between DNA methylation and gene expression within the samples were shown in Supplementary Figure S14. A negative correlation was observed in ±2 kb flanking regions and the gene body.

Common genes between DEGs and DMGs were screened, and their expression patterns were subjected to cluster and enrichment analysis to uncover the contribution of differential DNA methylations to transcript regulation. Results showed that 86.21, 7.89, and 85.83% DEGs in the MG-vs-FG, MG-vs-PMG, and PMG-vs-FG displayed differential DNA methylation (Supplementary Figure S15). Meanwhile, only 6.59, 9.69, and 7.55% DEGs from the MM-vs-FM, MM-vs-PMM and PMM-vs-FM were differentially methylated. In total, 11,087 and 278 overlapped genes were separately screened from the gonad and muscle (Supplementary Figure S15). Subsequently, negative and positive pairs were separately screened from M-vs-F, PM-vs-F, and M-vs-PM groups (Supplementary Figure S15B).

Nine-Quadrant Diagram for Module Genes

The DEGs from turquoise and brown modules were separately employed for quadrant diagram construction to understand the potential correlation between mRNA expression and DNA methylation levels. In the turquoise module, 1,420 genes from the MG-vs-FG group and 1,450 genes from the PMG-vs-FG group displayed upregulated transcripts and 2 kb upstream hypomethylation patterns in the female gonad tissues (Figure 7A, quadrant 1). Moreover, these genes were significantly classified into splicesome, DNA replication, ribosome biogenesis, and cell cycle pathways (q < 0.05) (Figure 7B). By analysis of 2 kb downstream methylation patterns, 1,551 and 1,568 genes exhibited upregulated transcripts and hypomethylated levels in MG-vs-FG and PMG-vs-FG groups. Among these genes enriched pathways, DNA replication and cell cycle stand out again (q < 0.05) (Figure 7B).

FIGURE 7

FIGURE 7. The nine quadrant diagram for turquoise module genes by combining transcripts and DNA methylation levels, and the functional classification. (A) The nine quadrant diagram for turquoise module genes and DMGs from 2 kb upstream and 2 kb downstream regions in MG-vs-FG and PMG-vs-FG. (B) The functional classification of quadrant 1 (upper left) genes by 2 kb upstream and 2 kb downstream.

In the brown module, 1,435 genes from MG-vs-FG and 1,422 genes from PMG-vs-FG demonstrated downregulated transcripts and gene body hypomethylation patterns in the female gonad (Figure 8A, quadrant 7). Subsequently, these genes were significantly enriched in hippo signaling pathway and cellular senescence (p < 0.05) (Figure 8B).

FIGURE 8

FIGURE 8. The nine quadrant diagram for brown module genes by combining transcripts and DNA methylation levels, and the functional classification. (A) The nine quadrant diagram for brown module genes and DMGs from gene body regions in MG-vs-FG and PMG-vs-FG. (B) The functional classification of quadrant 7 (upper down) genes by gene body.

The Recognition of Cell Cycle and Hippo Signaling Pathway from the Integration of DMR and DEGs

Based on the nine quadrant diagram results, 38 genes with upregulated transcripts were screened from the cell cycle of the female gonad and muscle tissues. The crucial effectors of cell cycle pathway, cyclin-dependent kinases (cdks)/cyclin, anaphase-promoting complex (apc)/cdcs, E2F Transcription Factor (e2f), gadd45, pcna, orcs, and mcms (Figure 9A), were also involved. Remarkably, their 2 kb upstream and downstream regions demonstrated hypomethylation levels in the female gonad and muscle tissues (Figure 9B).

FIGURE 9

FIGURE 9. The female-biased expression and hypomethylation of cell cycle genes in the female group. (A) Illustration of the cell cycle pathway and important genes. (B) The heatmap of 38 cell cycle-related genes in the gonad and muscle by the transcripts FPKM and mean differentially methylation levels in the 2 kb upstream or downstream (*) region. The color bar means expression or methylation levels from low (blue) to high (red).

In the hippo signaling pathway (Figure 10A), 28 genes exhibited downregulated transcripts in the female gonad tissues (Figure 10B), except for the yap1 (chr19) and transcriptional Enhancer Factor 5 (tef5) (also known as TEA domain transcription factor 3, tead3). A total of 30 genes comprised the core hippo signaling factors-yap1, mer (nf2a and nf2b), and tef5. Additionally, non-canonical hippo signaling pathway-bmps, wnts, smad, and tcf/lef (Figure 10A) also exhibited downregulated transcripts in the female gonad. Consistent with WGBS data, the expressions patterns of the hippo pathway-related genes of the gonad tissues positively correlated with DNA methylation in the gene body.

FIGURE 10

FIGURE 10. Inhibition of hippo signaling pathway and their gene body hypomethylation in the female gonad. (A) Illustration of hippo signaling pathway and essential genes. (B) The heatmap of 30 hippo pathway-related genes in the gonad by RNA transcripts FPKM and mean differentially methylation levels in the gene body region. The color bar means expression or methylation levels from low (blue) to high (red).

Discussion

SSD is a common phenomenon in many species and has received much attention on its evolutionary drivers. Although the Rensch’s Rule states that body size variation increases in male-biased SSD and decreases in female-biased SSD (Rensch, 1950), some studies have revealed that many organisms do not conform to this rule (Webb and Freckleton, 2007; Liao et al., 2013; Martin et al., 2017). To date, three hypotheses have been proposed to explain sexual size allometry: evolutionary constraints, natural selection, and sexual selection (Dale et al., 2007; McLean et al., 2018). However, the molecular mechanism of SSD in C. semilaevis is complex due to the reversal of female-to-pseudomale and the resembling growth performance between pseudomale and male. Herein, transcriptome and methylome integration was used to assess the female-biased SSD phenomenon in C. semilaevis.

Transcriptome and Methylome Patterns in the Male and Pseudomale Groups

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).

The Relationship Between Transcriptome and Methylome of C. semilaevis

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).

The Upregulated mRNA and Hypomethylated of Cell Cycle-Related Genes in the Female Group Might Cause Female-Biased SSD in C. semilaevis

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 Relationship Between Hippo Signaling and Negative Growth

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 Relationship Among Cell Cycle, Hippo Signaling, and Other Pathways

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 Possible Epigenetic Regulatory Mechanism in C. semilaevis

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.

Materials and Methods

Fish Sampling, Sex Identification, Tissue Collection, and RNA Isolation

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.

Library Construction and Sequencing

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).

Read Processing and Analysis of Differentially Expressed Genes

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.

Validation of DEGs by qRT-PCR

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.

The Weighted Gene Co-expression Network Analysis for DEGs

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).

Genomic DNA Library Construction, Sequencing, and Data Filtering

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).

Methylation Level Analysis

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.

The Identification of Differentially Methylated Regions and Functional Enrichment Analysis for Differentially Methylated Genes

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.

The Correlation of DNA Methylation and Gene Expression in Samples

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.

Correlation of DMGs and DEGs Between Groups

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 Nine Quadrant Diagram for Module Genes by Combining Their Transcripts and DNA Methylation Levels

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.

Data Availability Statement

The whole transcriptome and methylome data presented in the study are deposited in the NCBI SRA repository, accession number PRJNA743138 and PRJNA741107, respectively.

Ethics Statement

The animal study was reviewed and approved by The Animal Care and Use Committee at the Chinese Academy of Fishery Sciences.

Author Contributions

NW and SC conceived and designed the experiments. NW, WX, and YY performed fish tissue sampling. NW and JW conducted detection of pseudomale and male. NW, QY, JW, and RS carried out the qRT-PCR experiment. NW analyzed the data and wrote the paper. WX, ML, JG, and YC offered valuable suggestions in the article. We are grateful to Guangzhou Genedenovo Biotechnology Co., Ltd. for assisting in sequencing and bioinformatics analysis. All authors read and approved the final manuscript.

Funding

This work was supported by the National Key R&D Program of China (2018YFD0900205), the National Natural Science Foundation of China (31873037, 31730099), Central Public-interest Scientific Institution Basal Research Fund, CAFS (2020TD20) and the Taishan Scholar Project of Shandong Province.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.743722/full#supplementary-material

References

Agha, M., Ennen, J. R., Nowakowski, A. J., Lovich, J. E., Sweat, S. C., and Todd, B. D. (2018). Macroecological Patterns of Sexual Size Dimorphism in Turtles of the World. J. Evol. Biol. 31, 336–345. doi:10.1111/jeb.13223

CrossRef Full Text | Google Scholar

Ahmed, S. F., and Farquharson, C. (2010). The Effect of GH and IGF1 on Linear Growth and Skeletal Development and Their Modulation by SOCS Proteins. J. Endocrinol. 206, 249–259. doi:10.1677/joe-10-0045

PubMed Abstract | CrossRef Full Text | Google Scholar

Baker, J., Liu, J.-P., Robertson, E. J., and Efstratiadis, A. (1993). Role of Insulin-like Growth Factors in Embryonic and Postnatal Growth. Cell 75, 73–82. doi:10.1016/S0092-8674(05)80085-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartels, A., Han, Q., Nair, P., Stacey, L., Gaynier, H., Mosley, M., et al. (2018). Dynamic DNA Methylation in Plant Growth and Development. Ijms 19, 2144. doi:10.3390/ijms19072144

PubMed Abstract | CrossRef Full Text | Google Scholar

Bastian, M., Heymann, S., and Jacomy, M. (2009). “Gephi: an Open Source Software for Exploring and Manipulating Networks,” in International AAAI Conference on Weblogs and Social Media. doi:10.13140/2.1.1341.1520

CrossRef Full Text | Google Scholar

Brooks, A. J., Dai, W., O’Mara, M. L., Abankwa, D., Chhabra, Y., Pelekanos, R. A., et al. (2014). Mechanism of Activation of Protein Kinase JAK2 by the Growth Hormone Receptor. Science 344, 1249783. doi:10.1126/science.1249783

PubMed Abstract | CrossRef Full Text | Google Scholar

Cedar, H., and Bergman, Y. (2009). Linking DNA Methylation and Histone Modification: Patterns and Paradigms. Nat. Rev. Genet. 10, 295–304. doi:10.1038/nrg2540

PubMed Abstract | CrossRef Full Text | Google Scholar

Cedar, H., and Bergman, Y. (2012). Programming of DNA Methylation Patterns. Annu. Rev. Biochem. 81, 97–117. doi:10.1146/annurev-biochem-052610-091920

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., You, H., Li, Y., Xu, Y., He, Q., and Harris, R. C. (2018). EGF Receptor-dependent YAP Activation Is Important for Renal Recovery from AKI. Jasn 29, 2372–2385. doi:10.1681/asn.2017121272

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Wang, J., Yao, S.-F., Zhao, Y., Liu, L., Li, L.-W., et al. (2017). Effect of YAP Inhibition on Human Leukemia HL-60 Cells. Int. J. Med. Sci. 14, 902–910. doi:10.7150/ijms.19965

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S.-L., Tian, Y.-S., Yang, J.-F., Shao, C.-W., Ji, X.-S., Zhai, J.-M., et al. (2009). Artificial Gynogenesis and Sex Determination in Half-Smooth Tongue Sole (Cynoglossus Semilaevis). Mar. Biotechnol. 11, 243–251. doi:10.1007/s10126-008-9139-0

CrossRef Full Text | Google Scholar

Chen, S., Zhang, G., Shao, C., Huang, Q., Liu, G., Zhang, P., et al. (2014). Whole-genome Sequence of a Flatfish Provides Insights into ZW Sex Chromosome Evolution and Adaptation to a Benthic Lifestyle. Nat. Genet. 46, 253–260. doi:10.1038/ng.2890

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, E., Feng, Y., Rauskolb, C., Maitra, S., Fehon, R., and Irvine, K. D. (2006). Delineation of a Fat Tumor Suppressor Pathway. Nat. Genet. 38, 1142–1150. doi:10.1038/ng1887

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, Y. S., and Jiang, J. (2021). Hippo-Independent Regulation of Yki/Yap/Taz: A Non-canonical View. Front. Cel. Dev. Biol. 9, 658481. doi:10.3389/fcell.2021.658481

CrossRef Full Text | Google Scholar

Conerly, M. L., Teves, S. S., Diolaiti, D., Ulrich, M., Eisenman, R. N., and Henikoff, S. (2010). Changes in H2A.Z Occupancy and DNA Methylation during B-Cell Lymphomagenesis. Genome Res. 20, 1383–1390. doi:10.1101/gr.106542.110

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, R. M., Butler, M. A., and John-Alder, H. B. (2007). “The Evolution of Sexual Size Dimorphism in Reptiles,” in Sex, Size and Gender Roles: Evolutionary Studies of Sexual Size Dimorphism. Editors D. J. Fairbairn, W. U. Blanckenhorn, and T. Székely (Oxford University Press), 38–49. doi:10.1093/acprof:oso/9780199208784.003.0005

CrossRef Full Text | Google Scholar

Cox, R. M., Stenquist, D. S., and Calsbeek, R. (2009). Testosterone, Growth and the Evolution of Sexual Size Dimorphism. J. Evol. Biol. 22, 1586–1598. doi:10.1111/j.1420-9101.2009.01772.x

CrossRef Full Text | Google Scholar

Dale, J., Dunn, P. O., Figuerola, J., Lislevand, T., Székely, T., and Whittingham, L. A. (2007). Sexual Selection Explains Rensch's Rule of Allometry for Sexual Size Dimorphism. Proc. R. Soc. B. 274, 2971–2979. doi:10.1098/rspb.2007.1043

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, J., Feldmann, G., Huang, J., Wu, S., Zhang, N., Comerford, S. A., et al. (2007). Elucidation of a Universal Size-Control Mechanism in Drosophila and Mammals. Cell 130, 1120–1133. doi:10.1016/j.cell.2007.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, R., Kim, N.-G., and Gumbiner, B. M. (2013). Regulation of Hippo Pathway by Mitogenic Growth Factors via Phosphoinositide 3-kinase and Phosphoinositide-dependent Kinase-1. Proc. Natl. Acad. Sci. 110, 2569–2574. doi:10.1073/pnas.1216462110

PubMed Abstract | CrossRef Full Text | Google Scholar

Foellmer, M. W., and Moya-Laraño, J. (2007). “Sexual Size Dimorphism in Spiders: Patterns and Processes,” in Sex, Size and Gender Roles: Evolutionary Studies of Sexual Size Dimorphism. Editors D. J. Fairbairn, W. U. Blanckenhorn, and T. Székely (Oxford University Press), 71–82. doi:10.1093/acprof:oso/9780199208784.003.0008

CrossRef Full Text | Google Scholar

Fruman, D. A., Chiu, H., Hopkins, B. D., Bagrodia, S., Cantley, L. C., and Abraham, R. T. (2017). The PI3K Pathway in Human Disease. Cell 170, 605–635. doi:10.1016/j.cell.2017.07.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez-Nahm, S., Mendez, M. A., Benjamin-Neelon, S. E., Murphy, S. K., Hogan, V. K., Rowley, D. L., et al. (2018). DNA Methylation of Imprinted Genes at Birth Is Associated with Child Weight Status at Birth, 1 year, and 3 Years. Clin. Epigenet 10, 90. doi:10.1186/s13148-018-0521-0

CrossRef Full Text | Google Scholar

Goranov, A. I., Cook, M., Ricicova, M., Ben-Ari, G., Gonzalez, C., Hansen, C., et al. (2009). The Rate of Cell Growth Is Governed by Cell Cycle Stage. Genes Development 23, 1408–1422. doi:10.1101/gad.1777309

PubMed Abstract | CrossRef Full Text | Google Scholar

Halder, G., and Johnson, R. L. (2011). Hippo Signaling: Growth Control and beyond. Development 138, 9–22. doi:10.1242/dev.045500

PubMed Abstract | CrossRef Full Text | Google Scholar

Halvorsen, K. T., Sørdalen, T. K., Durif, C., Knutsen, H., Olsen, E. M., Skiftesvik, A. B., et al. (2016). Male-biased Sexual Size Dimorphism in the Nest Building Corkwing Wrasse (Symphodus Melops): Implications for a Size Regulated Fishery. ICES J. Mar. Sci. 73, 2586–2594. doi:10.1093/icesjms/fsw135

CrossRef Full Text | Google Scholar

Hamaratoglu, F., Willecke, M., Kango-Singh, M., Nolo, R., Hyun, E., Tao, C., et al. (2006). The Tumour-Suppressor Genes NF2/Merlin and Expanded Act through Hippo Signalling to Regulate Cell Proliferation and Apoptosis. Nat. Cel. Biol. 8, 27–36. doi:10.1038/ncb1339

PubMed Abstract | CrossRef Full Text | Google Scholar

Harvey, K. F., Pfleger, C. M., and Hariharan, I. K. (2003). The Drosophila Mst Ortholog, Hippo, Restricts Growth and Cell Proliferation and Promotes Apoptosis. Cell 114, 457–467. doi:10.1016/s0092-8674(03)00557-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hellman, A., and Chess, A. (2007). Gene Body-specific Methylation on the Active X Chromosome. Science 315, 1141–1143. doi:10.1126/science.1136352

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Wu, S., Barrera, J., Matthews, K., and Pan, D. (2005). The Hippo Signaling Pathway Coordinately Regulates Cell Proliferation and Apoptosis by Inactivating Yorkie, the Drosophila Homolog of YAP. Cell 122, 421–434. doi:10.1016/j.cell.2005.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ipsa, E., Cruzat, V. F., Kagize, J. N., Yovich, J. L., and Keane, K. N. (2019). Growth Hormone and Insulin-like Growth Factor Action in Reproductive Tissues. Front. Endocrinol. 10, 777. doi:10.3389/fendo.2019.00777

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, X.-S., Chen, S.-L., Jiang, Y.-L., Xu, T.-J., Yang, J.-F., and Tian, Y.-S. (2011). Growth Differences and Differential Expression Analysis of Pituitary Adenylate Cyclase Activating Polypeptide (PACAP) and Growth Hormone-Releasing Hormone (GHRH) between the Sexes in Half-Smooth Tongue Sole Cynoglossus Semilaevis. Gen. Comp. Endocrinol. 170, 99–109. doi:10.1016/j.ygcen.2010.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Jimenji, T., Matsumura, R., Kori, S., and Arita, K. (2019). Structure of PCNA in Complex with DNMT1 PIP Box Reveals the Basis for the Molecular Mechanism of the Interaction. Biochem. Biophysical Res. Commun. 516, 578–583. doi:10.1016/j.bbrc.2019.06.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, P. A. (2012). Functions of DNA Methylation: Islands, Start Sites, Gene Bodies and beyond. Nat. Rev. Genet. 13, 484–492. doi:10.1038/nrg3230

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, P. A., and Liang, G. (2009). Rethinking How DNA Methylation Patterns Are Maintained. Nat. Rev. Genet. 10, 805–811. doi:10.1038/nrg2651

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, P. A. (1999). The DNA Methylation Paradox. Trends Genet. 15, 34–37. doi:10.1016/s0168-9525(98)01636-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Jurkowska, R. Z., and Jeltsch, A. (2016). “Enzymology of Mammalian DNA Methyltransferases,” in DNA Methyltransferases - Role and Function. Editors A. Jeltsch, and R. Z. Jurkowska (Cham: Springer International Publishing), 87–122. doi:10.1007/978-3-319-43624-1_5

PubMed Abstract | CrossRef Full Text | Google Scholar

Justice, R. W., Zilian, O., Woods, D. F., Noll, M., and Bryant, P. J. (1995). The Drosophila Tumor Suppressor Gene Warts Encodes a Homolog of Human Myotonic Dystrophy Kinase and Is Required for the Control of Cell Shape and Proliferation. Genes Development 9, 534–546. doi:10.1101/gad.9.5.534

PubMed Abstract | CrossRef Full Text | Google Scholar

Kapoor, A., Yao, W., Ying, H., Hua, S., Liewen, A., Wang, Q., et al. (2014). Yap1 Activation Enables Bypass of Oncogenic Kras Addiction in Pancreatic Cancer. Cell 158, 185–197. doi:10.1016/j.cell.2014.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kent, B., Magnani, E., Walsh, M. J., and Sadler, K. C. (2016). UHRF1 Regulation of Dnmt1 Is Required for Pre-gastrula Zebrafish Development. Developmental Biol. 412, 99–113. doi:10.1016/j.ydbio.2016.01.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, W., Cho, Y. S., Wang, X., Park, O., Ma, X., Kim, H., et al. (2019). Hippo Signaling Is Intrinsically Regulated during Cell Cycle Progression by APC/CCdh1. Proc. Natl. Acad. Sci. USA 116, 9423–9432. doi:10.1073/pnas.1821370116

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinformatics 9, 559. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Levin, P. A., and Taheri-Araghi, S. (2019). One Is Nothing without the Other: Theoretical and Empirical Analysis of Cell Growth and Cell Cycle Progression. J. Mol. Biol. 431, 2061–2067. doi:10.1016/j.jmb.2019.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., He, F., Wen, H., Li, J., Si, Y., Liu, M., et al. (2017). Low Salinity Affects Cellularity, DNA Methylation, and mRNA Expression of Igf1 in the Liver of Half Smooth Tongue Sole (Cynoglossus Semilaevis). Fish. Physiol. Biochem. 43, 1587–1602. doi:10.1007/s10695-017-0395-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, D., Zhang, Z., Wu, H., Huang, C., Shuai, P., Ye, C.-Y., et al. (2014). Single-base-resolution Methylomes of Populus trichocarpa Reveal the Association between DNA Methylation and Drought Stress. BMC Genet. 15, S9. doi:10.1186/1471-2156-15-s1-s9

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, W. B., Zeng, Y., Zhou, C. Q., and Jehle, R. (2013). Sexual Size Dimorphism in Anurans Fails to Obey Rensch's Rule. Front. Zool. 10, 10. doi:10.1186/1742-9994-10-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindenfors, P., Gittleman, J. L., and Jones, K. E. (2007). “Sexual Size Dimorphism in Mammals,” in Sex, Size and Gender Roles: Evolutionary Studies of Sexual Size Dimorphism. Editors D. J. Fairbairn, W. U. Blanckenhorn, and T. Székely (Oxford: Oxford University Press), 16–26. doi:10.1093/acprof:oso/9780199208784.003.0003

CrossRef Full Text | Google Scholar

Lister, R., Pelizzola, M., Dowen, R. H., Hawkins, R. D., Hon, G., Tonti-Filippini, J., et al. (2009). Human DNA Methylomes at Base Resolution Show Widespread Epigenomic Differences. Nature 462, 315–322. doi:10.1038/nature08514

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, W.-J., Zhang, L.-Y., Shao, C.-W., Wang, N., Liu, K., Wen, H.-S., et al. (2014a). Molecular Characterization and Functional Divergence of Two Gadd45g Homologs in Sex Determination in Half-Smooth Tongue Sole (Cynoglossus Semilaevis). Comp. Biochem. Physiol. B: Biochem. Mol. Biol. 177-178, 56–64. doi:10.1016/j.cbpb.2014.09.001

CrossRef Full Text | Google Scholar

Liu, Y., Chen, S., Gao, F., Meng, L., Hu, Q., Song, W., et al. (2014b). SCAR-transformation of Sex-specific SSR Marker and its Application in Half-Smooth Tongue Sole (Cynoglossus Semilaevis). J. Agric. Biotechnol. 22.

Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 25, 402–408. doi:10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Q., Feng, W., Zhuang, Z., and Liu, S. (2017). Cloning, Expression Profiling and Promoter Functional Analysis of Bone Morphogenetic Protein 6 and 7 in Tongue Sole (Cynoglossus Semilaevis). Fish. Physiol. Biochem. 43, 435–454. doi:10.1007/s10695-016-0298-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Q., Liu, S.-F., Zhuang, Z.-M., Sun, Z.-Z., Liu, C.-L., and Tang, Q.-S. (2012b). The Co-existence of Two Growth Hormone Receptors and Their Differential Expression Profiles between Female and Male Tongue Sole (Cynoglossus Semilaevis). Gene 511, 341–352. doi:10.1016/j.gene.2012.09.082

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Q., Liu, S., Zhuang, Z., Lin, L., Sun, Z., Liu, C., et al. (2012a). Genomic Structure, Polymorphism and Expression Analysis of the Growth Hormone (GH) Gene in Female and Male Half-Smooth Tongue Sole (Cynoglossus Semilaevis). Gene 493, 92–104. doi:10.1016/j.gene.2011.11.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, S., Meng, Z., Chen, R., and Guan, K.-L. (2019). The Hippo Pathway: Biology and Pathophysiology. Annu. Rev. Biochem. 88, 577–604. doi:10.1146/annurev-biochem-013118-111829

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, O. Y., Michalczyk, Ł., Millard, A. L., Emerson, B. C., and Gage, M. J. G. (2017). Lack of Support for Rensch's Rule in an Intraspecific Test Using Red Flour Beetle (Tribolium castaneum) Populations. Insect Sci. 24, 133–140. doi:10.1111/1744-7917.12272

PubMed Abstract | CrossRef Full Text | Google Scholar

Mathews, K. W., Cavegn, M., and Zwicky, M. (2017). Sexual Dimorphism of Body Size Is Controlled by Dosage of the X-Chromosomal Gene Myc and by the Sex-Determining Gene tra in Drosophila. Genetics 205, 1215–1228. doi:10.1534/genetics.116.192260

PubMed Abstract | CrossRef Full Text | Google Scholar

McLean, C. J., Garwood, R. J., and Brassey, C. A. (2018). Sexual Dimorphism in the Arachnid Orders. 6, e5751. doi:10.7717/peerj.5751

PubMed Abstract | CrossRef Full Text | Google Scholar

Mei, J., and Gui, J.-F. (2015). Genetic Basis and Biotechnological Manipulation of Sexual Dimorphism and Sex Determination in Fish. Sci. China Life Sci. 58, 124–136. doi:10.1007/s11427-014-4797-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Michael, D. R., Banks, S. C., Piggott, M. P., Cunningham, R. B., Crane, M., MacGregor, C., et al. (2014). Geographical Variation in Body Size and Sexual Size Dimorphism in an Australian Lizard, Boulenger's Skink (Morethia Boulengeri). PLoS One 9, e109830. doi:10.1371/journal.pone.0109830

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizuno, T., Murakami, H., Fujii, M., Ishiguro, F., Tanaka, I., Kondo, Y., et al. (2012). YAP Induces Malignant Mesothelioma Cell Proliferation by Upregulating Transcription of Cell Cycle-Promoting Genes. Oncogene 31, 5117–5122. doi:10.1038/onc.2012.5

PubMed Abstract | CrossRef Full Text | Google Scholar

Poon, R. Y. C. (2016). Cell Cycle Control: A System of Interlinking Oscillators. Methods Mol. Biol. 1342, 3–19. doi:10.1007/978-1-4939-2957-3_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Reik, W., Dean, W., and Walter, J. (2001). Epigenetic Reprogramming in Mammalian Development. Science 293, 1089–1093. doi:10.1126/science.1063443

PubMed Abstract | CrossRef Full Text | Google Scholar

Reinecke, M., Schmid, A. C., Heyberger-Meyer, B., Hunziker, E. B., and Zapf, J. (2000). Effect of Growth Hormone and Insulin-like Growth Factor I (IGF-I) on the Expression of IGF-I Messenger Ribonucleic Acid and Peptide in Rat Tibial Growth Plate and Articular Chondrocytesin Vivo1. Endocrinology 141, 2847–2853. doi:10.1210/endo.141.8.7624

PubMed Abstract | CrossRef Full Text | Google Scholar

Rensch, B. (1950). Die Abhangigkeit der relativen sexualdifferenz von der Korpergroße. Bonner Zoologische Beiträge 1, 58–69.

Google Scholar

Rotwein, P. (2020). Regulation of Gene Expression by Growth Hormone. Mol. Cell Endocrinol. 507, 110788. doi:10.1016/j.mce.2020.110788

PubMed Abstract | CrossRef Full Text | Google Scholar

Sablowski, R., and Carnier Dornelas, M. (2014). Interplay between Cell Growth and Cell Cycle in Plants. J. Exp. Bot. 65, 2703–2714. doi:10.1093/jxb/ert354

CrossRef Full Text | Google Scholar

Schafer, K. A. (1998). The Cell Cycle: a Review. Vet. Pathol. 35, 461–478. doi:10.1177/030098589803500601

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, C., Li, Q., Chen, S., Zhang, P., Lian, J., Hu, Q., et al. (2014). Epigenetic Modification and Inheritance in Sexual Reversal of Fish. Genome Res. 24, 604–615. doi:10.1101/gr.162172.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, L., Lin, Q., and Su, B. (2015). Estrogen Regulation of Microcephaly Genes and Evolution of Brain Sexual Dimorphism in Primates. BMC Evol. Biol. 15, 127. doi:10.1186/s12862-015-0398-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Si, Y., He, F., Wen, H., Li, J., Zhao, J., Ren, Y., et al. (2016). Genetic Polymorphisms and DNA Methylation in Exon 1 CpG-Rich Regions of PACAP Gene and its Effect on mRNA Expression and Growth Traits in Half Smooth Tongue Sole (Cynoglossus Semilaevis). Fish. Physiol. Biochem. 42, 407–421. doi:10.1007/s10695-015-0147-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Stillwell, R. C., and Fox, C. W. (2007). Environmental Effects on Sexual Size Dimorphism of a Seed-Feeding Beetle. Oecologia 153, 273–280. doi:10.1007/s00442-007-0724-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Székely, T., Lislevand, T., and Figuerola, J. (2007). “Sexual Size Dimorphism in Birds,” in Sex, Size and Gender Roles: Evolutionary Studies of Sexual Size Dimorphism. Editors D. J. Fairbairn, W. U. Blanckenhorn, and T. Székely (Oxford University Press), 27–37. doi:10.1093/acprof:oso/9780199208784.003.0004

CrossRef Full Text | Google Scholar

Tapon, N., Harvey, K. F., Bell, D. W., Wahrer, D. C. R., Schiripo, T. A., Haber, D. A., et al. (2002). Salvador Promotes Both Cell Cycle Exit and Apoptosis in Drosophila and Is Mutated in Human Cancer Cell Lines. Cell 110, 467–478. doi:10.1016/s0092-8674(02)00824-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Teixeira, L. K., and Reed, S. I. (2013). Ubiquitin Ligases and Cell Cycle Control. Annu. Rev. Biochem. 82, 387–414. doi:10.1146/annurev-biochem-060410-105307

PubMed Abstract | CrossRef Full Text | Google Scholar

Thompson, B. J. (2020). YAP/TAZ: Drivers of Tumor Growth, Metastasis, and Resistance to Therapy. Bioessays 42, 1900162. doi:10.1002/bies.201900162

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, N., Wang, R., Wang, R., and Chen, S. (2018). Transcriptomics Analysis Revealing Candidate Networks and Genes for the Body Size Sexual Dimorphism of Chinese Tongue Sole (Cynoglossus Semilaevis). Funct. Integr. Genomics 18, 327–339. doi:10.1007/s10142-018-0595-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P., Zheng, M., Liu, J., Liu, Y., Lu, J., and Sun, X. (2016). Sexually Dimorphic Gene Expression Associated with Growth and Reproduction of Tongue Sole (Cynoglossus Semilaevis) Revealed by Brain Transcriptome Analysis. Ijms 17, 1402. doi:10.3390/ijms17091402

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Yi, Q., Ma, L., Zhou, X., Zhao, H., Wang, X., et al. (2014). Sequencing and Characterization of the Transcriptome of Half-Smooth Tongue Sole (Cynoglossus Semilaevis). BMC Genomics 15, 470. doi:10.1186/1471-2164-15-470

PubMed Abstract | CrossRef Full Text | Google Scholar

Webb, T. J., and Freckleton, R. P. (2007). Only Half Right: Species with Female-Biased Sexual Size Dimorphism Consistently Break Rensch's Rule. PLoS One 2, e897. doi:10.1371/journal.pone.0000897

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, S., Liu, Y., Zheng, Y., Dong, J., and Pan, D. (2008). The TEAD/TEF Family Protein Scalloped Mediates Transcriptional Output of the Hippo Growth-Regulatory Pathway. Developmental Cel. 14, 388–398. doi:10.1016/j.devcel.2008.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, K., Xu, C., Zhang, M., Wang, M., Min, L., Qian, C., et al. (2019). Yes-associated Protein Regulates Podocyte Cell Cycle Re-entry and Dedifferentiation in Adriamycin-Induced Nephropathy. Cell Death Dis. 10, 915. doi:10.1038/s41419-019-2139-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiu, Y., Shao, C., Zhu, Y., Li, Y., Gan, T., Xu, W., et al. (2019). Differences in DNA Methylation between Disease-Resistant and Disease-Susceptible Chinese Tongue Sole (Cynoglossus Semilaevis) Families. Front. Genet. 10, 847. doi:10.3389/fgene.2019.00847

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, T., Wang, W., Zhang, S., Stewart, R. A., and Yu, W. (1995). Identifying Tumor Suppressors in Genetic Mosaics: the Drosophila Lats Gene Encodes a Putative Protein Kinase. Development 121, 1053–1063. doi:10.1242/dev.121.4.1053

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, Y., Zhang, Y., Wang, B., Liu, X., Liu, Q., Song, X., et al. (2018). Leptin and Leptin Receptor Genes in Tongue Sole (Cynoglossus Semilaevis): Molecular Cloning, Tissue Distribution and Differential Regulation of These Genes by Sex Steroids. Comp. Biochem. Physiol. A: Mol. Integr. Physiol. 224, 11–22. doi:10.1016/j.cbpa.2018.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, J. L., Si, Y. F., He, F., Wen, H. S., Li, J. F., Ren, Y. Y., et al. (2015). Polymorphisms and DNA Methylation Level in the CpG Site of the GHR1 Gene Associated with mRNA Expression, Growth Traits and Hormone Level of Half-Smooth Tongue Sole (Cynoglossus Semilaevis). Fish. Physiol. Biochem. 41, 853–865. doi:10.1007/s10695-015-0052-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, H., Xiao, J., Chen, W., Zhou, Y., Tang, Z., Guo, Z., et al. (2014). DNA Methylation of Pituitary Growth Hormone Is Involved in Male Growth Superiority of Nile tilapia ( Oreochromis niloticus ). Comp. Biochem. Physiol. Part B: Biochem. Mol. Biol. 171, 42–48. doi:10.1016/j.cbpb.2014.03.006

CrossRef Full Text | Google Scholar

Zhou, L., Yang, A., Liu, X., Du, W., and Zhuang, Z. (2005). The Karyotype of the Tonguefish Cynoglossus Semilaevis. J. Fish. China 29, 417–419. CNKI:SUN:SCKX.0.2005-03-024.

Google Scholar

Zilberman, D., Coleman-Derr, D., Ballinger, T., and Henikoff, S. (2008). Histone H2A.Z and DNA Methylation Are Mutually Antagonistic Chromatin marks. Nature 456, 125–129. doi:10.1038/nature07324

PubMed Abstract | CrossRef Full Text | Google Scholar