Methylation calling and epigenomic statistical analyses
The RRBS reads were first trimmed using fastp software v0.19.7 (Chen et al., 2018), and quality filtered to keep only reads with a quality > 15. The BISMARK software v0.20.0 (Krueger & Andrews, 2011) was used for mapping reads on the masked reference genome with default parameters and a maximum of one allowed mismatch (-N 1). Methylation information for cytosine in a CpG context with sufficient coverage (≥10x) was extracted.
Similarly to genetic differentiation, an RDA was performed to describe epigenomic variation across location, habitat and sex. A partial RDA was also conducted to test for the habitat effect alone. Additionally, we investigated more finely whether individual methylation on CpG cytosines varied across location, habitat (urban vs forest), sex, and a location×habitat interaction, using an ANOVA, run on autosomes and Z chromosome separately.
We used the MethylKit R package to identify differentially methylated regions (DMRs) between groups of individuals (habitat or sex) for each location. With a logistic regression including sex as covariate (“calculateDiffMeth” function), we looked for 1000pb regions differing by at least 10% (as usually reported in the literature) of methylation between urban and forest individuals (regions present in at least 9 over 10 individuals). We kept only DMRs with a q-value <0.001, that were considered as significantly different between habitats. Similarly, we ran a logistic regression to identify DMRs between sexes, including habitat as covariable.