Intragenic DNA methylation regulates insect gene expression and reproduction through the MBD/Tip60 complex

Summary DNA methylation is an important epigenetic modification. However, the regulations and functions of insect intragenic DNA methylation remain unknown. Here, we demonstrate that a regulatory mechanism involving intragenic DNA methylation controls ovarian and embryonic developmental processes in Bombyx mori. In B. mori, DNA methylation is found near the transcription start site (TSS) of ovarian genes. By promoter activity analysis, we observed that 5′ UTR methylation enhances gene expression. Moreover, methyl-DNA-binding domain protein 2/3 (MBD2/3) binds to the intragenic methyl-CpG fragment and recruits acetyltransferase Tip60 to promote histone H3K27 acetylation and gene expression. Additionally, genome-wide analyses showed that the peak of H3K27 acetylation appears near the TSS of methyl-modified genes, and DNA methylation is enriched in genes involved in protein synthesis in the B. mori ovary, with MBD2/3 knockdown resulting in decreased fecundity. These data uncover a mechanism of gene body methylation for regulating insect gene expression and reproduction.


HIGHLIGHTS
Insect intragenic 5mC enhances gene expression through histone H3K27 acetylation MBD2/3 binds the intragenic 5mC and recruits Tip60 to promote H3K27 acetylation Intragenic 5mCs modify protein synthesis-related genes in insect ovaries The intragenic 5mC plays a role in insect reproduction

INTRODUCTION
Insects are the largest and most diverse group of organisms on Earth; these animals show fecundity, with large numbers of individuals and wide distributions (Mito et al., 2010). Metamorphosis is an insect-specific developmental process in which the organism undergoes changes in internal organs, external characteristics, and even living habitats from larval to adult stages. In addition to development, metamorphosis improves the adaptability of insects to the environment. Insects are classified into ametabolous, hemimetabolous, and holometabolous groups. During metamorphosis, sequential gene expression promotes ovarian development, and dysfunction in transcription causes dramatic defects in ovarian development and reproduction (Bownes, 1989;Swevers and Iatrou, 2003). For example, knockdown of DNA methyltransferase 1 (Dnmt1) in the hemimetabolous species Oncopeltus fasciatus resulted in abnormal ovarian development (Bewick et al., 2019). In addition, RNA interference (RNAi) of Dnmt1 in holometabolous Bombyx mori led to a decrease in the egg hatch rate (Xiang et al., 2013). Therefore, it is important to understand how insect DNA methylation regulates gene expression during ovarian development.
Indeed, DNA methylation is an important epigenetic modification in eukaryotes. Dnmt3 and Dnmt1, which function in de novo methylation and maintenance of existing methylation, respectively, mediate modifications at the fifth carbon of the pyrimidine ring of cytosines in a gene promoter or body, with 5-methylcytosine (5mC) being the most commonly methylated nucleotide in genomic DNA (Goll and Bestor, 2005;Kim et al., 2009). In vertebrates, studies on DNA methylation regulation are often focused on the CpG islands of gene promoter regions (Bird, 2002;Lister et al., 2009), which inhibits gene expression because methylcytosine-binding proteins (MeCP or MBDs) bind to cis-regulatory elements (CREs), blocking the binding of transcription factor(s) (Bird and Wolffe, 1999). MBDs also recruit histone deacetylase (Dobosy and Selker, 2001;Jones et al., 1998;Nan et al., 1998) and histone methyltransferase (Fuks et al., 2003) to silence gene transcription by assembling a repressive nucleosomal array. With regard to gene body methylation in mammals, several reports have demonstrated that gene body 5mC enhances gene expression by promoting In mammals, MBD family proteins consist of MeCP2 and MBD1-6 isoforms. MBDs bind to DNA via a highly conserved MBD domain; MBD3 cannot directly bind to DNA because it lacks the MBD domain. MBD2 and MBD3 repress gene transcription as a part of the nucleosome remodeling and histone deacetylation (NuRD) complex; both proteins directly compete with one another for binding to NuRD components through their intrinsically disordered regions and coiled-coil domains (Baubec et al., 2013;Le Guezennec et al., 2006). Sequence comparison in insects has revealed that the orthologs of mammalian MBD2 and MBD3 are the same protein; this MBD protein with MBD and coiled-coil domains is called MBD2/3 (Hendrich and Tweedie, 2003). Nevertheless, it is unknown whether insect MBD2/3 functions in regulating gene expression in vivo.
Here, we report that insect-specific MBD2/3 binds to 5mC in the gene body of ovarian developmentrelated genes and enhances gene expression to ensure egg development and hatching. We found that in the holometabolous insect B. mori, enrichment of 5mC in ovarian gene bodies promotes transcription via an increase in H3K27ac mediated by the complex of MBD2/3 and acetyltransferase Tip60. We further demonstrate the role of MBD2/3 in the reproduction of holometabolous and hemimetabolous insects. Our study reveals a mechanism for the regulation of insect ovarian development through MBD2/3-mediated intragenic DNA methylation.

RESULTS
5mC is enriched in 5 0 untranslated regions of gene bodies in the B. mori ovary To investigate whether and how intragenic methylation regulates ovarian gene expression, 5mC in the ovarian genome was examined. The presence of 5mC in B. mori ovaries was confirmed by immunofluorescence staining with anti-5mC antibodies and by electrophoretic mobility shift assays (EMSAs) with an intragenic methylated DNA fragment from the B. mori ovary ( Figure 1A). Analyses of spatial expression profiles of both mRNA and protein during pupal stages showed that BmDnmt1, the only DNA methyltransferase in B. mori that introduces 5mC ( Figure S1), was highly expressed in the ovary compared with in the other tissues tested ( Figure 1B). DNA methylation in the B. mori ovary was also analyzed by methylation sequencing. Higher levels of CG DNA methylation, versus at CHG and CHH, were detected. Only 0.52% (303,813 cytosines) of genomic cytosines were methylated, among which 83.61% of mCs were at CG sites and 13.47% and 2.93% were at CHG and CHH sites, respectively ( Figure 1C and Table S1). In particular, more 5mC was detected near the transcription start site (TSS) of gene bodies when compared with other regions. Overall, the sequence with the highest frequency of methylation among methylated genes in the genome was TATCGAT (methyl-preferred fragment) ( Figure 1D). These results suggest that 5mC mainly occurs at CG sites of 5 0 untranslated regions of gene bodies in the B. mori ovary.

Gene body 5mC enhances gene expression
To investigate whether there is a relationship between gene 5mC levels and gene expression, 5mC and gene expression in the ovarian transcriptome were analyzed in the B. mori ovary. The results showed that ll OPEN ACCESS iScience 24, 102040, February 19, 2021 iScience Article most genes with methylation rates higher than 1% also had high expression levels (Figures 2A, left panel  and Tables S2, S3, S4, and S5), similar to the previous results in insects (Suzuki and Bird, 2008;Zemach et al., 2010;Sarda et al., 2012;Glastad et al., 2014). Indeed, intragenic DNA methylation levels correlated positively with gene expression levels, regardless of the DNA methylation rate in the promoter region ( Figure 2A, right panel). To demonstrate whether intragenic DNA methylation influences gene expression, we evaluated mRNA levels in DNA methylation inhibitor-treated ovary-derived Bm12 cells. According to qPCR, the mRNA levels of 14 of 32 genes randomly selected from those with approximately 10% DNA methylation in the ovary (Table S6) were reduced by 5-aza-2 0 -deoxycytidine (5-aza-dC) ( Figure 2B), a DNA methylation inhibitor.
To confirm that upregulation of gene expression is mediated by intragenic DNA methylation, and not through promoter methylation, we determined the transcriptional activity of luciferase reporters in which 0, 2, 4, or 6 repeats of methyl-preferred fragments (TATCGAT) were inserted between the luciferase gene body and actin gene core promoter ( Figure S2 and Table S7). Although luciferase activity increased in the presence of four or more repeated methyl-preferred fragments ( Figure 2C), the luciferase activity of constructs in which these methyl-preferred fragments were inserted in the promoter region or downstream of the luciferase gene body did not increase ( Figure 2D). Additionally, the transcriptional activity of a construct carrying the intragenic methyl-preferred fragment but the mutated methyl-preferred fragment Figure 1. Analyses of DNA methylation in the Bombyx mori ovary (A) 5mC DNA methylation was detected by immunofluorescent staining in the nuclei of B. mori ovaries using an anti-5mC antibody (scale bar, 4 mm) (left) (n = 9). An EMSA shows the specific binding of the anti-5mC antibody to mC probes (right). (B) qRT-PCR (top) and western blot (bottom) analyses of the expression patterns of DNA methyltransferase 1 (BmDnmt1) in different tissues of silkworms during the pupal stage (n = 3). (C) Analyses of whole-genome DNA methylation patterns and distribution in B. mori ovaries of 7-day-old pupae (n = 2). Fraction of mCs identified in each sequence context (CG, CHG, and CHH) based on methylated cytosine percentage in the whole genome (left). Methylation profile in the gene region. The gene region consists of 2 kb upstream and downstream of the TSS and TTS regions and the gene body from the TSS to TTS. TSS: transcription start site; TTS: transcription termination site (middle). Methylation level at different regions (exon and intron) in the gene body (right). (D) Sequence context of CG methylation. CG bases were counted to obtain a base composition of 7 bp near the mCG base in the whole genome, and a Logo map was drawn using seqLogo software. For all measurements, ''n'' represents the number of biological replicates. Data are presented as the mean G SEM. ***p < 0.001 by ANOVA.

OPEN ACCESS
iScience 24, 102040, February 19, 2021 3 iScience Article was enhanced by an increase in DNA methylation due to BmDnmt1 overexpression ( Figure 2E). Conversely, knockdown of BmDnmt1 by RNAi ( Figure 2F) or treatment with a DNA methylation inhibitor (5-aza-dC) blocked the increase in luciferase activity ( Figure 2G). These data suggest that DNA methylation in the gene body enhances gene expression. Figure 2. The relationship between intragenic DNA methylation and gene expression in B. mori ovarian cells (A) Analyses of transcriptome and genomic 5mC levels in B. mori ovaries reveal a positive correlation between DNA methylation and gene expression levels (left), whereby higher gene expression levels are related to intragenic 5mC but not to promoter 5mC (right) (n = 2). R 2 : goodness of fit of the trend lines. *p < 0.05, **p < 0.01, ***p < 0.001 by the F test.
(B) Changes in mRNA levels of 14 genes in Bm12 cells, a cell line of the B. mori ovary, in response to the 5mC inhibitor 5aza-dC. Relative mRNA levels are presented as the gene mRNA level in 5-aza-dC-treated samples divided by the control, which was normalized to the expression level of the housekeeping gene ribosomal protein 49 (RP49) (n = 3).
(C) Effect of the different lengths of methyl-CpG sequences in the 5 0 -terminus of the intragenic region on luciferase activity. The different repeats of the methyl-CpG sequence, ATCGAT, in which methyl-CpG most frequently appears in this sequence of the ovarian genome, were constructed into the luciferase vector between the actin core promoter and the luciferase ORF and cotransfected with the BmDnmt1 overexpression vector for the determination of luciferase activity (n = 3).
(D) Effect of the methyl-CpG sequence at positions of the promoter and 5 0 -or 3 0 -terminus of the gene body on luciferase activity in Bm12 cells (n = 3).
(F) Effect of BmDnmt1 RNAi on intragenic 5mC-regulated luciferase activity in Bm12 cells (n = 3). The BmDnmt1-EGFP vector or dsBmDnmt1 was cotransfected with the vector including six repeats of the methyl-CpG sequence in the gene body for the determination of luciferase activity. In the mutated CpG sequence, ''C'' is changed to ''T.'' (G) Effect of 5-aza-dC (one microgram) on intragenic 5mC-regulated luciferase activity in Bm12 cells transfected with the construct containing six repeats of methyl-CpG sequence in the gene body; ddH 2 O was used as a control (n = 3). Data are presented as the mean G SEM. *p < 0.05, **p < 0.01, ***p < 0.001 by ANOVA. iScience Article Gene body 5mC enhances gene expression through MBD2/3 and Tip60-mediated H3K27 acetylation In mammals, MeCP2 or MBDs link DNA methylation to histone modification by recruiting histone deacetylase (Jones et al., 1998;Nan et al., 1998). Genes encoding two MBDs, BmMBD2/3 and BmMBD4, are found in the B. mori genome. To investigate which MBD recruits histone modification enzymes, binding assays of MBD with DNA were performed. The results showed that some proteins from ovarian nuclear extracts ( Figure 3A, left panel) or purified BmMBD2/3 (GenBank accession no.: XP_004929675.1) ( Figure 3A, middle panel) bound to the methyl-preferred fragment, and this binding was competed by excess unlabeled probe ( Figure 3A). No band was detected in control assays with the unmethylated fragment. Moreover, the complex with the proteins from the ovarian nuclear extract was supershifted by a specific To determine whether BmMBD2/3 and BmTip60 interact with each other to enhance H3K27ac, far-western blot and protein pull-down analyses were performed. The purified recombinant BmTip60-His protein (approximately 68 kDa) was transferred to a nitrocellulose membrane, incubated with the purified BmMBD2/3-GST protein, and then probed with an anti-GST-tagged antibody. BmTip60-His was able to interact with BmMBD2/3-GST, as indicated by the presence of a band recognized by the anti-GST-tagged antibody ( Figure 4A, left panel). Similarly, BmMBD2/3-GST (approximately 54 kDa) was transferred to a membrane and incubated with BmTip60-His, and the anti-His-tag antibody recognized BmTip60-His bound to BmMBD2/3-GST ( Figure 4A, second to left panel), indicating interaction between these two proteins. A similar result was obtained in the GST pull-down experiment: BmTip60-EGFP (approximately 87 kDa) as pulled down by BmMBD2/3-GST and recognized by an anti-EGFP antibody ( Figure 4A, right panel). Furthermore, immunofluorescence showed that the BmMBD2/3 and BmTip60 proteins colocalized in the nucleus of Bm12 cells ( Figure 4B, left panel). Taken together, these data demonstrate that BmMBD2/ 3 and BmTip60 interact and form a complex.
To determine whether interaction between BmMBD2/3 and BmTip60 regulates gene transcription, the luciferase activity of a construct containing a fragment of six repeat methyl-preferred sequences was determined. We found that luciferase activity decreased after BmMBD2/3 or BmTip60 knockdown by RNAi ( Figure 4B, middle panel). Additionally, luciferase activity increased when both BmMBD2/3 and BmTip60 were overexpressed simultaneously in Bm12 cells ( Figure 4B, right panel), suggesting that a complex of the two proteins upregulates gene expression by binding to the methylated gene body region and enhancing H3K27ac.
The aforementioned results demonstrate that B. mori MBD2/3 recruits acetyltransferase instead of deacetylase, as observed in mammals. To confirm the function of BmMBD2/3, the MBD and coiled-coil domains of BmMBD2/3 and Mus musculus MBD3 (MmMBD3) (GenBank accession no.: NP_038623.1) with approximately 50% amino acid sequence similarity ( Figure 4C, top panel) were exchanged to create recombinant Bm:Mm hybrid MBD2/3 proteins ( Figure 4C, left and middle panels), after which luciferase activity under control of the actin promoter and the B. mori methyl-preferred fragment was determined. Luciferase activity was enhanced by overexpressing BmTip60 and MBD2/3 containing the coiled-coil domain from B. mori but not that from MmMBD3 ( Figure 4C, middle panel). In addition, GST pull-down assays confirmed that only MBD2/3 containing the B. mori coiled-coil domain interacted with BmTip60 ( Figure 4C, right panel). These data show that the B. mori MBD2/3 coiled-coil domain is essential for recruiting HAT Tip60. Furthermore, to investigate why B. mori MBD2/3 recruits acetylase, insect coiled-coil domains, which are involved in protein-protein interactions, were analyzed. According to the results, insect coiled-coil domains are not highly conserved with those of mammals, but they are highly similar among insects ( Figure S8), suggesting that the insect-specific coiled-coil domains are responsible for the interaction of MBD and Tip60.
To further demonstrate that the BmMBD/BmTip60 complex binds to mCpG regions in target genes, Chromatin immunopreciptation assays were performed with four DNA-methylated genes, BGIBGMA001106, iScience Article BGIBGMA005486, BGIBGMA009840, and BGIBGMA010975, the mRNA levels of which were reduced by 5aza-dC treatment ( Figure 2B). Bm12 cells were transfected with the recombinant plasmid BmMBD2/3-33FLAG or BmTip60-33FLAG (EGFP vector as a control), and expression of both proteins in the cells was confirmed by western blot ( Figure 5A). In these BmMBD2/3-33FLAG or BmTip60-33FLAG-overexpressing cells, anti-FLAG antibodies, but not control IgG, recognized BmMBD2/3-33FLAG or BmTip60-33FLAG, precipitating and enriching the mCpG regions of four test genes ( Figures 5B and 5C). The enriched mCpG region sequences were amplified by PCR and examined by DNA sequencing ( Figure 5D). These experiments demonstrated that the BmMBD/BmTip60 complex bound to the intragenic mCpG regions of genes in vivo. All these results suggest that BmMBD2/3 binds to intragenic methyl-CpG through its Gene body DNA methylation affects egg production and hatch rates The above-mentioned results show MBD2/3-mediated intragenic DNA methylation in insects. To investigate whether MBD2/3 plays a role in insect reproduction similar to Dnmt1, the expression patterns of MBD2/3 from a set of randomly selected insect species were analyzed. The B. mori MBD2/3 expression pattern ( Figure 6A) was similar to that of BmDnmt1 ( Figure 1B), which was highly expressed in ovaries compared to other tissues. Similar expression patterns of MBD2/3 were also observed in Periplaneta americana, Apis cerana, Spodoptera litura and Apis mellifera ( Figure 6A), suggesting that insect MBD2/3 may iScience Article have a physiological function in egg production and embryonic development. Next, we examined the holometabolous insects B. mori and S. litura and the hemimetabolous insect Locusta migratoria to study the function of DNA methylation. Injection of the 5mC inhibitor 5-aza-dC into B. mori, S. litura, and L. migratoria decreased the number of eggs oviposited ( Figure 6B). Microcomputed tomographic analysis of ovarian morphology in 4-day-old pupae (the middle stage of B. mori pupae) after 5-aza-dC treatment revealed that ovarian development was retarded ( Figure 6C). These results suggest that DNA methylation enhances oocyte development and egg production and subsequent embryonic development. To determine the involvement of Dnmt1 and MBD2/3 in hatching, the dsRNAs targeting Dnmt1 and MBD2/3 were injected into eggs on the first day of the embryonic development. Knockdown of these genes significantly reduced hatch rates in B. mori, S. litura, and L. migratoria compared with the control injected with dsGFP ( Figure 6D). Furthermore, we analyzed the methylated genes in the methylation genome of the B. mori ovary and found that the top 20 enriched methylated pathways include protein synthesis-related RNA transport, ribosome, and protein processing in the endoplasmic reticulum (ER) and spliceosome (Figure S9). These results suggest that intragenic DNA methylation may be involved in the upregulation of genes coding for proteins that play key roles in protein synthesis that function in the production of a large number of eggs during insect reproduction.

DISCUSSION
DNA methylation in the promoter of genes has been extensively investigated for its role in inhibiting gene expression (Jones, 2012). Additionally, 5mC enrichment in the gene body was demonstrated to promote gene expression in cancer cells by regulating alternative splicing via H3K36m3 (Yang et al., 2014), and 5mCs enrichment in gene bodies has been detected in insects (Zemach et al., 2010;Glastad et al., 2014;Hunt et al., 2013;Bonasio et al., 2012). However, the role of intragenic DNA methylation in regulating gene expression and the underlying mechanism remain unknown. In the present study, we discovered the function of gene body DNA methylation in enhancing gene transcription and examined the underlying mechanisms. Intragenic 5mCpG near the TSS enhanced the transcriptional activity of genes ( Figure 2). Treatment with a 5mC inhibitor resulted in a decrease in the mRNA levels of target genes that are intragenically methylated in B. mori ovarian cells ( Figure 2). Furthermore, MBD2/3 bound to methyl-CpG in the intragenic regions of genes via its methyl-CpG-binding domain ( Figure 3A) and then interacted with the acetyltransferase Tip60 through its coiled-coil domain. The MBD2/3-Tip60 complex increased H3K27ac ( Figures  3 and 4). Similar to BmDnmt1 knockdown, RNAi-mediated knockdown of BmMBD2/3 decreased the egg number and hatch rates of B. mori ( Figure 6D). Taken together, these data conclusively show that intragenic DNA methylation upregulates gene expression to ensure a high reproductive rate in insects to enable them to survive challenging environments, raising many intriguing questions for future research.
In insects such as B. mori, analyses of the methylated genome and transcriptome have revealed that intragenic DNA methylation is enriched in highly expressed genes (Xiang et al., 2010). However, experimental evidence for the ability of intragenic 5mC to regulate gene expression is lacking. Our results of promoter activity and mRNA levels in 5mC inhibitor-treated cells showed that 5mC in the 5 0 UTR region of the gene body enhances gene transcription ( Figure 2). In addition, a DNA fragment with more than two CpGs was able to promote gene transcription ( Figure 2). These data suggest that in insects, 5 0 UTR DNA methylation enhances gene transcription and explains why the insect intragenic DNA methylation rate correlates positively with gene expression levels. Nevertheless, not all methylated genes were affected by 5mC inhibitor treatment (Figure 2). Several reasons may explain this result. First, hypermethylation is more likely to occur in highly expressed housekeeping genes, which may be regulated by multiple pathways or mechanisms, including DNA methylation. Therefore, compared with their high background expression levels, it is difficult to detect small changes induced by 5mC inhibitor treatment on the expression levels of these genes. Second, Dnmt1 RNAi or 5mC inhibitor treatment may not affect the intragenic methylation complex (Cook et al., 2019), leading to unaltered gene expression. Hence, it will be interesting to investigate differences in intragenic methylation-regulated gene expression by comparing 5mC inhibitor treatment-effective and -ineffective gene expression in future studies.
We discovered that 5 0 UTR 5mC enhances the promoter activity of genes via H3K27ac. H3K27ac usually functions as an enhancer of gene expression and is enriched around the TSS (Creyghton et al., 2010;Karli c et al., 2010). We found that DNA methylation affected the level of H3K27ac but not that of H3K4me3 and H3K36me3 in B. mori. H3K4me3 is a mark for transcription initiation (Liu et al., 2016), and H3K36me3 is a mark of transcriptional elongation (Yang et al., 2014). We also found that B. mori methyl-CpGs and ll OPEN ACCESS iScience 24, 102040, February 19, 2021 iScience Article H3K27ac are enriched in the region near the TSS (Figure 1), suggesting the relationship between H3K27ac and intragenic DNA methylation. Therefore, this work reveals an intragenic methylation regulation, in contrast to intragenic 5mC, which enhances transcription elongation in the methylated area by regulating alternative splicing (Maunakea et al., 2013;Yang et al., 2014). As 5mC has been found to be enriched in the 5 0 -region of gene bodies in other invertebrates (Zemach et al., 2010;Hunt et al., 2013;Bonasio et al., 2012;Lyko et al., 2010), intragenic DNA methylation might enhance gene transcription in other insects and invertebrates via H3K27ac.
We found that insect intragenic DNA methylation enhances H3K27ac through the methyl-CpG reader MBD2/3 and its partner Tip60. In mammals, MBD2 binds to mCG via the MBD domain and associates with corepressors and other protein complexes through its transcriptional repression domains to inhibit gene expression (Nan et al., 1998;Baubec et al., 2013). MBD3 also recruits deacetylase with its coiledcoil domain to inhibit gene expression (Le Guezennec et al., 2006). These functions of MBD2 and MBD3 appear to be different from those of insect MBD2/3. In fact, MBD2/3 is the only MBD family protein in insects with the same protein structure as mammalian MBD3 (Hendrich and Tweedie, 2003), although the sequence similarity of the coil-coiled domain in insect MBD2/3 and mammalian MBD3 is not highly conserved ( Figure S8). Our results demonstrate that B. mori MBD2/3 forms a complex with ATCGAT and BmTip60 using its MBD and coiled-coil domains (Figures 3 and 4), respectively. However, the coiled-coil domain of mouse MuMBD3 is not able to recruit BmTip60 ( Figure 4). Thus, we identified an insect-specific regulatory mechanism of MBD-mediated DNA methylation. However, several questions about this mechanism need to be addressed in future studies. First, although we demonstrated that the insect MBD coiledcoil domain can recruit HAT Tip60, there is no typical coiled-coil domain in Tip60; thus, further studies are needed to understand how MBD2/3 interacts with Tip60. Second, MBD2/3 is the only MBD protein that binds to 5mC in insects, but it is unknown whether it can recruit deacetylase to inhibit transcription factor binding to the gene promoter in addition to recruiting Tip60 to enhance gene expression. As the NuRD complex is present in D. melanogaster (Zhang et al., 2016), it needs to be determined how MBD2/3 can switch binding between deacetylase and acetyltransferase. Third, it remains to be explored how MBD2/ 3-mediated intragenic DNA methylation regulates housekeeping gene expression to affect insect development. With these questions in mind, the present study, therefore, opens a research perspective for insect epigenetics.
Insects are the largest group of animals on Earth and have high reproduction rates. We found that intragenic DNA methylation regulates insect reproduction. MBD2/3 RNAi affected the reproduction rates of the holometabolous insects B. mori and S. litura as well as the hemimetabolous insect L. migratoria (Figure 6). The effect of MBD2/3 RNAi was similar to that of Dnmt1 RNAi (Figure 6), suggesting that Dnmt1 regulates insect ovarian development, probably through MBD2/3. The effect of Dnmt on egg production and hatch rates has been reported for other insects. For instance, Nlu-Dnmt1 and Nlu-Dnmt3 knockdown led to fewer offspring in N. lugens (Zhang et al., 2015) and O. fasciatus Dnmt1 RNAi affected the number and development of eggs and damaged the nuclear structure of the follicular epithelium (Bewick et al., 2019). In Nasonia, Dnmt1 loss resulted in embryonic lethality during the onset of gastrulation (Zwier et al., 2012). MBD2/3 likely regulates reproduction in other insects in addition to those assessed in this study. Thus Dnmt and MBD2/3 may constitute target genes for pest control.
Why does insect intragenic DNA methylation affect reproduction? We analyzed genes that are methylated in the B. mori ovary and found that methyl-CpGs mostly appeared in the gene bodies of genes encoding proteins involved in protein synthesis, including genes participating in transcription and translation ( Figure S9). Treatment with a 5mC inhibitor at the beginning of ovarian development resulted in retardation of ovary development and a decrease in mature eggs and their hatch rate ( Figure 6). These results suggest that intragenic DNA methylation may enhance the protein synthesis required for egg maturation, including chorion synthesis and metabolism, to support high fecundity and hatch rates. By enhancing gene expression through MBD2/3-mediated H3K27ac, intragenic DNA methylation in the ovary and during early embryo development is necessary for a high reproduction rate in insects. Regardless, how intragenic DNA methylation regulates development-related gene expression to affect ovarian and fertilized egg development needs further investigation. Additionally, whether intragenic DNA methylation affects other tissues in addition to the ovary and fertilized egg needs to be clarified in the future. Although higher 5mC levels appear in the gene bodies of ovaries (this study), silk glands ( iScience Article with 0.52% in ovary, 0.6% in midgut and fat body, and 0.11% in silk gland. Of this 83.61%, 81%, and 99.2% of mCs were at CG sites in the ovary, midgut, and silk gland, respectively, whereas only 13.47% and 2.93% were at CHG and CHH sites in the ovary and 19% mC and 0.8% mCs were at both sites in the midgut and the silk gland, respectively. Unlike in the ovary, 5mCs in CHH sites were higher than at CHG sites in the midgut. These data suggest that there may be different functions and regulatory mechanisms of 5mC in different tissues of B. mori.
In this study, we show a regulatory mechanism of intragenic DNA methylation in insects. MBD2/3 binds to methyl-CpG in intragenic DNA and recruits Tip60, leading to elevated histone H3K27ac levels and enhanced gene expression. This regulatory mechanism of intragenic DNA methylation facilitates a high rate of reproduction.

Limitations of the study
This study reveals a regulatory mechanism that intragenic DNA methylation enhances gene expression and demonstrates the role of intragenic DNA methylation in insect reproduction. However, several questions about this mechanism and function need to be addressed in future studies. First, the complex of methylated DNA and proteins needs to be explored in depth. It needs to understand how MBD2/3 interacts with Tip60 to enhance H3K27 acetylation. Second, it needs to explain why the expression levels of most methylated genes are not changed after demethylation? Whether there are regulatory pathways to compensate for the loss of DNA methylation and to promote these gene expression needs to be clarified. Third, how intragenic DNA methylation regulates development-related gene expression to affect insect ovarian and fertilized egg development needs further investigation.

Resource availability Lead contact
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Sichun Zheng (sczheng@scnu.edu.cn). This study did not generate new unique reagents.

Materials availability
All reagents used in this study are available from the Lead Contact with a completed Materials. All unique reagents generated in this study are available from the Lead contact with a completed Materials Transfer Agreement.

Data and code availability
The raw Illumina sequencing data from the BS-seq and RNA-seq generated during this study are deposited in SRA under accession number SRP217273 at NCBI.

METHODS
All methods can be found in the accompanying Transparent Methods supplemental file.

Library construction and mapping
Genomic DNA was extracted from the ovaries of seven-day-old pupae using the phenol extraction method. After genomic DNA was extracted, DNA concentration and integrity were detected by a NanoDrop spectrophotometer and agarose gel electrophoresis, respectively. DNA libraries for bisulfite sequencing were prepared by Gene Denovo Biotechnology Co (Guangzhou, China). Briefly, genomic DNA was fragmented into 100-300-bp fragments by sonication (Covaris, Woburn, USA) and purified with a MiniElute PCR Purification Kit (QIAGEN, Dusseldorf, Germany).
The fragmented DNAs were end-repaired, and a single "A" nucleotide was added to the 3' end of the blunt fragments. The genomic fragments were ligated to methylated sequencing adapters and fragments with adapters were bisulfite converted using a Methylation-Gold kit (ZYMO, Irvine, USA). Finally, the converted DNA fragments were PCR amplified and sequenced using an Illumina HiSeqTM 2500. Short reads generated by Illumina sequencing were aligned to the B. mori reference genome (Wang et al., 2005;Xia et al., 2004).

Methylation level analysis
To obtain high-quality clean reads, raw reads were filtered by removing reads containing more than 10% unknown nucleotides and low-quality reads containing more than 40% low-quality (Q-value ≤ 20) bases. The obtained clean reads were mapped to the species reference genome using BSMAP software (Xi and Li, 2009) (version: 2.90) by default. Then, a custom Perl script was used to call methylated cytosines, and the methylated cytosines were tested with the correction algorithm described in Lister et al. (Lister et al., 2009). The methylation level was calculated based on the methylated cytosine percentage of each chromosome and in different regions of the genome for each sequence context (CG, CHG, and CHH). To assess different methylation patterns in different genomic regions, the methylation profile at flanking 2-kb regions and gene bodies (or transposable elements) was plotted based on the average methylation levels for each window.

RNA-seq analyses
Total RNA was extracted from the ovaries of 7-day-old pupae using TRIzol reagent (TaKaRa, Dalian, China) according to the manufacturer's instructions.
RNA-seq was performed using the Illumina HiSeq TM 4000 platform by Gene Denovo Biotechnology, and paired reads with an average length of 150 bp were generated.
The clean reads that were filtered from the raw reads were used for mapping to

ChIP-seq analyses
H3K27ac ChIP-seq data were obtained from NCBI with the BioProject accession number PRJNA450142 (Cheng et al., 2018). We used Bowtie2 to map paired-end clean reads to the silkworm reference genome. Highly enriched peaks were obtained by MACS2 using standard settings.

Quantitative real-time PCR (qRT-PCR)
Total RNA was extracted using TRIzol (TaKaRa, Dalian, China), and cDNAs were synthesized using First-Strand cDNA Synthesis Kit (TaKaRa, Dalian, China) following the manufacturer's protocol. qPCR was performed using 2×SYBR Premix EXTaq TM Kit (TaKaRa, Dalian, China). The relative mRNA level of gene expression was calculated by the 2 -△△Ct method (Livak and Schmittgen, 2001). Data were normalized to the housekeeping gene ribosomal protein 49 (Rp49) and analyzed by the 2 -△△Ct method. All data consisted of three biological replicates. The sequences of the qRT-PCR primers used are shown in Table S9.

Construction of the reporter luciferase vector
Genomic  Table S9.

Cell transfection and transcriptional activity determination
The

Expression and purification of recombinant proteins
Open reading frames of select genes were amplified using cDNA or synthesized by TSINGKE (Beijing, China). The cDNAs were subcloned into the pPET-28a, pGEX-6P-1 or EGFP vector with a 6×His, GST tag or GFP tag to generate recombinant expression vectors, respectively. The recombinant proteins were expressed in Escherichia coli (BL21) or Bm12 cells. The specific primers and restriction sites used for each protein expression are shown in Table S9.

Chromatin immunoprecipitation (ChIP) assay
ChIP was performed in Bm12 cells following the instructions of Pierce™  Table S10.

RNA interference (RNAi)
For RNAi in the Bm12 cell line, a 400~600-bp unique fragment of the ORF of target genes was chosen as the template for synthesizing gene-specific dsRNA using  Table S9.

Western blot and far-western blot analyses
Tissues or Bm12 cells were homogenized in Cell Lysis Buffer (Beyotime, Shanghai, China). For western blot, proteins (30~100 μg) were mixed with 5× loading buffer and then separated by 12% SDS-PAGE, followed by transfer to a nitrocellulose blot membrane. The membrane was blocked with Tris-buffered saline (pH 9.0) with 0.5% Tween-20 (TBST) containing 3% (w/v) BSA, followed by hybridization overnight in TBST with 1% BSA and the primary antibody. After the membranes were washed in TBST three times, a horseradish peroxidase (HRP)-conjugated antibody (Dingguo Biotechnology, Beijing, China) was applied. An anti-tubulin antibody (Dingguo Biotechnology, Beijing, China) was used to verify equal loading of the protein samples on the gel. The antibodies used in this study are listed in Table   S10.
For far-western blot, the purified recombinant BmTip60-His or BmMBD2/3-GST protein (3 mg) was separated by 12% SDS-PAGE and transferred to nitrocellulose membranes. The membranes were washed in TBST to remove SDS, blocked with 3% (wt/vol) BSA in TBST, and then probed with another purified protein (2 mg/mL), BmMBD2/3-GST or BmTip60-His, in TBST containing 0.1% BSA. The blots were incubated overnight; after washing, the membranes were probed with antibodies against the GST or His tag, followed by a horseradish peroxidase (HRP)-conjugated antibody.

Immunofluorescence
Bm12 cells at the logarithmic growth phase were subcultured on glass coverslips The nuclei were stained with DAPI (Beyotime, Shanghai, China) for 10 min. Wing tissues stained with anti-5mC and DAPI were observed and imaged using an FV3000 confocal microscope (Olympus). The antibodies used in this study are listed in Table   S10.

Bioinformatics analysis of protein sequences
Phylogenetic tree analysis was conducted using MEGA7.

Egg number and hatch rate statistics
To record egg number, B. mori and S. litura were injected with 5-aza-dC on the first day of pupation, and L. migratoria was injected with 5-aza-dC on the first day of adulthood. Ovaries of B. mori, S. litura or L. migratoria were dissected from adult females on the first, third or tenth day. Eggs were collected from females treated with 5-aza-dC or ddH2O (control) and counted.
To calculate the hatch rate, dsRNA (20-50 ng) was injected into B. mori, S. litura and L. migratoria embryos within 2 h after fertilization. The larvae hatched from the resulting eggs of B. mori and L. migratoria were recorded at 9-13 days after dsRNA injection; the hatched larvae of S. litura were recorded at 3-5 days after dsRNA injection. The hatch rate was calculated as follows: hatch rate=larval number/egg number.

Statistical analysis
Data are presented as the mean ± SEM. "n" represents the number of biological replicates. P values for the purpose of group comparisons were calculated using ANOVA (*P < 0.05, **P < 0.01, ***P < 0.001). Figure Table S2.  Table S9.   "n" represents the number of biological replicates. Data are presented as mean ± SEM.