DNA Methylation Is Correlated with Gene Expression during Diapause Termination of Early Embryonic Development in the Silkworm (Bombyx mori)

DNA modification is a naturally occurring DNA modification in prokaryotic and eukaryotic organisms and is involved in several biological processes. Although genome-wide methylation has been studied in many insects, the understanding of global and genomic DNA methylation during insect early embryonic development, is lacking especially for insect diapause. In this study, we analyzed the relationship between DNA methylomes and transcriptomes in diapause-destined eggs compared to diapause-terminated eggs in the silkworm, Bombyx mori (B. mori). The results revealed that methylation was sparse in this species, as previously reported. Moreover, methylation levels in diapause-terminated eggs (HCl-treated) were 0.05% higher than in non-treated eggs, mainly due to the contribution of CG methylation sites. Methylation tends to occur in the coding sequences and promoter regions, especially at transcription initiation sites and short interspersed elements. Additionally, 364 methylome- and transcriptome-associated genes were identified, which showed significant differences in methylation and expression levels in diapause-destined eggs when compared with diapause-terminated eggs, and 74% of methylome and transcriptome associated genes showed both hypermethylation and elevated expression. Most importantly, Kyoto Encyclopaedia of Genes and Genomes (KEGG) analyses showed that methylation may be positively associated with Bombyx mori embryonic development, by regulating cell differentiation, metabolism, apoptosis pathways and phosphorylation. Through analyzing the G2/M phase-specific E3 ubiquitin-protein ligase (G2E3), we speculate that methylation may affect embryo diapause by regulating the cell cycle in Bombyx mori. These findings will help unravel potential linkages between DNA methylation and gene expression during early insect embryonic development and insect diapause.


Introduction
DNA methylation is an important epigenetic modification in plants and animals [1][2][3][4]. It has been shown to play key roles in a broad range of biological processes, including embryonic development [5], 3 of 22 The silkworm is an ideal model for studying relationships between DNA methylation, embryonic development and insect diapause.
Currently, functional analyses of DNA methylation during insect embryogenesis, especially insect diapause, are limited. Our previous study provided functional insights into BmDnmt1 and BmDnmt2 in the regulation of silkworm embryonic development; we showed that the expression of BmDnmt1 and BmDnmt2 was elevated in diapause-terminated eggs that were HCl-treated when compared with diapause-destined eggs [40].
In this study, we investigate DNA methylation patterns using WGBS and RNA sequencing (RNA-seq) technologies in diapause-terminated eggs compared with diapause-destined eggs in B. mori. We explore different methylation patterns and methylation modification genes and pathways associated with embryonic development. Based on our findings, DNA methylation appears to be essential during B. mori diapause and embryogenesis.

Transcriptional Profiles of Silkworm Diapause-Destined Eggs and Diapause-Terminated Eggs
To explore gene expression mechanisms during early embryonic development in the silkworm, we compared the transcriptional profiles of diapause-terminated eggs treated with HCl (HCl-treated) and diapause-destined eggs (non-HCl-treated) using RNA-seq, and the raw sequencing data are available from the repository of NCBI Sequence Read Archive (SRA) with the accession No. PRJNA598980. More than 48 million clean reads were obtained from each sample. The CG content of each of the six libraries was approximately 44%, and CycleQ30 was >94% for each library. The proportion of total reads that mapped to the reference genome reached approximately 95% coverage. In addition, 1732 differentially expressed genes (DEGs) were identified, with 856 genes upregulated and 876 genes downregulated ( Table 1). The details of DEGs are shown in supplementary  Table S2. Gene Ontology (GO) analysis showed that common DEGs were mainly distributed to binding, transporter activity, locomotion, catalytic activity ( Figure 1A). Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analyses showed that DEGs were involved in cellular processes, environmental information processing, genetic information processing and metabolism and organismal systems pathways ( Figure 1B). It is noteworthy that the most upregulated genes were enriched for signal transduction and carbohydrate metabolism pathways, with 39 and 30 DEGs, respectively. In addition, the most downregulated genes were enriched for lipid metabolism and amino acid metabolism pathways, with 24 and 19 DEGs, respectively. Table 1. Statistical analysis of RNA-seq data obtained of diapause-destined eggs (non-treated) and diapause-terminated eggs (HCl-treated).  The expression of BmDnmt1 (and BmDnmt2 to a lesser extent) was upregulated in diapause-terminated eggs after HCl-treatment. This was consistent with our previous findings [40] and consistent with RNA-seq data generated in this study ( Figure 2). These changes in DNA The expression of BmDnmt1 (and BmDnmt2 to a lesser extent) was upregulated in diapause-terminated eggs after HCl-treatment. This was consistent with our previous findings [40] and consistent with RNA-seq data generated in this study ( Figure 2). These changes in DNA methyltransferase expression suggest that genomic DNA methylation patterns may be associated with diapause initiation and termination in silkworm embryonic development.

Overview of Methylation Landscapes of Silkworm Diapause-Destined Eggs and Diapause-Terminated Eggs
To investigate methylation patterns in silkworm embryonic development, WGBS was conducted in diapause termination eggs treated with HCl (HCl-treated) and diapause-destined eggs (non-HCl-treated). Three repeats were performed for each of the two groups, and the raw sequencing data are available from the repository of NCBI Sequence Read Archive (SRA) with the accession No. PRJNA598995. Up to ~9000 million sequencing reads were generated for each replicate, yielding ~12 Gb of data representing >30× the NCBI reference genome for silkworm. Nearly 98% of the reads (bases) were retained and were high-quality clean reads (bases) ( Table 2). The non-HCl sample (control) genome presented ~0.21% mC on the total sequenced C sites, which reflected the methylation level percentage in the genome. It contained ~98.55% (mCG), 0.27% (mCHG) and 1.18% (mCHH) on the total sequenced mC sites. Accordingly, the HCl-treated sample presented ~0.26% of the total sequenced C sites, and 98.94%, 0.21% and 0.85% of the CG, CHG and CHH sites, respectively (Supplementary Table S3). It was observed that the mC rate of the total sequenced C sites was significantly increased in diapause-terminated eggs when compared with diapause-destined eggs ( Figure 3A). In addition, our results showed that mCG sites were significantly increased in diapause-terminated eggs when compared with diapause-destined eggs, while others were not significantly changed ( Figure 3B).

Figure 2.
Analysis of methyltransferase expression. The expression fold change of BmDnmt1 and BmDnmt2 by qRT-PCR and RNA-seq, respectively, comparing diapause-terminated eggs with diapause-destined eggs.

Overview of Methylation Landscapes of Silkworm Diapause-Destined Eggs and Diapause-Terminated Eggs
To investigate methylation patterns in silkworm embryonic development, WGBS was conducted in diapause termination eggs treated with HCl (HCl-treated) and diapause-destined eggs (non-HCl-treated). Three repeats were performed for each of the two groups, and the raw sequencing data are available from the repository of NCBI Sequence Read Archive (SRA) with the accession No. PRJNA598995. Up to~9000 million sequencing reads were generated for each replicate, yielding~12 Gb of data representing >30× the NCBI reference genome for silkworm. Nearly 98% of the reads (bases) were retained and were high-quality clean reads (bases) ( Table 2). The non-HCl sample (control) genome presented~0.21% mC on the total sequenced C sites, which reflected the methylation level percentage in the genome. It contained~98.55% (mCG), 0.27% (mCHG) and 1.18% (mCHH) on the total sequenced mC sites. Accordingly, the HCl-treated sample presented~0.26% of the total sequenced C sites, and 98.94%, 0.21% and 0.85% of the CG, CHG and CHH sites, respectively (Supplementary  Table S3). It was observed that the mC rate of the total sequenced C sites was significantly increased in diapause-terminated eggs when compared with diapause-destined eggs ( Figure 3A). In addition, our results showed that mCG sites were significantly increased in diapause-terminated eggs when compared with diapause-destined eggs, while others were not significantly changed ( Figure 3B).

Validation of RNA-seq and WGBS
For the validation of RNA-seq results, three genes that were significantly upregulated and three genes that were significantly downregulated were selected for mRNA quantification using RT-qPCR. The results were consistent with RNA-seq ( Figure 4A). McrBC is a restriction enzyme that cuts methylated, but not unmethylated, DNA. McrBC-qPCR has been routinely used for DNA methylation testing and marker validation [41,42]. For WGBS validation, we assayed DNA methylation levels of two genomic regions using methylation-sensitive McrBC-qPCR. The gene LOC101739208 is hypermethylated, and we determined that its methylation levels increased by 28% in HCl-treated samples when compared with non-treated samples. In another experiment using the gene LOC110385781, which is hypomethylated, we observed that methylation levels decreased by 25% in HCl-treated samples when compared with non-treated samples. Methylated DNA can be digested by McrBC; thus, higher qPCR signals indicate lower methylation levels, suggesting our McrBC-qPCR assay was consistent with our WGBS data ( Figure 4B).

Validation of RNA-seq and WGBS
For the validation of RNA-seq results, three genes that were significantly upregulated and three genes that were significantly downregulated were selected for mRNA quantification using RT-qPCR. The results were consistent with RNA-seq ( Figure 4A). McrBC is a restriction enzyme that cuts methylated, but not unmethylated, DNA. McrBC-qPCR has been routinely used for DNA methylation testing and marker validation [41,42]. For WGBS validation, we assayed DNA methylation levels of two genomic regions using methylation-sensitive McrBC-qPCR. The gene LOC101739208 is hypermethylated, and we determined that its methylation levels increased by 28% in HCl-treated samples when compared with non-treated samples. In another experiment using the gene LOC110385781, which is hypomethylated, we observed that methylation levels decreased by 25% in HCl-treated samples when compared with non-treated samples. Methylated DNA can be digested by McrBC; thus, higher qPCR signals indicate lower methylation levels, suggesting our McrBC-qPCR assay was consistent with our WGBS data ( Figure 4B).

DNA Methylation Patterns in Genomic Regions of Diapause-Destined Eggs and Diapause-Terminated Eggs
From our data, we observed that DNA methylation patterns varied at different genomic intervals. To understand these patterns in silkworms, we analyzed methylation profiles within genes, including Coding sequence (CDS), downstream, exons, intergenic areas, introns, lnc-RNAs, miRNAs, mRNAs, transcripts and upstream ( Figure 5A). When compared with the non-treated group, CG methylation levels increased in the HCl-treated group for all types of genomic intervals, with the highest methylation levels found in transcript elements. Methylation profiles were also analyzed in transposon elements (TEs) and repeat elements (REs) using RepeatMasker software (GIRI, v4.0.4). Transposons include DNA segments, long interspersed nuclear elements (LINE), short interspersed nuclear elements (SINE), retrotransposons (RC) and long terminal repeats (LTRs), and repeat elements include simple_repeat, low_complexity, satellites and RNArepeats ( Figure 5B); in particular, satellites of HCl treated were not detected. The results suggest that CG context methylation sites are mainly concentrated in SINE repeats, DNA transposons and RC transposons.
To understand the relationship between DNA methylation profiles and gene expression, DNA methylation profiles were divided into~2 kb upstream and downstream transcribed regions to study methylation changes ( Figure 6). The distribution of CG methylation sites in three distinct genomic functional regions showed strong bias, and CG methylation showed high enrichment levels in~2 kb upstream, transcription initiation sites (TSS) downstream and transcription termination sites (TTS) downstream regions. In addition, CG methylation levels in HCl-treated samples were elevated among gene features when compared with non-treated samples, except for regions near TSS and TTS. . Three genes that were significantly upregulated (101742044, 101739131 and 101739208) and three genes that were significantly downregulated (110385781, 101740063 and 101741941) from RNA-seq analyses were selected for mRNA quantification using RT-qPCR (A). DNA methylation levels of the hypermethylated gene, (101739208) and the hypomethylated gene, (110385781) from our WGBS were analyzed using methylation-sensitive McrBC-qPCR (B). Methylated DNA was digested by McrBC; thus, higher qPCR signals indicate lower methylation levels, and lower qPCR signals indicate higher methylation levels. Significant differences are indicated by ** p < 0.01 and * p < 0.05.

DNA Methylation Patterns in Genomic Regions of Diapause-Destined Eggs and Diapause-Terminated Eggs
From our data, we observed that DNA methylation patterns varied at different genomic intervals. To understand these patterns in silkworms, we analyzed methylation profiles within genes, including Coding sequence (CDS), downstream, exons, intergenic areas, introns, lnc-RNAs, miRNAs, mRNAs, transcripts and upstream ( Figure 5A). When compared with the non-treated group, CG methylation levels increased in the HCl-treated group for all types of genomic intervals, with the highest methylation levels found in transcript elements. Methylation profiles were also analyzed in transposon elements (TEs) and repeat elements (REs) using RepeatMasker software (GIRI, v4.0.4). Transposons include DNA segments, long interspersed nuclear elements (LINE), short interspersed nuclear elements (SINE), retrotransposons (RC) and long terminal repeats (LTRs), and repeat elements include simple_repeat, low_complexity, satellites and RNArepeats ( Figure 5B); in particular, satellites of HCl treated were not detected. The results suggest that CG context methylation sites are mainly concentrated in SINE repeats, DNA transposons and RC transposons.
To understand the relationship between DNA methylation profiles and gene expression, DNA methylation profiles were divided into ~2 kb upstream and downstream transcribed regions to study methylation changes ( Figure 6). The distribution of CG methylation sites in three distinct genomic functional regions showed strong bias, and CG methylation showed high enrichment levels in ~2 kb upstream, transcription initiation sites (TSS) downstream and transcription termination sites (TTS) downstream regions. In addition, CG methylation levels in HCl-treated   different genetic elements, including CDS, downstream, exon, gene, intergenic, intron, inc-RNA, miRNA, mRNA, transcript, and upstream (A). Methylation percentages of transposon elements (TEs) and repeat elements (REs) (B). DNA, long interspersed nuclear elements (LINE), short interspersed nuclear elements (SINE), retrotransposons (RC) and long terminal repeats (LTR) are classified as TE, and Simple_repeat, Low_complexity, satellite and RNA are classified as RE.

Differentially Methylated Regions (DMRs) and Differentially Methylated Promoters (DMPs) Responding to Embryonic Development
To investigate the influence of embryonic development and insect diapause on methylation, we analyzed DMRs in diapause-terminated eggs when compared with diapause-destined eggs. In total, 43,781 DMRs were identified (methylation difference ≥ 10%, Q value ≤ 0.05) containing 42,877 hypermethylated DMRs and 904 hypomethylated DMRs. In addition, 4449 different methylation genes (DMGs) located in DMRs were identified. To investigate the biological functions of DMGs, a KEGG pathway enrichment analysis was performed. As shown in Figure 7A, DMGs were mainly assigned to pathways related to protein processing in endoplasmic, cell cycle, RNA transport, etc. functions.
The promoter region of a gene is closely related to transcriptional regulation. Differentially methylated modifications located in the promoter region may change a gene's transcriptional expression. To characterize such methylation alterations in promoter regions during early embryonic development, 99 differentially methylated promoters (DMPs) located in DMRs were identified (Supplementary Table S4). KEGG pathway analyses showed that DMPs were involved in hippo-signalling, Wnt-signalling, mTOR-signalling, etc. pathways. (Figure 7B). assigned to pathways related to protein processing in endoplasmic, cell cycle, RNA transport, etc. functions.
The promoter region of a gene is closely related to transcriptional regulation. Differentially methylated modifications located in the promoter region may change a gene's transcriptional expression. To characterize such methylation alterations in promoter regions during early embryonic development, 99 differentially methylated promoters (DMPs) located in DMRs were identified (Supplementary Table S4). KEGG pathway analyses showed that DMPs were involved in hippo-signalling, Wnt-signalling, mTOR-signalling, etc. pathways. ( Figure 7B).

Correlations Between DNA Methylation and Gene Expression during Silkworm Early Embryonic Development
To investigate whether changes in CG methylation observed during silkworm early embryonic development altered gene expression, we synthetically analyzed DEG and DMG data. As a result, there were 364 methylome-and transcriptome-associated genes (MTGs) with significant differences

Correlations Between DNA Methylation and Gene Expression during Silkworm Early Embryonic Development
To investigate whether changes in CG methylation observed during silkworm early embryonic development altered gene expression, we synthetically analyzed DEG and DMG data. As a result, there were 364 methylome-and transcriptome-associated genes (MTGs) with significant differences in methylation and expression levels in diapause-terminated eggs when compared with diapause-destined eggs (Supplementary Table S5). More than 98% of genes showed hypermethylation, and more than 74% of genes showed both hypermethylation and upregulated protein expression. To be specific, 270 hypermethylated genes were positively correlated with expression changes in the first quartile, and 88 hypermethylated genes were negatively correlated with expression changes in the second quartile ( Figure 8A), indicating a strong positive correlation between DNA methylation and gene expression for first quadrant genes. Meanwhile, percentage methylation levels in gene body, 2 kb upstream and downstream regions for positively and negatively correlated genes, were detected, to understand the relationship between DNA methylation profiles and gene expression (Supplementary Figure S1). The results showed that CG methylation sites of 270 hypermethylated-upregulated genes showed strong bias. The CG methylation levels in HCl-treated samples were strongly elevated in 2 kb downstream and weakly elevated near the TSS of gene body when compared with non-treated samples. By contrast, there is no obvious change trend for the 88 hypermethylated-downregulated genes. Moreover, GO analysis showed that MTGs were highly enriched for phosphatidylinositol dephosphorylation, G1/S transition in the mitotic cell cycle, RNA polymerase II transcription corepressor activities etc., (Figure 8B), and related genes with changes in methylation levels or methylation patterns may be closely correlated with embryonic development. KEGG pathway analyses showed that MTGs were involved in signalling pathways regulating stem cell pluripotency, phosphatidylinositol signalling system, different types of N−glycan biosynthesis, insulin-signalling pathways etc. (Figure 8C), suggesting that DNA methylation variations appear to regulate gene expression in multiple pathways during early embryonic development.

MTGs Are Involved in Metabolic Biosynthesis
In embryonic development, when the diapause eggs have terminated the diapause, energy metabolism and material synthesis begin to be activated. Here, KEGG pathway analyses showed that 23 MTGs related to metabolic biosynthesis were identified for pathways involved in N-glycan biosynthesis, glycosaminoglycan biosynthesis, glutathione metabolism, regulation of lipolysis in adipocytes, biosynthesis of amino acids and glyoxylate and dicarboxylate metabolism (Table 3). These pathways were associated with metabolism and synthesis of glycogenolysis, protein synthesis and lipolysis. Moreover, the most gene-enriched pathway was the protein-synthesis-related pathway, in which 16 MTGs were identified. Previous studies have shown that glycogen synthesis and tricarboxylic acid cycle dynamics are enhanced during embryonic development [43]. Here, four MTGs were involved in glycogen synthesis and three MTGs in the tricarboxylic acid cycle. Three MTGs involved in lipid metabolism were also identified; it is worth noting that a homologue of adenylate cyclase type 2 (Gene ID: LOC101737606) had a sevenfold increase in protein expression. Thus, KEGG analyses demonstrated that methylation modification was related to embryonic metabolic biosynthesis.

MTGs Are Involved in Phosphorylation
A previous study has shown that some key enzymes in silkworm embryo development regulate their activities through protein phosphorylation and dephosphorylation to sustain embryonic metabolism [44]. Our KEGG pathway analyses showed that six MTGs were involved in phosphatidylinositol-signalling pathways (Table 4), with all genes showing increased methylation and expression levels in diapause-terminated eggs when compared with diapause-destined eggs. It is noteworthy that three INPP-family genes of the six MTGs were identified, and that INPP genes regulate many cellular activities, including vesicle transport, cytoskeletal dynamics, protein synthesis, cell proliferation and survival [45]. Based on these observations, our data indicate that methylation may play a role in embryonic development by regulating phosphorylation.

MTGs are Involved in the Cell Cycle, Apoptosis and Stem Cell Pluripotency
Silkworm embryonic cells are arrested in G2 at diapause, and cell cycles become slower in proportion to an increasing length of G1. Once diapause terminates, the embryos are competent to resume development at 25 • C, cells enter the M phase, and cell division then resumes in the embryos [38]. Based on this evidence, we hypothesized that embryonic development and diapause were related to cell cycle regulation. Here, KEGG pathway analyses showed that 16 MTGs were involved in G1/S transition of the mitotic cell cycle, regulation of apoptotic processes and signalling pathways regulating stem cell pluripotency (Table 5). Moreover, six MTGs were enriched in the G1/S transition of the mitotic cell cycle, and we observed that 101739208, a homologue of the G2/M phase-specific E3 ubiquitin-protein ligase (G2E3), was upregulated and hypermethylated. Alternately, diapause-destined embryos are arrested at gastrulation, which occurs two days after the germ-band formation stage (24 h after oviposition). At this stage, cell differentiation gradually ceases, and the embryo gradually forms [46]. KEGG analyses revealed that six MTGs were enriched in signalling pathways regulating stem cell pluripotency, which increased gene expression. Those MTGs can promote cell differentiation and organ formation in diapause-terminated eggs. Based on this evidence, these data suggest that diapause-destined eggs may regulate the cell cycle, apoptosis and stem cell pluripotency to terminate diapause via methylation.

Discussion
DNA methylation is one of the most widespread epigenetic markers in the genome and has been linked to insect development, specifically embryonic development [16,47]. Our previous study showed that BmDNMTs were highly expressed during embryonic development, especially at early embryonic stages; its expression was significantly upregulated in diapause-terminated eggs after HCl-treatment [40]. However, DNA methylation regulation in silkworm embryonic development and diapause still remains unclear. This study provided insights into the potential cytosine methylation in diapause-destined eggs and diapause-terminated eggs using WGBS in relation to methylation and embryonic development.
Overall, the levels of DNA methylation in most insects is <1% [28,48,49]. This research has shown~0.21%-0.26% mC methylation levels on the total sequenced C sites, and more than~98.5% mC methylation in CG sequences. Based on previous methylome studies in the silkworm, the average methylation levels in the mid-gut and fat body at mC methylation in the genome were approximately 1% [29], and~0.11% for silk glands [28]. When comparing eggs, the mid-gut, fat body and silk gland, it has been shown that DNA methylation levels have tissue specificity in the silkworm. Moreover, mC rates in sequenced C sites were significantly increased in diapause-terminated eggs (HCl-treated) when compared with diapause-destined eggs (non-treated), and change in methylation levels mainly comes from CG sites but not from CHG and CHH sites. It is suggested that when compared with other plants and animals which have higher methylation levels at CHG and CHH sites [50,51], DNA methylation works primarily through CGs site in the silkworm.
DNA methylation generally functions as a repressive transcriptional signal, but it is also known to activate gene expression [52,53]. Here, we calculated methylation levels in the context of gene regions and~2 kb upstream and downstream regions. Consistently, CG sites showed more regularity than CHG and CHH sites. CG methylation, but not CHG or CHH methylation, exhibited a characteristic peak in the body of protein-coding genes and showed significant increases after TSS and TES, which maintain low methylation levels. These results are similar to those previously reported for other species [54]. Moreover, CG methylation levels in the HCl-treated group were generally higher than those in the non-treated group. Thus, we speculate that CG methylation occurs in gene body regions, enhancing gene transcription for embryonic development.
To investigate the specific modification regions of DNA methyl groups in the silkworm genome, methylation profiles of genetic elements were analyzed, including CDS, downstream, upstream, exons, etc. Notably, methylation levels in all genetic elements were increased in HCl-treated samples when compared with non-treated samples, and methylation mainly occurred in coding sequences such as transcript, CDS and exons. Previous studies in silkworms have shown that CG methylation is substantially enriched in gene bodies and is positively correlated with gene expression levels, suggesting positive roles in gene transcription [26,28]. Therefore, hypermethylation of coding sequences may be positively associated with the terminating of diapause in silkworm embryos.
DNA methylation patterns are prevalent in TEs and REs [17]. Our results show that methylation in TEs, especially SINEs, DNA transposable, and RC are higher than in others. SINEs are an abundant class of TEs found in a wide variety of eukaryotes, mainly integrate into hypomethylated DNA regions and are targeted by methylases for de novo methylation [55][56][57]. Previous studies have shown that SINEs promote the amplification of gene enrichment region fragments or stress-induced gene expression, which increase gene expression [58]. Based on these analyses, hypermethylated SINEs may promote gene transcription and expression during embryo development.
Differentiation and organogenesis are two key phases of embryogenesis in the silkworm. Metabolic biosynthesis is also an essential activity during this period. Metabolic pathways were significantly activated after HCl-treatment, especially for glycogenolysis, protein synthesis and lipolysis. It is worth noting that HS6ST2 homologous gene (GeneID: LOC101737390) regulates energy metabolism and promoted embryonic development in Drosophila [59,60]. Moreover, the Pi3k60 homologous gene (GeneID: LOC100158253), which was upregulated and hypermethylated, reportedly promotes the decomposition of lipids into energy for embryonic development [61,62]. Thus, these results suggest that methylation modification may influence embryo development by regulating metabolic biosynthesis phases.
Enzyme inactivation via protein phosphorylation regulates embryo development during embryogenesis [44]. Here, six MTGs were identified as being involved in phosphatidylinositol signalling. These genes showed increased methylation and expression level in diapause-terminated eggs when compared with diapause-destined eggs. Ponnuvel et al. reported that enzyme inactivation via protein phosphorylation during early silkworm embryogenesis was followed by dephosphorylation in later stages [44]. On the other hand, INPP genes are important in regulating many cellular activities, including vesicle transport, cytoskeletal dynamics, protein synthesis, cell proliferation and survival and phosphatidylinositol signalling [45]. Here, Inpp5j, Inpp5e and Inpp3a, three of the six MTGs involved in the phosphatidylinositol dephosphorylation pathway, showed increased expression and methylation levels. It was also reported that Inpp5e promotes fundamental metabolism in fat body via phosphorylation [63]. These results indicate that phosphorylation pathways are activated after gene methylation, thereby promoting embryonic signal transduction, material synthesis and other key activities during embryo development.
Our previous data showed that embryonic cells are arrested in G2 at diapause, and cell cycles become slower in proportion to increasing G1 length. Once diapause terminates, the embryos resume development at 25 • C, cells enter M phase, and cell division resumes in the embryos [38]. GO analysis revealed that six MTGs were enriched in the G1/S transition of the mitotic cell cycle. It was interesting that gene 101739208, one of the six MTGs, was a homologue of the G2/M phase-specific E3 ubiquitin-protein ligase (G2E3) and was upregulated and hypermethylated by HCl-treatment. BS-PCR results showed that non-treated and HCl-treated sample methylation levels were 33.3% and 48.1%, respectively for G2E3 (Figure 9), consistent with WGBS data. Previous research has shown that G2E3 is essential for early embryonic development in preventing apoptotic death, and it strengthens the synthesis and transport of genetic material and energy in the M phase [64]. This indicates that hypermethylation is associated with embryonic diapause and may be regulated through the cell cycle and apoptosis inhibition. Pi3k60 homologous gene (GeneID: LOC100158253), which was upregulated and hypermethylated, reportedly promotes the decomposition of lipids into energy for embryonic development [61,62]. Thus, these results suggest that methylation modification may influence embryo development by regulating metabolic biosynthesis phases. Enzyme inactivation via protein phosphorylation regulates embryo development during embryogenesis [44]. Here, six MTGs were identified as being involved in phosphatidylinositol signalling. These genes showed increased methylation and expression level in diapause-terminated eggs when compared with diapause-destined eggs. Ponnuvel et al. reported that enzyme inactivation via protein phosphorylation during early silkworm embryogenesis was followed by dephosphorylation in later stages [44]. On the other hand, INPP genes are important in regulating many cellular activities, including vesicle transport, cytoskeletal dynamics, protein synthesis, cell proliferation and survival and phosphatidylinositol signalling [45]. Here, Inpp5j, Inpp5e and Inpp3a, three of the six MTGs involved in the phosphatidylinositol dephosphorylation pathway, showed increased expression and methylation levels. It was also reported that Inpp5e promotes fundamental metabolism in fat body via phosphorylation [63]. These results indicate that phosphorylation pathways are activated after gene methylation, thereby promoting embryonic signal transduction, material synthesis and other key activities during embryo development.
Our previous data showed that embryonic cells are arrested in G2 at diapause, and cell cycles become slower in proportion to increasing G1 length. Once diapause terminates, the embryos resume development at 25 °C, cells enter M phase, and cell division resumes in the embryos [38]. GO analysis revealed that six MTGs were enriched in the G1/S transition of the mitotic cell cycle. It was interesting that gene 101739208, one of the six MTGs, was a homologue of the G2/M phase-specific E3 ubiquitin-protein ligase (G2E3) and was upregulated and hypermethylated by HCl-treatment. BS-PCR results showed that non-treated and HCl-treated sample methylation levels were 33.3% and 48.1%, respectively for G2E3 (Figure 9), consistent with WGBS data. Previous research has shown that G2E3 is essential for early embryonic development in preventing apoptotic death, and it strengthens the synthesis and transport of genetic material and energy in the M phase [64]. This indicates that hypermethylation is associated with embryonic diapause and may be regulated through the cell cycle and apoptosis inhibition.

GO and KEGG Enrichment Analysis
Gene Ontology (GO) enrichment analyses of DMGs and DEGs were implemented by the GO-seq R package [68], where gene length biases was corrected. GO terms with corrected p-values < 0.05 were considered significantly enriched by DMGs and DEGs. The Kyoto Encyclopaedia of Genes and Genomes (KEGG: http://www.genome.jp/kegg/) [69] is a database that provides high-level function annotation for a diverse range of biological systems. The KOBAS software (version 3.0, Chinese Academy of Sciences, Beijing, China) was used to test for statistical enrichment of DMGs and DEGs in KEGG pathways [70]. KEGG pathway annotation and enrichment analyses were performed similar to GO analysis.

RT-qPCR and McrBC-qPCR Analysis
Total RNA for RNA-seq analyses was also used for RT-qPCR. The primers for RT-qPCR are shown (Supplementary Table S1). RT-qPCR reactions were performed in 25 µL reaction volumes containing SYBR Premix Ex Taq (TaKaRa, Dalian, China) according to manufacturer's instructions. The reaction was performed in a CFX96TM Real-Time System (Bio-Rad, Hercules, CA, USA). Thermal cycling profiles consisted of an initial denaturation at 95 • C for the 30s and 40 cycles at 95 • C for 5 s and 60 • C for 30 s. GO and KEGG pathway analysis identified some important pathways that may be related to diapause; we then selected genes from these pathways for verification from the list of DEGs presented in Supplementary Table S2 . The negative control was performed without GTP for digestion of the standard. qPCR was performed using 20 ng DNA as template. The qPCR process was the same as described earlier, and primers used are shown (Supplementary Table S1). Three independent experiments, with three technical replicates, were performed. All assays were performed in triplicate, and relative expression levels were calculated using the 2 −∆∆Ct method [71]. Statistical analyses were conducted using analysis of variance (ANOVA) and a least significant difference (LSD) posteriori test using SPSS software (version 19.0, IBM, Chicago, USA).

Bisulfite Sequencing Validation of Different Methylation Genes
Genomic DNA was bisulfite-converted using the EZ DNA Methylation-Gold kit (ZYMO, Tustin, CA, USA) according to manufacturer's instructions. Bisulfite-converted DNA was amplified by PCR using Premix Ex Taq™ Hot Start Version (TaKaRa, Dalian, China) according to manufacturer's instructions. Primers for BS-PCR were designed using the online MethPrimer program (http://www. urogene.org/cgibin/methprimer2/Meth Primer.cgi) (Supplementary Table S1). PCR products were purified and cloned into the pMD19-T vector (TaKaRa, Dalian, China), and six clones, selected from each sample, were sent to Sangon Biotech Co. Ltd. (Shanghai, China) for sequencing. The sequencing results were analyzed using the online Quma software (http://quma.cdb.riken.jp/) for individual CG sites and DNA methylation levels.

Conclusions
In conclusion, epigenetic modification is a significant dimension to embryonic development in insects. Our findings revealed that methylation levels in diapause-terminated eggs were higher than in diapause-destined eggs in silkworm and that methylation sites in diapause-terminated eggs tended to be in coding sequences and promoter regions when compared with diapause-destined eggs. Most importantly, our result indicated that methylation was positively associated with embryonic development in silkworm, by regulating cell differentiation, cell cycle, metabolism, apoptosis and phosphorylation. These data will be useful in understanding molecular mechanisms underlying DNA methylation and gene expression changes during embryonic development and insect diapause.

Conflicts of Interest:
The authors declare no conflict of interest.