Environmental stress during larval development induces DNA methylation shifts in the migratory painted lady butterfly (Vanessa cardui)

Seasonal environmental fluctuations provide formidable challenges for living organisms, especially small ectotherms such as butterflies. A common strategy to cope with harsh environments is to enter diapause, but some species avoid unsuitable conditions by migrating. Despite a growing understanding of migration in the life cycles of some butterfly species, it remains unknown how individuals register and store environmental cues to determine whether and where to migrate. Here, we explored how competition and host plant availability during larval development affect patterns of DNA methylation in the migratory painted lady (Vanessa cardui) butterfly. We identify a set of potentially functional methylome shifts associated with differences in the environment, indicating that DNA methylation is involved in the response to different conditions during larval development. By analysing the transcriptome for the same samples used for methylation profiling, we also uncovered a non‐monotonic relationship between gene body methylation and gene expression. Our results provide a starting point for understanding the interplay between DNA methylation and gene expression in butterflies in general and how differences in environmental conditions during development can trigger unique epigenetic marks that might be important for behavioural decisions in the adult stage.

tions by migrating. Despite a growing understanding of migration in the life cycles of some butterfly species, it remains unknown how individuals register and store environmental cues to determine whether and where to migrate. Here, we explored how competition and host plant availability during larval development affect patterns of DNA methylation in the migratory painted lady (Vanessa cardui) butterfly. We identify a set of potentially functional methylome shifts associated with differences in the environment, indicating that DNA methylation is involved in the response to different conditions during larval development. By analysing the transcriptome for the same samples used for methylation profiling, we also uncovered a non-monotonic relationship between gene body methylation and gene expression. Our results provide a starting point for understanding the interplay between DNA methylation and gene expression in butterflies in general and how differences in environmental conditions during development can trigger unique epigenetic marks that might be important for behavioural decisions in the adult stage.

K E Y W O R D S
DNA methylation, environmental stress, migration, painted lady butterfly, Vanessa cardui machinery underlying migratory behaviour is somewhat less known, with research focused on a few species of moths and the monarch (Danaus plexippus) butterfly (Chowdhury et al., 2022;Jiang et al., 2011;Reppert & de Roode, 2018).
Painted ladies undertake an annual migratory cycle that can span ~6-10 generations. The most well-studied circuit runs from tropical Africa to northern Europe and back (Menchetti et al., 2019;Talavera et al., 2018;Talavera & Vila, 2017). Since the migratory movement of the painted lady is multi-generational, the direction and duration of migratory flight cannot be 'learned' by conspecifics. Instead, individuals must use environmental cues to determine whether and where to migrate. It is not known at what life stage the decisions are made.
An important environmental cue is likely host plant abundance (Hu et al., 2021;López-Mañas et al., 2022;Stefanescu et al., 2021), which hypothetically could be sensed already in the early larval stages.
In the armyworm (Mythimna separata) moth, for example, periodic starvation during development increases migration tendencies in the adults (Cao et al., 1996). Both the abundance of food and the density of competing conspecifics could inform the developing larvae about the expected future availability of host plants. Higher larval rearing densities seem to enhance flight phenotypes in most, but not all, Lepidoptera and other insects showing migration tendencies (Du et al., 2022;Jiang et al., 2011;Lü et al., 2022;Parker & Gatehouse, 1985).
Information about environmental signals such as density of conspecifics or food scarcity needs to be stored, processed and transmitted during metamorphosis to influence migratory behaviour of the adult butterfly. One category of molecular mechanisms that fulfils the requirements of information storage and transmission is epigenetic modifications, especially 5′-CpG-3′ DNA methylation (hereafter 'methylation'), which can be transmitted to daughter strands during DNA replication with high fidelity (Catania et al., 2020). In plants and vertebrates, methylation of promoter sequences has been linked to transcriptional repression (Bird, 2002;Sarkies, 2022;Suzuki & Bird, 2008;Yi, 2017), while the functional role of gene body methylation is still under debate (Jjingo et al., 2012;Neri et al., 2017;Zilberman, 2017). In arthropods, patterns of methylation are less characterized and potential functional relationships with transcription have been investigated to some extent but results remain inconclusive (Bonasio et al., 2012;Hunt et al., 2013;Jones et al., 2018;Lewis et al., 2020;Liu et al., 2019;Oldroyd & Yagound, 2021;Xiang et al., 2010Xiang et al., , 2013Xu et al., 2018Xu et al., , 2021. A study in the crustacean Daphnia pulex revealed that DNA methylation patterns induced by a stressful environment can be transgenerationally inherited (Harney et al., 2022). It remains to be shown if this is a general pattern in arthropods. Recent initiatives suggest that the genome-wide levels of methylation can vary from moderate (~ 20%) in some arthropod species to almost undetectable (~ 0%) in others (Bewick et al., 2017;Lewis et al., 2020;Oldroyd & Yagound, 2021). An analysis of silk gland tissue of silk moth (Bombyx mori) larvae revealed a low genome-wide methylation level (0.11%) and a positive relationship between gene expression and gene body methylation (Xiang et al., 2010). An investigation of potential relationships between methylation and expression in adult cotton bollworm moths (Helicoverpa armigera) with different flight performance phenotypes pinpointed some differentially methylated CpG sites (Jones et al., 2018). However, the experiment covered only 16 genes and the data were not sufficient to establish if differentially expressed genes were enriched with differentially methylated CpGs. A general aim of our study was therefore to cover this gap by studying genome-wide associations between gene expression and DNA methylation in general. We combine that with specific analyses of methylation and expression differences between individuals reared under different environmental conditions.
Here we investigate the DNA methylome of the painted lady using whole-genome bisulphite sequencing of experimental cohorts exposed to different environmental cues that potentially can affect the propensity of individuals to migrate. Our study consists of three components. First, we map general patterns of methylation among genomic features. Second, we combine the methylation data with gene expression profiling of the same samples to uncover the relationship between methylation and gene expression at an unprecedented scale in butterflies. Third, we test if DNA methylation can be a mechanism to transfer environmental signals experienced during larval development to the imaginal stage in female butterflies by splitting larvae into three treatment groups with different larval densities and food availabilities. Since we hypothesise that neurogenetics has a major role in 'decision-making', we analyse heads (including antennae), where the nervous system of butterflies is con- 2 | MATERIAL S AND ME THODS

| Experimental setup
Mated female painted lady (Vanessa cardui) butterflies were collected in Catalonia, Spain, in spring 2020 and kept in individual cages at 25°C, 18:6 h light: dark regime, with access to a host plant (Malva sylvestris) and 10% sugar water. F 1 offspring were raised on ad libitum (AL) M. sylvestris, and individually marked as adults and released in a monitored common cage with other F 1 adults. F 1 females were separated after first mating and put in individual flasks with M. sylvestris for egg-laying. We used a split-brood design were F 2 eggs from each mated female were randomly split into three different treatment groups representing larval environmental conditions: LDAL (low density, ad libitum food), HDAL (high density, ad libitum food) and HDLI (high density, limited food). LD represented larvae kept individually, and HD represent larvae raised in groups of 10. In the LI treatment, food was only replaced every other day, leading to a 1:1 day food: starvation regime. Female F 2 adults were flash frozen in liquid nitrogen immediately after eclosure as adults and stored at −80°C until nucleic acid-extraction. The heads (including antennae) of nine individuals (total n = 9, Table S1) were homogenized individually on ice and separated in aliquots for RNA and DNA extraction.
Sample sizes were n = 3 per treatment group, except for RNA-seq of HDLI (n = 2) due to inadequate RNA quality in one sample.

| DNA extraction and bisulphite sequencing
DNA was extracted using a standard phenol-chloroform extraction protocol. Library preparation and whole-genome bisulphite sequencing were performed by the SciLife SNP&SEQ Technology Platform in Uppsala, Sweden. Sequencing libraries were prepared from 4 to 200 ng of DNA according to the SPLAT method (Raine et al., 2017). Prior to bisulphite conversion, 0.5-1% of lambda phage was spiked-in. The bisulphite conversion error rate was low and varied from 0.2 to 0.3% between samples. Paired-end sequencing (150 bp read length) was performed using the NovaSeq 6000 technology and v1.5 sequencing chemistry in a single S4 flow cell.

| RNA extraction and sequencing
The RNeasy Mini Kit (Qiagen) was used to extract RNA following the general instructions in the kit manual. Library preparations (Illumina TruSeq Stranded mRNA polyA selection kit) and RNA sequencing were performed by the National Genomics Infrastructure (NGI) in Stockholm. Libraries were sequenced on two lanes on one S4 flow cell on the NovaSeq S6000 with 150 bp paired-end reads (in multiplex with additional samples used for another study, data not shown).
Reads were trimmed based on a maximum allowed error rate of 10% and bases with a Phred quality score of at least 20 were removed.
Further clipping of reads was performed using the Epignom profile (8 bp from both 5′ and 3′ ends of both reads in a pair). Trimmed reads were aligned using Bismark v0.22.3 (Krueger & Andrews, 2011) with Bowtie 2 as aligner using directional alignment to a Vanessa cardui reference genome (Lohse et al., 2021). The range of the percentage of mapped reads was 35.5-42.3% (Table S1). Duplicated reads were removed and the median read depth after deduplication ranged from 9 to 28X (Table S1). Bismark was also used for methylation calling with (−-no_overlap) specified. Quality control of output from the Methylseq pipeline was performed using MultiQC v1.8 (Ewels et al., 2016).

| Gene expression quantification and differential expression analysis
To assess potential relationships between gene expression and methylation, we sequenced the transcriptome of 8 of the same samples for which methylation levels were scored and quantified the gene expression levels in each individual. All processing steps from adapter filtering to read mapping and transcript quantification were performed using the reproducible Nextflow v21.10.6 nf-core rnaseq v3.8.1 pipeline (Di Tommaso et al., 2017). Reads were trimmed using Trim Galore! v0.6.7, employing Cutadapt v3.4 (Martin, 2011). Reads were mapped to the reference using STAR v2.7.10a (Dobin et al., 2013) and a previously available gene annotation consisting of 13,161 genes (Shipilina et al., 2022). Read quantification was performed using salmon v1.5.2 (Patro et al., 2017), and we used transcripts per million reads (TPM) values as a measure of gene expression level. Quality control of the RNAseq pipeline was performed using MultiQC v1.11 (Ewels et al., 2016). Differential expression (DE) analyses were performed to investigate if the larval treatments resulted in transcriptional profile shifts and whether they were associated with differentially methylated regions. To identify DE genes, we used deseq2 v1.28.0 in R v4.2.1 taking into account sibling relationships using the design: ~Family+Treatment (Love et al., 2014). A false discovery rate of 0.1 was used to correct for multiple testing.

| Methylation level
When using bulk samples of DNA, we need to measure methylation as a proportion, since chromatids differ in methylation status dependent on cell type and methyltransferase dynamics. To be precise, it is the combined level of 5′-CpG-3′ methylation and hydroxymethylation that is measured since bisulphite sequencing cannot distinguish these marks (Yong et al., 2016). Here, we measured methylation level by summarizing the number of methylated x mCpG reads mapping to both strands of a reference CpG dinucleotide and dividing by the total number of reads x uCpG + x mCpG . To measure the methylation level covering n CpG dinucleotides, we used the average proportion across the individual dinucleotides (i), To decrease the influence of measurement error, we only included CpG dinucleotides with at least six and less than 200 reads (i.e., 5 < x uCpG + x mCpG < 200). This measure of the methylation level was used for all methylation analyses except for the inference of differentially methylated regions (see section 2.8).

| Gene profile
We selected genes having at least one transcript with 5′-and one with a 3′ untranslated region (UTR), such that an average TSS and TTS could be defined with more confidence. This retained 3685 out of 13,161 annotated genes. Regions 5 kb up-and downstream of each gene were split into 100 bp non-overlapping segments. Gene bodies were also split into 100 bp segments, then normalized jointly in 99 windows based on rank along the gene profile. For the analysis of the relationship between gene expression and DNA methylation, genes with 5′ and 3′ UTRs were divided into three rank categories based on expression level: low (0-10th percentile), medium (10-90th percentile) and high (90-100th percentile). DNA methylation along the gene profile was plotted for these three expression categories. This was done separately for each individual, meaning that the same gene could be in different categories based on the expression level in a specific sample.
We assessed statistical significance between expression levels by fitting local polynomial regression curves for each sample and expression category combination and then observing whether their 95% (CI) confidence intervals overlap (Cumming, 2009). We also performed Spearman rank correlations along the gene profile and assessed where the 95% CI did not overlap ρ = 0. To check for potential confounding effects on the relationship between DNA methylation and gene expression, we measured GC and CpG content along the gene profile.

| Differential methylation analysis
To identify candidate regions associated with environmental conditions, we performed a genome-wide search for differentially methylated regions (DMRs). We identified DMRs using BSmooth (Hansen et al., 2012) in the bsseq R package. The assumption with this method is that DNA methylation varies continuously approximated by a second-degree polynomial. In the analysis, BSmooth fits a local-likelihood smoother on at least 70 CpGs in a window of at least 2 kb to fit a curve of the methylation level. Thus, information on the smoothed methylation level at a particular CpG is based on a larger region. We defined DMRs as regions at the top two percentile of genome-wide absolute methylation differences, keeping only CpG sites where at least two samples each, for both groups, had a coverage >= 2, following the standard protocol for the programme (Hansen, 2022). This is done to ensure that the smoothed curve has some read information at the focal point of the DMR. We further required DMRs to contain at least three CpGs with a mean methylation difference >= 10 percentage points. In addition, gaps between CpGs inside the same DMR were not allowed to be larger than 300 bp (Hansen et al., 2012). We investigated two contrasts in the DMR analysis: (1) between the treatments HDAL and LDAL, to find regions associated with larval density in ad libitum food conditions, and (2) between the treatments HDAL and HDLI, to find regions associated with larval food limitation in high-density conditions.
To test for enrichment of DMRs at annotated features, we used a resampling method to shuffle the DMRs n = 1000 times across the genome to obtain an expected amount of overlap with annotated features. We calculated a two-tailed empirical p-value as 2r/n for r/n <= 0.5 and 2(1 − r/n) for r/n > 0.5, where r is the number of replicates with an overlap greater than or equal to the overlap for the observed data (Ewens, 2003). Enrichment was defined as the following odds ratio: Gene and transposable element (TE) annotation were obtained from Shipilina et al. (2022).
We tested whether DE genes were enriched in DMRs, which is expected if DMRs in cis (locally) are the main regulatory cause of DE genes or correlated with a change in expression. We did this by investigating if there was a significant overlap between DE and gene body as well as promoter DMRs.

| Gene ontology enrichment analysis and candidate gene identification
We performed a gene ontology enrichment (GO) analysis to investigate whether DMRs overlapping gene features and promoters were associated with genes involved in certain biological processes. We also performed GO analyses on DE genes and compared with GO analyses on DMRs. A list of all genes and associated gene ontology terms was used as reference set (Shipilina et al., 2022). All GO analyses were performed using the 'weight01' algorithm in topGO (Alexa et al., 2006).
We used the Benjamini-Hochberg method to correct for multiple testing to a family-wise error rate (FWER) of 0.1. To find candidate genes, we focused on genes with multiple DMRs, which are more likely to be biologically relevant (i.e., there is a low probability that many DMRs appear in the same gene by chance). Under the conservative assumption that all identified DMRs overlap CDSs, the expected frequency is ~0.2 DMRs per gene. We selected genes with four or more DMRs overlapping the coding sequence in each comparison.
Odds ratio = Overlap between DMRs and annotation Total length of DMRs ∕ Total length of annotation Genome length .

| DNA methylation is concentrated in genic regions
The genome-wide CpG methylation level in V. cardui in all nine samples was low (~ 0.8%). However, a detailed analysis of variation in methylation across genomic regions revealed a considerably higher level in genic regions; both coding sequences (CDS), 3′ untranslated regions (UTRs) and promoters had ~10% methylation (Figure 1). We found that the levels of methylation also varied within gene regions with a considerable dip in methylation around the transcription start site (TSS; Figure 2). Most TEs had low levels of methylation (~0.6%;

| Methylation levels along the gene body vary with gene expression levels
All genes, regardless of expression level, were unmethylated close to the TSS (Figure 2). Genes with a higher expression level also had on average higher methylation levels in the promoters and the gene bodies. However, just downstream of the TSS, genes with medium expression levels showed the highest degree of methylation. Local polynomial regressions were fit to each combination of sample and expression category, which confirmed that these patterns were significantly different ( Figure S1). Spearman rank correlations along the gene profile showed a slight negative relationship at the 5′ end of the gene body shifting through a section without a monotonic relationship (corresponding to the region where genes with medium expression showed a higher methylation level), towards a positive relationship at the 3′ end ( Figure S2). To check for a potential confounding variables, we plotted the CpG density and GC content gene profiles for the three expression categories ( Figure S3). Both variables where on average highest for genes with the lowest expression indicating that neither can explain the relationship between expression and methylation observed here (Figure 2).

| Most differentially methylated regions show higher methylation levels under stress
We identified DMRs to pinpoint genomic regions involved in the response to environmental stress during larval development. Most of the ~3000 DMRs in each contrast showed a higher methylation level in the more stressful condition (Binomial tests p < .0005; Table S2; Figure 3). We also analysed where in the genome DMRs were localized and found a significant enrichment (Monte-Carlo p < .003, 1000 replicates) within different gene regions ( Figure S4), with 3′ UTRs, CDSs and promoters particularly enriched. This is expected to some extent given the higher average methylation level in these regions which translates to more opportunities for substantial changes in methylation levels.

| Functional analysis of genes enriched in differentially methylated regions between treatments reveals candidate pathways and genes involved in environmental response
We performed GO analysis to investigate what genes and biological processes were involved in response to environmental stress during larval development. For DMRs overlapping CDSs in both comparisons, we found enrichment of genes involved in peptide transport (Tables S3 and S4). In addition, for the HDAL-LDAL contrast, we also found enrichment of genes involved in several processes related to translation (Table S4). To find candidate genes, we focused on genes with four or more DMRs (Tables S5 and S6). We found several genes with known expression in head, antennae, and central nervous system in Drosophila (Tables S5 and S6). In the HDAL-HDLI comparison, we found two candidate genes involved in flight behaviour: rigour mortis, which regulates mRNA splicing, and mettl3, an enzyme with the ability to methylate adenosines in mRNA (Table S6).

| Differentially expressed genes are not enriched in differential methylation
We performed DE analysis to investigate the relationship between methylome-and transcriptome profile shifts between treatments.

| DISCUSS ION
In this study, we characterized the DNA methylome of a nearly cosmopolitan migratory butterfly species and tested whether methylation could be a mechanism to store or react to distinct environmental signals during larval development.

| Diversification of methylation within Lepidoptera revealed by novel relationship with expression
Our results showed that the association between methylation and gene expression in V. cardui is different from what has previously been observed in B. mori (Xiang et al., 2010), indicating that DNA methylation patterns have diversified within Lepidoptera.
Methylation across the gene body has been shown to be positively associated with gene expression in B. mori, while we observed a more complex pattern in V. cardui, with intermediately expressed genes having the highest methylation at the 5′ end of genes and a positive relationship at the 3′ end. Our results support that the gene body methylation level in V. cardui can be determined by the occupancy level of RNA polymerase, as has been suggested in humans (Jjingo et al., 2012). According to this model, the opportunity of methylation should be greatest when chromatin is open and RNA polymerase occupancy is low, which likely reflects the situation in genes with intermediate expression. For highly expressed genes, the chromatin is probably also highly accessible, but RNA polymerase occupancy levels are high enough to prevent DNA methyltransferase activity.
Evidence from several studies indicates that RNA polymerase occupancy is especially high at the 5′ end of expressed genes (Adelman & Lis, 2012). This could explain the observed shifting patterns of methylation levels between 5′ and 3′ ends of the gene body for different levels of gene expression in V. cardui. Several hypotheses could explain why the previous studies on B. mori may have failed F I G U R E 2 Methylation profiles of all genes with UTR annotation (n = 3685). The three categories represent genes with low (blue; 0-10th percentile), medium (yellow; 10-90th percentile) and high (red; 90-100th percentile) expression levels. TSS and TTS are the inferred transcription start-and termination sites respectively.

F I G U R E 3
Differentially methylated regions in the two contrasts HDAL-LDAL and HDAL-HDLI respectively. Treatment groups are sorted based on level of stress, with the most stressful condition to the left (HDLI) and the least stressful (LDAL) to the right. Note that the most stressful condition in each contrast shows a significant (***p < .0005) excess of differentially methylated regions (DMRs). HDAL is used as the reference for each contrast.

TA B L E 1 Overlap between DMRs and DE genes
to detect the complex pattern here described for V. cardui. (i) There is a difference between these lineages in the mechanism regulating methylation. (ii) The mechanism varies according to the tissue and/ or developmental stage analysed. (iii) The pattern varies according to the different gene sets used. Here we restricted the analysis to the 3685 genes with annotated 5′ and 3′ UTRs to increase accuracy in the inference of TSS and TTS. As previously mentioned, the function of gene body methylation is still debated (Neri et al., 2017;Zilberman, 2017). Targeted partial demethylation at three gene body loci in a B. mori embryonic cell line leads to slightly increased levels of gene expression (1.2-1.7 fold) compared to controls (Liu et al., 2019).
While only at a few loci, these results challenge other results from B. mori (Xiang et al., 2010(Xiang et al., , 2013Xu et al., 2021), but in line with the relationship between DNA methylation and gene expression, we observed here among genes. It is also possible that the relationship among genes does not reflect the relationship between treatments at a single gene. More data are needed before any decisive conclusions can be drawn on the function(s) of gene body methylation in Lepidoptera. We also observed both a higher level as well as variation between expression categories of methylation in promoters of V. cardui genes compared to B. mori genes which appear to predominantly be unmethylated (Xiang et al., 2010). Unlike vertebrates (e.g., Kaluscha et al., 2022), V. cardui has a higher promoter methylation in genes with intermediate and high expression levels. This pattern is also seen in the postman butterfly (Heliconius melpomene) but its promoter methylation is considerably lower regardless of expression level (Lewis et al., 2020). Consequently, our results show that also promoter methylation levels differ significantly between species/ tissues in Lepidoptera. The functional importance of this difference remains to be determined.

| Stress during development shifts the methylation profile in adults but evidence for a relationship with differential expression remains elusive
Our results showed that environmental signals during larval development can be stored as methylation patterns in the adult. Further confidence for this conclusion comes from the fact that several genes showing multiple DMRs are expressed in homologous structures of the head in Drosophila melanogaster (Gramates et al., 2022).
Both treatments that induced stress, especially the HDLI treatment, resulted in an increase in methylation in the identified DMRs.
This could be an effect of a prolonged development time induced by stress (Kelly & Debinski, 1999), or a more direct effect on, for example, gene expression causing correlated changes in methylation. The identified DMRs could represent regions that either affect regulation at the time of sampling, or at any previous point during development. In contrast, the level of gene expression is likely a much more ephemeral signal (Bird, 2002;Chan et al., 2018). Thus, an absence of enrichment in overlap between DE and DMRs does not mean that differential methylation has not played a regulatory role. The relationship between gene expression and methylation could alternatively be dominated by trans changes or changes at distal regulatory sites. Effects of methylation at cis-regulatory elements are transcription-factor dependent in humans (Yin et al., 2017), and largely unexplored in arthropods (Xu et al., 2018). Both the DMRs and DE analyses revealed a concordance with enrichment of genes associated with translational regulation, albeit in two different contrasts that include the HDAL treatment group (see section 4.4).
Interestingly, methylated genes in B. mori, D. pulex, honey bee (Apis mellifera) and fire ant (Solenopsis invicta) are also enriched in ontologies related to translation (Harney et al., 2022;Hunt et al., 2013;Xiang et al., 2010). Whether or not methylation has a causal role in determining gene expression levels, the concordance on the level of GO between DMRs and DEs, indicates that they are co-regulated to some extent.

| Differential methylation in an RNA methyltransferase gene: mettl3
An interesting candidate gene for further study focusing on the interplay between expression and DNA methylation is mettl3, in which we found four DMRs overlapping CDS in the HDAL-HDLI contrast.
A recent study (Deng et al., 2022) using human cell cultures showed that mettl3-mediated RNA methyladenosine formation leads to demethylation of DNA in neighbouring genomic loci of methylated transcripts, which affects chromatin accessibility and gene expression. If a similar mechanism acts in butterfly cells, it is possible that mettl3 autoregulates by methylating its own transcript, which in turn leads to increased DNA demethylation at the mettle3 locus.
However, mettl3 was not DE in the HDAL-HDLI contrast, indicating that differences in mettle3 activity caused by food availability probably occur prior to the adult stage.

| Caveats
Separation of the marginal effects of density and food limitation was not possible in our experimental design. For example, in both the HDAL-LDAL and HDAL-HDLI comparisons we observed that peptide transport function was enriched for genes with CDS DMRs. This could be due to similar molecular systems being involved in responses to density (in ad libitum food conditions) as well food limitation (in high-density conditions), but to confidently assess this we would need to investigate marginal effects of each treatment which requires a 2 × 2 factorial design. We also cannot exclude the effects of the same HDAL samples being reused in both comparisons.
For example, eight of 35 genes were DE in both contrasts as well as 513 of ~3000 DMRs (when requiring >80% reciprocal overlap).
These subsets could either represent true biological signals or correlated false positives in both comparisons. However, it is reassuring for our conclusions that this set represents a minority of DE genes and DMRs. In addition, while the results presented here point to methylation being part of the response to larval treatments, they do not provide causal evidence. Shifts in methylation may instead be downstream effects of changes in gene expression and/or other epigenetic mechanisms, such as histone modifications (Näsvall et al., 2023;Xu et al., 2021).

| Outlook
Whether the methylation and expression differences between treatment groups we observe here influence migratory behaviour remains to be proven. However, we have shown, as a first step, that high density and food limitation induce molecular programmes in the heads of larvae that echoes throughout development. We did uncover two interesting candidate genes, both involved in flight behaviour in Drosophila melanogaster. Future studies could investigate their expression and methylation profiles across ontogenetic stages and assess their potential functions associated with migratory behaviour using, for example, RNAi or CRISPR-Cas9 methodologies (Gudmunds et al., 2022).

| CON CLUS ION
Here, we provide a detailed characterization of genome-wide methylation profiles in the painted lady butterfly, a key species for insect migration research. Besides giving a broad sense of how methylation is distributed in a butterfly, we also show that gene-body methylation is variably associated with gene expression in Lepidoptera. Most importantly, we contrast multiple treatment groups and identify differentially methylated regions that provide a first glimpse into how environmental cues can alter the epigenetic landscape and be transmitted from larval to adult stages. The identified candidate regions will provide important starting points for investigating potential mechanistic links between epigenetic modifications and regulation of migration.

AUTH O R CO NTR I B UTI O N S
JB wrote and revised the original manuscript, and carried out experimental design and data analysis. NB was involved in writing, experimental design, funding and laboratory work. LH and YZ were involved in writing, laboratory work and experimental design. GT carried out writing, funding and data acquisition. RV carried out writing and funding. Näsvall for insightful discussions about the results.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Bisulphite and RNA-seq data are available at the European Nucleotide Archive under study id ERP142326. Scripts are available on Github in the following repository: https://github.com/Jespe rBoma n/DNA-methy latio n-of-the-Paint ed-Lady.