Caste- and sex-specific DNA methylation in a bumblebee is associated with codon degeneracy

Social insects display extreme phenotypic differences between sexes and castes even though the underlying genome can be almost identical. Epigenetic processes have been proposed as a possible mechanism for mediating these phenotypic differences. Using whole genome bisulfite sequencing of queens, males and reproductive female workers we have characterised the sex- and caste-specific methylome of the bumblebee Bombus terrestris. We have identified a potential role for DNA methylation in histone modification processes which may influence sex and caste phenotypic differences. We also find differentially methylated genes generally show low levels of DNA methylation which may suggest a separate function for lowly methylated genes in mediating transcriptional plasticity. Unlike highly methylated genes which are usually involved in housekeeping functions. We also examined the relationship between the underlying genome and the methylome using whole genome re-sequencing of the same queens and males. We find DNA methylation is enriched at zero-fold degenerate sites. We suggest DNA methylation may be acting as a mutagen at these sites thereby providing substrate for selection via changes in gene transcription mediated by the underlying genotype. However, we did not see any relationship between DNA methylation and rates of positive selection in our samples. In order to fully assess a possible role for DNA methylation in adaptive processes a specifically designed study using natural population data is needed.


Introduction
In many organisms, the generation of sexual dimorphic traits is mediated via genetically di erent sex DNA methylation calls (Wulfridge et al., 2019). To account for this we used whole genome re-sequencing data from the parents to call SNPs and create N-masked genomes per replicate colony. 122 Whole genome re-sequencing data of the parents were checked using fastqc v.0.11.5 (Andrews,  non-conversion rate determined from the lambda spike. CpGs were classed as methylated when the 145 false-discovery rate corrected p-value < 0.05. CpG sites were then filtered to remove any site that 146 did not return as methylated in at least one sample.

147
Di erential methylation was assessed at the CpG level in pair-wise comparisons (queen-male, was labelled between 0-4 depending on how many nucleotide substitutions would be synonymous. 162 We then determined the proportion of zero-fold and four-fold degenerate sites which were classed 163 as methylated (determined via the binomial test above). We also checked the methylation level

182
Genome-wide sex-and caste-specific DNA methylation 183 Here we examine the first genome-wide DNA methylation profiles of B. terrestris males and queens 184 and compare these profiles to reproductive workers. We find low genome-wide DNA methylation   (Fig.1b). 199 We also segregated genes into categories of di ering levels of DNA methylation to explore   We next classed a gene as di erentially methylated if a given exon contained at least two 223 di erentially methylated CpGs and had an overall weighted methylation di erence of at least 15%. 224 We find 161 genes are di erentially methylated between males and workers and males and queens  We carried out a GO enrichment analysis on all di erentially methylated genes and on 233 hypermethylated genes for each sex/caste per comparison (supplementary 1.0.6). Whilst most terms 234 are involved in core cellular processes, we specifically find di erentially methylated genes between 235 queens and workers are enriched for chromatin-related terms (e.g. "histone H3-K27 acetylation" 236 (GO:0006338) and "chromatin remodeling" (GO:0097549)) and reproductive terms (e.g. "germ cell 237 development" (GO:0007281) and "sexual reproduction" (GO:0019953)).

238
Di erentially methylated genes between males and workers were also enriched for a large 239 number of histone modification related terms (e.g. "regulation of histone H3-K9 methylation" 240 (GO:1900112), "regulation of histone deacetylation" (GO:0031063)) as well as reproductive related 241 terms (e.g. "gamete generation" (GO:0007276), "oviposition" (GO:0018991) and "spermatogenesis" 242 (GO:0007283)). Multiple histone related terms and reproductive terms were also found for 243 di erentially methylated genes between males and queens, as well as the above we also found "histone  Around 66% of all di erentially methylated genes occur uniquely between comparisons, with 251 69/205 (⇠ 33%) genes occurring in multiple comparisons (Fig.1c). Specifically, we find 21 genes are 252 hypermethylated in queens and workers when compared to males and 20 genes are hypermethylated 253 in males when compared to queens and workers. We carried out a GO enrichment on these genes 254 using all di erentially methylated genes from all comparisons as a background set. We find general

284
In this study we have explored the sex-and caste-specific DNA methylation profiles of B. terrestris 285 and examined the relationship between the methylome and underlying genome. We find that sexes  modifications. An exploration of the functional relationship between DNA methylation and histone 306 modifications is needed across a greater diversity of insect species in order to understand how these 307 processes may interact to produce downstream gene expression and thus phenotypic di erences. 308 We also find di erentially methylated genes between queens and reproductive workers are  Reciprocally, previous work suggested that low methylation of genes allows for more variation in gene 320 expression levels increasing phenotypic plasticity (Roberts and Gavery, 2012). We find evidence in 321 support of this idea in that all of the genes found to be di erentially methylated between castes and sexes show overall low levels of DNA methylation, with no highly methylated genes di erentially 323 expressed. However, it is worth noting that there are proportionally more lowly methylated genes 324 in this species compared to highly methylated genes which may also account for this trend. In 325 addition, whilst the di erential DNA methylation analysis done here is particularly stringent for 326 this field, requiring a percentage di erence in DNA methylation levels will allow for more lowly 327 methylated genes to be identified, compared to highly methylated genes which would require overall provide substrate for selection, not through directly altering gene expression levels but through its 347 mutagenic e ect on the underlying genome, which may then cause transcriptional variation.
In order to begin to explore this idea in B. terrestris, we examined the relationship between 349 DNA methylation and codon degeneracy. We found a higher proportion of zero-fold degenerate 350 sites (those exposed to selection) are methylated compared to two-fold, three-fold and four-fold sites.

351
This trend is also observed in humans (Chuang and Chen, 2014). Given that the presence of DNA 352 methylation has been shown, in humans, to increase mutation rates by over ten-fold compared to Nasonia. Again, this supports the idea that methylation is present in constrained genes, potentially 372 acting as a stabilising factor. However, as mentioned above this needs to be assessed fully in terms 373 of highly methylated and lowly methylated genes rather than on a genome-wide scale in order to 374 fully understand this relationship.

375
Conclusion 376 We have characterised the sex-and caste-specific methylome of B. terrestris identifying a potential We thank the editors and reviewers for their helpful feedback on earlier versions of this manuscript. 388 We thank Biobest N.V. for providing us with the bumblebee colonies and Dr Ben Hunt and Dr