Complete Metamorphosis in Manduca sexta Involves Specific Changes in DNA Methylation Patterns

The transition between morphologically distinct phenotypes during complete metamorphosis in holometabolous insects is accompanied by fundamental transcriptional reprogramming. Using the tobacco hornworm (Manduca sexta), a powerful model for the analysis of insect evolution and development, we conducted a genome-wide comparative analysis of gene expression and DNA methylation in caterpillars and adults to determine whether complete metamorphosis has an epigenetic basis in this species. Bisulfite sequencing indicated a generally low level of DNA methylation with a unimodal CpGO/E distribution. Expression analysis revealed that 24 % of all known M. sexta genes (3.729) were upregulated in last-instar larvae relative to the adult moth, whereas 26 % (4.077) were downregulated. We also identified 4.946 loci and 4.960 regions showing stage-specific differential methylation. Interestingly, genes encoding histone acetyltransferases and histone deacetylases were differentially methylated in the larvae and adults, indicating there is crosstalk between different epigenetic mechanisms. The distinct sets of methylated genes in M. sexta larvae and adults suggest that complete metamorphosis involves epigenetic modifications associated with profound transcriptional reprogramming, involving approximately half of all the genes in this species.


INTRODUCTION
Transformation from one morphologically distinct phenotype to another during development is known as metamorphosis. The complete metamorphosis of holometabolous insects involves a pupal stage, during which larval tissues and organs are remodeled and/or replaced by the structures required in the imagoes. A fundamental question in developmental biology is how the same genome can generate morphologically and ecologically distinct phenotypes, such as the caterpillars, pupae, and imagoes in the order Lepidoptera. Many if not all developmental processes involve the orchestrated large-scale transcriptional reprogramming of genes. We therefore postulate that the metamorphosis of holometabolous insect larvae into adults involves fundamental changes in gene expression, and that epigenetic mechanisms may be involved in this process. Epigenetics can explain the phenotypic diversity of insect larvae, pupae, and adults sharing a common genome because such mechanisms globally modulate gene expression rather than altering the DNA sequence (Belles, 2017;Glastad et al., 2019).
There are three major epigenetic mechanisms, one of which is the expression of microRNAs. These short, non-coding RNAs (18-24 bp) operate at the post-transcriptional level to inhibit the translation of specific mRNAs by base-pairing with the untranslated regions or occasionally the coding region (Asgari, 2013;Hussain and Asgari, 2014). Individual microRNAs can regulate the expression of single genes or hundreds of genes. The role of microRNAs in the regulation of complete metamorphosis in holometabolous insects is well-established (Belles, 2017;Ylla et al., 2017). The first microRNAs that are differentially expressed during lepidopteran metamorphosis have been identified in the greater wax moth Galleria mellonella (Mukherjee and Vilcinskas, 2014).
The two other principal epigenetic mechanisms regulate transcriptional initiation. The first is the acetylation and deacetylation of histones by enzymes with opposing activities, thus controlling the ability of transcription factors to access chromatin and initiate gene expression. The addition of acetyl groups to histones is mediated by histone acetyltransferases (HATs), which enhance access to DNA by loosening the chromatin, whereas the removal of acetyl groups by histone deacetylases (HDACs) has the opposite effect and therefore causes gene silencing (Marks et al., 2003). We have previously shown that histone acetylation/deacetylation modulates gene expression during the complete metamorphosis of G. mellonella (Mukherjee et al., 2012). The involvement of both histone modification and microRNAs in this process provides evidence for cross-talk between different epigenetic mechanisms (Mukherjee and Vilcinskas, 2014).
The final major epigenetic mechanism is DNA methylation (Glastad et al., 2019). The addition of a methyl group to a cytosine residue in the dinucleotide sequence CpG results in the formation of 5-methylcytosine, which retains the base-pairing specificity of the unmodified nucleoside but influences its interactions with regulatory proteins. The transfer of methyl groups to DNA is mediated by evolutionarily-conserved enzymes collectively known as DNA methyltransferases (DNMTs). These can be divided further into maintenance methyltransferases (DNMT1), which complete the symmetrical methylation marks on newly-replicated DNA by recognizing the hemimethylated sequences inherited from each parent, and de novo methyltransferases (DNMT3), which establish new methylation marks on unmethylated DNA (Bestor, 2000;Klose and Bird, 2006). In insects, DNMT1 is evolutionarily conserved and widely distributed whereas DNMT3 is restricted to some species (Bewick et al., 2017). Despite the apparent loss of DNMT3 in the order Lepidoptera, these insects have functional methylation systems (Provataris et al., 2018). Indeed, the DNA of holometabolous insects is methylated, although typically only in exon sequences (Provataris et al., 2018).
Here, we tested the hypothesis that metamorphosis-related transcriptional reprogramming in the lepidopteran model species M. sexta involves DNA methylation. We carried out a genome-wide comparative analysis of gene expression and DNA methylation in caterpillars and adults, and identified correlations between stage-specific genes and methylation marks, including marks affecting genes involved in other epigenetic processes.
Our data provide insight into the potential role of epigenetic modifications of complete metamorphosis of holometabolous insects.

Methylation in Manduca sexta
Comparing DNA methylation in different developmental stages in M. sexta, the mean mapping efficiency for Bismark mapping for paired end mapping was 45.6 % and could be increased to 57.2 % due to the following single end mapping (Supplementary Table 1). The HISAT2 mean mapping for transcriptome sequencing results was 67.6 % (Supplementary Table 1).
The overall level of DNA methylation in each sample (L1-3, A1-3) was 0.18-0.48 %, but when focusing on specific positions (found in two out of three larvae or adults) the level of DNA methylation was much lower (0.05 %), as shown in Supplementary Material (Supplementary Table 2). To gain more insight into the methylated sites in the M. sexta genome, we compared the positions based on the sequence context (CpG, CHG, or CHH) and the sample type (larvae, adults, or both). This showed that the CpG context was the most prominent (62 % of all methylated sites) compared to CHG at 6 % and CHH at 32 % ( Figure 1A and Supplementary Table 3). The positions were almost equally distributed in all three samples: 28 % in larvae, 37 % in adults, and 35 % in both ( Figure 1B and Supplementary Table 3). The normalized CpG content (CpG O/E ) revealed a unimodal distribution (Figure 2). Given the relatively large number of methylated positions in the sequence contexts CHG and CHH ( Figure 1A), we carried out a global analysis of differential methylation in M. sexta (CX) rather than focusing on CpG, CHG or CHH methylation alone. We identified 4946 differentially methylated loci (DMLs) and 4960 differentially methylated regions (DMRs) in total. These were classified according to their genomic location as belonging to the gene body (gb), regulatory region (rr), or intergenic region (ir). Interestingly, metamorphosis was accompanied by a significant shift in the location of DMLs and DMRs, resulting in an increase in the number of methylated sites in gene bodies at the expense of sites especially in the intergenic regions, whereas the number of methylated sites in regulatory regions was almost the same in both developmental stages (Figure 3 and Supplementary Table 4).

Differential Gene Expression
The annotated M. sexta genome sequence contains 15.542 genes (Kanost et al., 2016a,b,c). We compared the expression profiles of all these genes in larvae and adults, and found that about half of them were either not expressed at detectable levels during either stage (n.d., 35 %) or showed no evidence of significant differential expression (n.s., 15 %). The remaining genes were either   Table 5). Gene ontology enrichment analysis revealed eight functional categories with stage-specific over-representation, seven of which were related to reproduction and one involved in sensory perception (Supplementary Table 6).

Correlation Between Differential Methylation and Gene Expression
To investigate the correlation between differential methylation and gene expression, we selected all genes expressed at detectable levels at one or both developmental stages that also contained DMLs or DMRs (without distinguishing between the type of methylated site). Given the possibility that genes might overlap with multiple DMLs or DMRs, we also distinguished between genes that were exclusively methylated in the larvae or adults (L_5mC or A_5mC), methylated more extensively but not exclusively in the larvae or imagoes (pL_5mC or pA_5mC), and those methylated to the same extent in larvae and adults (mixed). As shown in Table 1, most of the differentially expressed genes showed no evidence of methylation (74 %). After this, the most abundant categories were the exclusively differentially methylated genes (L_5mC 14 % and A_5mC 10 %). Only a small proportion of the genes were assigned to the mixed, pL_5mC, and pA_5mC categories.
For the investigation of a potential correlation between methylation pattern and expression status, we analyzed the log fold change for all differential expressed genes in dependency of their methylation status. With an exception of pL_5mC all medians for log fold changes are negative and indicate a down regulation of the gene expression. For pL_5mC this trend is reversed and the median of the log fold changes is positive and represents an up regulation (Figure 4).

Genes of Interest
To understand the functional implications of the correlation between gene expression and DNA methylation during metamorphosis, we focused on 850 M. sexta genes representing eight different functional categories: chitin metabolism, detoxification, immunity, epigenetic modification, juvenile hormone/ecdysone signaling, chorion proteins, and vitellogenin (Supplementary Table 7). The genes coding for chorion proteins and vitellogenin showed no evidence of differential methylation. Where differential methylation was detected, in general there was no correlation between gene expression levels and the degree of DNA methylation at each developmental stage ( Table 2 and  Supplementary Table 8). For example, seven genes encoding chitin-binding proteins were higher expressed in larvae, three of which were hypermethylated in larvae and four hypermethylated in imagoes. However, one exceptional case was the immunityrelated gene Atg8 (Msex2.12227), which produces two separate transcripts encoding protein isoforms involved in autophagy. Transcript Msex2.12227-RA was upregulated in larvae relative to adults, whereas transcript Msex2.12227-RB was downregulated. The genomic region for that gene contain three differentially methylated regions: two hypermethylated in adults covering both transcripts and one hypermethylated in larvae exclusive covering the transcript Msex2.12227-RA. We also identified genes with roles in histone acetylation/deacetylation that were upregulated in imagoes and generally hypermethylated in the larvae, with the exception of one HDAC gene that hypermethylated in adults. 1 | Classification of gene expression based on differential methylation analysis.

Up
Down Not significant All Genes are either upregulated (higher expression in larvae than adults), downregulated (higher expression in adults than larvae) or not significant (expressed at similar levels at both stages). For each gene, the methylation status was assigned as (probably) higher methylation in adults (p)A_5mC, (probably) higher methylation in larvae (p)L_5mC, mixed methylation signature with equal methylation counts higher in adults and larvae, or not differentially methylated.
FIGURE 4 | Correlation between gene expression (log2 fold change) and the different possible groups based on methylation status per gene. The methylation status of each gene was defined as mixed, for equal numbers of methylated loci and regions in larvae and adults, as higher methylation in larvae (L_5mC) or higher methylation in adults (A_5mC) where methylation was exclusively state-specific, or as probably higher methylation in larvae (pL_5mC) or probably higher methylation in adults (pA_5mC) for loci/regions methylated in both samples but with a preference for one stage over the other.

DISCUSSION
DNA methylation is involved in the regulation of many physiological and developmental processes in insects, but its role in metamorphosis has not been investigated in detail. Previous studies have shown that the transcriptional reprogramming that accompanies metamorphosis is associated with other forms of epigenetic regulation, specifically microRNA expression and the acetylation/deacetylation of histones (Mukherjee et al., 2012;Mukherjee and Vilcinskas, 2014).  , and not significantly differentially expressed (n.s.) genes are listed according to the type of differential methylation: higher methylation in larvae (L_5mC), probably higher methylation in larvae (pL_5mC), higher methylation in adults (A_5mC), probably higher methylation in adults (pA_5mC), or equal amount of differential methylated regions (mixed).
Therefore, we investigated the relationship between genomic DNA methylation and developmentally-regulated gene expression during metamorphosis in the model lepidopteran species M. sexta. Initially we measured the general prevalence of DNA methylation in the M. sexta genome and found that the overall level of methylation was low. This is in line with recent findings in other holometabolous insects, which tend to have lower levels of DNA methylation in exons in comparison to hemimetabolous insects (Provataris et al., 2018). Accordingly, CpG O/E analysis revealed a unimodal distribution which indicates a limited degree of C to T transitions over evolutionary timescales, consistent with a low level of CpG methylation (Provataris et al., 2018). Similar unimodal distributions have been reported in the silkworm Bombyx mori, the red flour beetle Tribolium castaneum, and the large milkweed bug Oncopeltus fasciatus (Bewick et al., 2017) as well as the mosquito Anopheles gambiae (Bewick et al., 2017). Two different distributions in one specimen were found in the ant species Camponotus floridanus and Harpegnathos saltator, which also show low DNA methylation levels. While CpG dinucleotides have a bimodal distribution, the CHH and CHG context showed a unimodal one (Bonasio et al., 2012). In contrast, the honeybee Apis mellifera features a bimodal CpG O/E distribution reflecting DNA methylation associated with differential gene expression in different bee castes (Elango et al., 2009). However, positive correlation between gene expression and cytosines located in CpG context in species with a unimodal CpG O/E and low methylation levels, like it was observed for T. castaneum (Song et al., 2017). Consequently, the sparse DNA methylation in M. sexta does not exclude the possibility of a role during complete metamorphosis. We therefore expanded our comparative analysis of DNA methylation in larvae and adults beyond the CpG O/E distribution to include the separate analysis of sites located in gene bodies, regulatory regions, and intergenic regions. This revealed that both stages showed preponderance for DMLs and DMRs in the gene body. Nevertheless, a significant shift during metamorphosis was detected with a greater proportion of gene body sites in the adult moths primarily at the expense of intergenic sites, which were more prevalent in the larvae. The number of sites in gene regulatory regions remained almost steady throughout development. These observations are in contrast with the report of Song et al. (2017) showing that more than 80 % of the methylated cytosine residues in T. castaneum are almost entirely restricted to intergenic regions (Song et al., 2017).
We were particularly interested to probe the relationship between differential methylation and differential gene expression during metamorphosis in M. sexta. We found that almost half of all genes were differentially expressed, the remainder either are expressed at similar levels in both larvae and adults or not expressed at detectable levels in either stage. There was little overall correlation between the methylation patterns and gene expression profiles, but there were a small number of striking correlations for highly relevant genes. One of the differentially methylated genes in M. sexta encoded chitin synthase 2. In B. mori a gene from the same class is also differentially methylated during development, with the methylation mark being erased from the promoter in the adult allowing its expression in the wings (Xu et al., 2018). We also observed a correlation between expression and methylation for two different transcripts of the atg8, which is involved in autophagy. Isoform RA was upregulated in larvae and the corresponding genomic region was potentially hypermethylated in imagoes, whereas isoform RB was downregulated in larvae and the corresponding genomic region was exclusively hypermethylated in the adults. Similarly, G. mellonella Atg8 is induced during metamorphosis, with the increasing abundance of the transcript and protein suggesting a role in midgut transformation (Khoa and Takeda, 2012). Combined with our findings, the data suggest that DNA methylation contributes to the regulation of Atg8 during metamorphosis in the Lepidoptera. Juvenile hormone and ecdysone are particularly important regulators in insect development, and the corresponding genes were recently found to be regulated by histone acetylation and deacetylation (Roy and Palli, 2018). We therefore focused on the corresponding M. sexta genes, but only two of these were differentially methylated genes and there was no relationship with the observed expression profiles. The absence of a link between DNA methylation and juvenile hormone has also been demonstrated in honeybees (Cardoso-Junior et al., 2018). The latter study also showed that DNA hypomethylation led to an increase in vitellogenin expression, which is important for reproduction. In contrast, DNA hypomethylation in the brown planthopper Nilaparvata lugens caused a loss of fecundity (Zhang et al., 2015). In M. sexta, we found no evidence that genes encoding chorion proteins or vitellogenin were differentially methylated, suggesting they are not regulated directly by this epigenetic mechanism.
The expression of differentially methylated genes related to the RNA interference pathway declined during the transition from larvae to imagoes, suggesting that RNA-dependent regulation may play a more prominent role in M. sexta larvae than adults. We identified few differentially methylated genes involved in detoxification and immunity-related functions, but interestingly we found that some genes with roles in histone acetylation/deacetylation were upregulated in adults and generally hypermethylated in the larvae, with the exception of one HDAC gene that hypermethylated in imagoes. This indicates there is crosstalk between different epigenetic mechanisms and that histone modification may work in concert with DNA methylation during metamorphosis in M. sexta. Similarly, we found that the specific inhibition of HATs in G. mellonella delays pupation whereas the inhibition of HDACs leads to precocious pupation and accelerated development (Mukherjee et al., 2012). Given that the opposing activities of HATs and HDACs are tightly regulated, we conclude that differential methylation of the corresponding genes has the potential to contribute to the epigenetic regulation of transcriptional reprogramming during metamorphosis in M. sexta and other lepidopterans. Furthermore, histone acetylation and deacetylation in T. castaneum influence the production of juvenile hormone (Roy and Palli, 2018). Thus, our findings suggest that DNA methylation could indirectly regulate the juvenile hormone system. Metamorphosis in the oyster Crassostrea gigas is accompanied by changes in the methylation status of exons, introns, promoters, repeats and transposons, including the increased methylation of exons at the expense of introns, which contrasts with our findings in M. sexta (Riviere et al., 2017). The same study revealed that, regardless of the level of methylation, the pattern of methylation within genes is associated with transcriptional regulation. Another recent study revealed that the model cephalochordate Amphioxus lanceolatum has also low levels of CpG methylation, but this nevertheless is involved in the regulation of metamorphosis (Marletaz et al., 2018). Changes in DNA methylation patterns therefore have the capability to accompany metamorphosis in all organisms that contain methylated cytosines, but there is immense diversity in terms of the methylation targets: gene bodies, regulatory regions, and intergenic regions. Even low levels of methylation can help to regulate the transcriptional reprograming that occurs during metamorphosis in holometabolous animals, potentially in an indirect manner involving other epigenetic mechanisms. We have provided the first evidence that DNA methylation is potentially involved in the complete metamorphosis of M. sexta, a model lepidopteran.

Sample Preparation
Samples were taken from our in-house M. sexta stock, which was reared as previously described (Gegner et al., 2019). We selected three female L5 during their feeding stage (2-3 days after molting) and three female adults 2-3 days after their eclosion. Individual insects were frozen and ground in liquid nitrogen. RNA was isolated from 50 mg of ground tissue of one animal using the Direct-zol RNA MiniPrep Kit (ZymoResearch) including DNA digestion according to the manufacturer's recommendations. Sample quantity and purity were assessed using a NanoDrop ND-1000 UV/Vis spectrophotometer (Thermo Fisher Scientific) and the integrity was verified using an Agilent 2100 Bioanalyzer and a RNA 6000 Nano Kit (Agilent Technologies). DNA was extracted from equal amounts of tissue from the same individual using the Genomictip 20/G and Genomic DNA Buffer Set (Qiagen) according to the manufacturer's protocol. RNA-Seq using poly(A)+ enriched RNA fragmented to an average of 250 bp and bisulfite sequencing were carried out using the llumina 2000 HiSeq2500 platform at GATC Biotech (Eurofins Genomics) based on single samples of 2 µg RNA or 1 µg DNA, respectively. Sample preparation for Bisulfite and RNA sequencing was performed by GATC Biotech using their in house protocols.
The M. sexta draft genome sequence (Kanost et al., 2016b) was used as a reference (Kanost et al., 2016a,b,c). Reads were mapped by paired-end mapping in Bismark v0.15.0 (Krueger and Andrews, 2011) using Bowtie2 v2.3.4.1 with default settings (Langmead and Salzberg, 2012). The software was combined into a Bismark Docker image (Förster, 2018a) tagged v0.15.0. Unmapped reads were remapped separately by single-end mapping, and all mapped reads were deduplicated using the Perl script deduplicate_bismark_alignment_output.pl from Bismark software package. Positions of one biological sample with coverage > 10 and ≥ 10% methylated cytosine residues were defined as methylated . If positions were called methylated in two out of three biological samples per developmental stage, they were defined as methylated within that stage. Differential methylation analysis was processed using DSS v2.26.0 (Wu et al., 2013(Wu et al., , 2015Feng et al., 2014;Park andWu, 2016) in R v3.4.4 (R Development Core Team, 2020) with default settings and an unspecified cytosine context (CX). The p-value threshold was set to 0.001 for differentially methylated loci and to 0.01 for differentially methylated regions. The DSS installation was combined into a DSS Docker image (Förster, 2018c) and tagged v2.28.0.

Characterization of Methylated Regions
The original genome annotation file contained only some annotations for 5′ and 3′ untranslated regions (UTRs). Therefore, the complete file was re-analyzed and annotations for 5′ and 3′ UTRs were added to the transcript annotation, if the coding region was shorter than the complete transcript. The genome was sectioned into three sequence categories: gene bodies (gb: exons, introns, 5′ and 3′ UTRs), regulatory regions (rr: sequences 2 kbp beyond the UTRs), and intergenic regions (ir). DMLs and DMRs were characterized by their location in the genome and the expression of any genes. These results were summarized per gene and the methylation status of each gene was defined as mixed, for equal numbers of methylated loci and regions in larvae and imagoes, as higher methylation in larvae (L_5mC) or higher methylation in adults (A_5mC) where methylation was exclusively state-specific, or as probably higher methylation in larvae (pL_5mC) or probably higher methylation in adults (pA_5mC) for loci/regions methylated in both samples but with a preference for one stage over the other.

RNA Raw Read Processing and Differential Expression Analysis
Quality control measures, including the filtering of high-quality reads based on fastq file scores, the removal of reads containing primer/adapter sequences, and trimming of the read length, were carried out using CLC Genomics Workbench v8.1 (Qiagen Bioinformatics, 2018). Predicted transcript sequences in the M. sexta draft genome (Kanost et al., 2016a) sequence were downloaded from (Kanost et al., 2016c) and matched to our RNA-Seq data using HISAT2 v2.1.0 (Kim et al., 2015). The software was combined into a HISAT2 Docker image (Förster, 2018d) and tagged v2.1.0. The resulting SAM files were converted into BAM files using samtools v1.9 (Li et al., 2009). The BAM files were processed by stringtie v1.3.4d (Pertea et al., 2015(Pertea et al., , 2016 to generate count information for each transcript and sample. Stringtie includes the script prepDE.py which was used to prepare the mapping counts for the differential expression analysis. Docker images are available containing samtools (Förster, 2018g) and stringtie (Förster, 2018h), tagged v1.9 and v1.3.4d, respectively. Differential expression analysis was performed using DESeq2 v1.14.1 (Love et al., 2014) for GNU R v3.4.4, which is available as a Docker image (Förster, 2018b) tagged v1.14.1. All transcripts with normalized count sums of < 1.000 across all samples were defined as not detected (n.d.). Transcripts with adjusted p ≥ 0.05 showed no significant differential expression profiles (n.s.). Where the adjusted p-values were < 0.05, a log2 fold change < 0 reflected the upregulation of the corresponding gene in larvae compared to imagoes whereas a log2 fold change > 0 reflected the downregulation of that gene in larvae compared to adults. Genes were classified as modulated when the expression of at least one transcript was upregulated or downregulated. If the same gene was simultaneously represented by upregulated and downregulated transcripts, the expression class was defined by the mRNA with the highest raw mean value.

CpG O/E Analysis
The normalized CpG dinucleotide content was used as a proxy for the presence of DNA methylation because methylated cytosine residues are prone to spontaneous deamination into thymine, leading to a gradual decline in the prevalence of CpG dinucleotides, a phenomenon known as CpG depletion (Provataris et al., 2018). Based on the published gene set, CpG O/E analysis was performed according to Bird (1980). We used the script cpgoe.pl with the gene annotation file and the draft genome scaffolds as inputs.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA517460.

AUTHOR CONTRIBUTIONS
JG and FF contributed to analysis and interpretation of data. HV, AB, and AV designed the study. AV inspired the study and acquired funding. All authors contributed to the drafting of the manuscript and gave their final approval.