Integrative methylome and transcriptome analysis of Japanese flounder (Paralichthys olivaceus) skeletal muscle during development

DNA methylation is an important epigenetic modification in vertebrate and is essential for epigenetic gene regulation in skeletal muscle development. We showed the genome-wide DNA methylation profile in skeletal muscle tissue of larval 7dph (JP1), juvenile 90dph (JP2), adult female 24 months (JP3) and adult male 24 months (JP4) Japanese flounder. The distribution and levels of methylated DNA within genomic features (1stexons, gene body, introns, TSS200, TSS1500 and intergenic) show different developmental landscapes. We also successfully identified differentially methylated regions (DMRs) and different methylated genes (DMGs) through a comparative analysis, indicating that DMR in gene body, intron and intergenic regions were more compared to other regions of all DNA elements. A gene ontology analysis indicated that the DMGs were mainly related to regulation of skeletal muscle fiber development process, Axon guidance, Adherens junction, and some ATPase activity. Methylome and transcriptome clearly revealed a exhibit a negative correlation. And integration analysis revealed a total of 425, 398 and 429 negatively correlated genes with methylation in the JP2_VS_JP1, JP3_VS_JP1 and JP4_VS_JP1 comparison groups, respectively. And these genes were functionally associated with pathways including Adherens junction, Axon guidance, Focal adhesion, cell junctions, Actin cytoskeleton and Wnt signaling pathways. In addition, we validated the MethylRAD results by bisulfite sequencing PCR (BSP) in some of the differentially methylated skeletal muscle growth-related genes (Myod1, Six1 and Ctnnb1). In this study, we have generated the genome-wide profile of methylome and transcriptome in Japanese flounder for the first time, and our results bring new insights into the epigenetic regulation of developmental processes in Japanese flounder. This study contributes to the knowledge on epigenetics in vertebrates. Author summary Epigenetic mechanisms like DNA methylation have recently reported as vital regulators of some species skeletal muscle development through the control of genes related to growth. To date, although genome-wide DNA methylation profiles of many organisms have been reported and the Japanese flounder reference genome and whole transcriptome data are publically available, the methylation pattern of Japanese flounder skeletal muscle tissue remains minimally studied and the global DNA methylation data are yet to be known. Here we investigated the genome-wide DNA methylation patterns in Japanese flounder, throughout its development. These findings help to enrich research in molecular and developmental biology in vertebrates.


Introduction
As to genetic regulations, a growing number of studies have reported that epigenetic modifications play critical roles in gene expression. Epigenetics refers to heritable changes that modify DNA or associated proteins but without changing the fundamental DNA sequence itself [1]. The epigenome is a dynamic entity influenced by predetermined genetic programs and external environmental cues [2]. DNA methylation is an imporant epigenetic modification of the genome found in most eukaryotes and plays a key role in muscle development. It occurs at the C5 position of cytosine within CpG and non-CpG in the genome. The regulation and mechanisms of the DNA methylation still remain enigmatic, although it is essential for normal development and crucial in many biological processes, such as gene expression regulation,embryogenesis, cellular differentiation, genomic imprinting, Xchromosome inactivation, maintenance of genomic stability by transposon silencing [3][4][5][6]. Now, DNA methylation has attracted much attention owing to its broad impact, reversibility, heritability and genetic characteristics. The profile of DNA methylation across the genome is important to understand DNA methylation dynamics during different developmental muscle development. The genome-wide DNA methylation profiles and functional analysis of many organisms, such as human [7], rat [8], Arabidopsis [9] has been reported. However, little is known about the DNA methylation patterns in Japanese flounder.
Japanese flounder is one of the commercially important marine fish in China and has been widely cultured in recent years. Skeletal muscle represents the most abundant tissue in the body and its features have a direct impact on meat quality.
Understanding the growth and development of skeletal muscle is important. Skeletal muscle development is a very complicated but precisely regulated process, which contains four steps: determination of myoblasts, proliferation of myoblasts, differentiation and fusion of myoblasts into myotubes and myofibers, and growth and maturation until postnatal [10][11]. The research of complex mechanism underlying skeletal muscle development is helpful to genetic improvement for meat quality. In transcriptional silencing and suppresses the corresponding protein products in most eukaryotes [6,[12][13][14]. Thus, DNA methylation plays critical roles in cellular processes and the development of skeletal muscle tissue.
There are many approaches to decipher a genome-wide DNA methylation profile, including Methylation-dependent restriction-site associated DNA sequencing (MethyIRAD), MeDIP-seq and whole-genome bisulfite sequencing (WGBS). The gold standard to determine the DNA methylome is genome-wide bisulfite sequencing, which firstly converts all the unmethylated cytosines into uracil while left the methylated cytosines unchanged by sodium bisulfite under denaturing conditions, which can be distinguished subsequently by sequencing [15]. whereas genome-wide bisulfite sequencing is highly expensive and time-consuming. Now, studies have shown that MethylRAD is a suitable method for high-throughput sequencing to analyze the DNA methylation status of methylated genome regions at a fraction of the cost and time of genome-wide bisulfite sequencing. MethylRAD uses methylation-dependent restriction enzymes which can specifically discriminate methylated cytosines between CG and non-CG methylation. These enzymes have the unique ability to produce 32-base-long fragments around fully methylated restriction sites, which are suitable for high-throughput sequencing to profile cytosine methylation on a genomic scale [16,17]. Many recent studies have shown that MethyIRAD can reflect the relative genome-wide DNA methylation profile [18,19].
Therefore, we chose MethyIRAD to analyze genome-wide profiles of DNA methylation in Japanese flounder in our study. In this study, we have performed the first integrated genome-wide analysis of DNA methylation, and mRNA transcriptional activity, using the transcriptome and MethylRAD (a simple genomic methylation site detection method) [16,20]

Global mapping of DNA methylation in Japanese flounder
The MethylRAD analysis was used to study the global mapping of DNA methylation pattern in the Japanese flounder skeletal muscle tissues of JP1, JP2, JP3 and JP4. We generated about 439 -751 million raw reads from each samples. After low-quality data filtration, about 42 to 54 million reads assessed as clean data were analyzed and mapped (S1 Table). Base distribution and quality distribution maps of clean reads were plotted (S1 Fig). Of the high-quality methylation tag libraries in the twelve samples, 71.76-81.08% were comparable to unique positions with a highquality read alignment against the Japanese flounder reference genome using SOAP software (version 2.21) (S1 Table).
The percentage of the DNA methylation sites (CCGG sites and CCWGG sites) in We found a substantial amount of CCGG methylation and a small amount of CCWGG methylation. Therefore, we analyzed the genome coverage of the CCGG, CCWGG sites under different sequencing depth (S2 Table).
The sequencing depths of the DNA methylation sites (CCGG sites and CCWGG sites) in each sample are shown in a box plot (the number of methylation sites in each sample had a depth higher than 3) (Fig 1). The JP1, JP2, JP3 and JP4 methylation sites (CCGG/CCWGG) were identified on the chromosomes of Japanese flounder.
Methyl-RAD reads were detected in most chromosomal regions (chromosomes  in each group ( Figure 2).

Distribution of DNA methylation sites of different functional regions
The distribution of MethyIRAD reads that were aligned on a unique locus in different genome regions represents a genome-wide methylation pattern. We obtained the DNA methylation site annotation of the Japanese flounder genome and the comparison of average methylation sites showed that there were differential methylation site distribution in different components of the genome. The distribution patterns of most methylated sites at the different elements of genomes were similar in the four groups. We found that the major proportion of DNA methylation sites were mainly enriched in the intergenic regions followed by the regions at the gene body, of methylation site than those in the JP2 group, while JP3 group showed lower number of methylation site than those in the JP4 group ( Fig 3B). The distribution of methylation sites on different gene elements in each sample indicates that the skeletal muscle growth difference during different developmental stages might be associated with global methylation.

Relative quantification of DNA methylation levels around the Gene body
We found that the DNA methylation site distribution curve had TSS representing an upstream sequence centered on the transcription initiation site, and TTS representing a downstream sequence centered on the transcription termination site. The DNA methylation levels at either the CCGG sites or the CCWGG sites in the gene regions were similar for four groups (Fig 4).

Differentially methylated regions (DMRs) analysis
To characterize the differences of DNA methylation levels among samples, DMRs were detected. For assessing the methylation level of differential methylation sites between four groups for the three biological replicates, the cluster heat map was shown to further show the changes in CCGG/CCWGG methylation levels among the The number of hypermethylation DMRs is less than hypomethylation DMRs.
The number of DMRs in CCGG sites is lower than that in CCWGG sites. DMRs that are unique or shared among the four groups are shown ( Fig 5). The results of a boxplot analysis of DMRs showed that the methylation level of the JP2 group was the lowest among four groups and the JP3 group is lower than that in JP4 group (S4 Fig).
The pie map distribution of differential methylation sites with differential methylation

Analysis of differential methylation site-related genes
To investigate the differential methylation site-related genes regulatory role, the function of a gene was described by the GO and KEGG enrichment analysis of the gene where the differential DNA methylation sites were located. We found that these The bisulfite sequencing results showed a high degree of consistency with the MethylRAD data (Fig 6). These results indicated that our genome-wide methylation results obtained by MethylRAD are reliable.

Transcriptome assembly and annotation
Using RNA-Seq, this study compared the transcriptomic landscapes of skeletal muscle from the larval 、 juvenile and adult ( female and male ) stages used to construct mRNA libraries. All the samples sequenced on the Illumina HiSeq X Ten  Table ).

Comparative and enrichment analysis of differentially expressed genes.
We defined genes with fold changes > 2 and P-values < 0.05 were recognized as significantly differentially expressed. Different expression genes that are unique or

RNA-Seq data validation
To examine the reliability of the RNA-seq results, three DEGs (MyoD1, Six1 and Ctnnb1) involved in the development of skeletal muscle were selected for validation using qRT-PCR. The mRNA expression levels of these key genes, such as

Association analysis of methylRAD and the transcriptome (RNA-Seq).
The association analyses between the transcriptome and methylation were based on RNA-Seq and MethylRAD sequencing data. We calculated each gene's methylation level and expression level in larval, juvenile and adult female and adult male period Japanese flounder in terms of DNA methylation and the mRNA transcriptome. We observed that methylation levels correlates negatively with expression levels in both the CCGG and CCWGG pattern in larval, juvenile and adult Japanese flounder skeletal muscle tissue ( Fig 9B). Furthermore, to explore the relationship between these DMGs and the DEGs found at the transcriptome level, an association analysis was performed. We found that a lot of genes that both different  Table). We speculate that DNA methylation mainly affects skeletal muscle development in the larval period, and differences in skeletal muscle between juvenile and adult period Japanese flounder may be affected by other factors rather than DNA methylation, and that require further research.
Among these DEGs which also exhibited differential methylation levels, a large  Table).
To further investigate the signaling pathway associated with negatively correlated genes with methylation and expression levels in juvenile and adult compare to larval period, we performed KEGG enrichment analysis of these genes. The results showed that there were significantly enriched KEGG signaling pathway (P <0.05) between larval and juvenile, between larval and female male adult, and between larval and male adult Japanese flounder, respectively (S10 Fig). In our study, KEGG enrichment analysis screens criteria for pathway entries with the number of differential genes greater than 2, however, there are too few differential genes to show neural axons follow specific, predictable paths to reach their target locations [37].
Differential methylation changes in this pathway were used as a focus to identify how epigenetic changes during aging could potentially associate with the well-known decrease of skeletal muscle function with increasing age [38]. In addition, we also found the signaling networks that guide diverse cell behaviours and functions are connected to tight junctions transmitting information to and from the cytoskeleton on Japanese flounder growth still require further study in the future. This study expands the Japanese flounder methylated genes and could initiate further study in the muscle development of Japanese flounder.
In addition, We discovered some differentially DNA methylated genes involved in skeletal muscle development in larval、juvenile and adult Japanese flounder. For example, we found that MyoD1, a master regulatory gene of skeletal muscle differentiation [46]. Another well-known gene named myf6 (MRF4) is involved in inducing fibroblasts to differentiate into myoblasts and affects skeletal muscle development [47]. Six1 has been shown to play a pivotal role in skeletal muscle development [48][49][50] which is a transcription factor essential for embryonic myogenesis and also regulates MyoD1 expression in muscle progenitor cells. In addition, a recent study shows that Six1 contributes to the regeneration of adult muscle by enhancing and maintaining MyoD1 expression in adult muscle satellite cells in addition to its role in embryonic muscle formation [51]. MyoD1 is able to promote the transformation of multipotent stem cells to skeletal muscle by binding and activating the expression of a subset of pre-myogenic mesoderm genes, including Six1 [52]. And gene Ctnnb1 modulates skeletal muscle development by acting on transcription factors controlling myogenesis such as MyoD [53].
We believed that the methylation of these genes might partially contribute to the growth, which contributes to muscle growth-related genes, is still unclear and require further study in the future.

Conclusions
We

Ethics statement
All experimental procedures and sample collection were conducted according to the guidelines and were approved and supervised by the respective Animal Research and Ethics Committees of Ocean University of China. The field studies did not involve endangered or protected species. The fish were all euthanized by tricaine methanesulfonate (MS-222).

Experimental fish and data collection
The experimental animals were collected from Donggang District Institute of marine treasures in Rizhao of Shandong province, and were temporary reared in a 500L bucket in seawater in Ocean University of China within the same environment. in the CmCGG and mCCWGG sites, and generate a double-stranded DNA break on the 30 side of the modified cytosine at a fixed distance (N12/N16).Accordingly, symmetrical DNA methylated sites were bidirectionally cleaved by FspEI to generate 32-base-long fragments. Then, two adaptors were added to the digested DNAs by T4 DNA ligase (NEB, USA), and the ligation products were amplified in 20 μl reactions by specific primers. PCR products were purified using a MinElute PCR Purification Kit (Qiagen) and pooled for sequencing using the Illumina X-ten PE 150 sequencing platform [16]. Base quality values were calculated using a Phred quality score (Q sanger=-10log10p). Input sequencing data before operation and computing were called raw reads. Raw reads were first subject to quality filtering and adaptor trimming. After operation, the data, including adapter reads and low-quality sequences, were removed from raw reads as clean reads.

DNA methylation data analysis
To improve the accuracy in the following analysis,filtering pair-end sequencing paired clean reads according to the following terms: (i) remove low quality reads significant. The function of the gene was described by a GO and KEGG function enrichment analysis of the gene where the differential methylation site was located.
The number of genes included in each GO entry and KEGG pathway was counted and the significance of gene enrichment for each GO entry and KEGG pathway was calculated using the hypergeometric distribution test [57], GO entries with the number of corresponding genes greater than 2 in three categories were screened and the GO enrichment analysis results. Differences were considered significant at P < 0.05.

RNA library construction and high-throughput sequencing
Total Clean reads were generated by filtering the low-quality reads, reads containing adapter, reads containing ploy-N from the raw reads of fastq form by NGS QC Toolkit [58]. All clean reads with high quality were annotated and classified by mapping them to the Japanese flounder reference genome by Tophat (http://tophat.cbcb.umd.edu/).
The expressed genes were confirmed based on the annotation information of the clean reads. The expression level of each gene was calculated and normalized by the fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) [59] using bowtie2 [60] and eXpress (v1.5.1) software [61].
Differential expression analysis of the genes was performed by using the DESeq R package (2012). The NB (negative binomial distribution test) was used to test the difference in the number of reads. The transcript expression was estimated by the basemean value. The significantly DEGs between the two arbitrary samples were identified based on the following thresholds: fold changes > 2 and P-values < 0.05.
The assembled transcripts were annotated by Genomes (KEGG) and Gene Ontology (GO). The relevant biological process、cellular component and molecular function of the GO categories and KEGG biological pathways were identified through gene enrichment analyses [62]. The hypergeometric test was conducted to identify the significantly enriched GO terms and KEGG pathway (corrected p-value < 0.05).

Quantitative RT-PCR
The differential expression patterns of the genes detected by transcriptome data were validated by qRT-PCR analysis. The specific primer pairs were designed for the detection of corresponding genes (S6 Table). The 18S gene from Japanese flounder recommendations [63]. Three sample from each developmental stage were used, and three technological replicates were performed to ensure the reliability of quantitative analysis.