Natural epigenetic variation within and among six subspecies of the house sparrow, Passer domesticus

Epigenetic processes can regulate gene expression in response to developmental and environmental stimuli (Richards, 2011). DNA methylation is one of several epigenetic processes (such as histone modification and chromatin structure) that can influence an individual's phenotype. Many recent studies show the relevance of DNA methylation in shaping phenotypic variation within an individual's lifetime (e.g. Herrera et al., 2012), such as the effect of larval diet on the methylation patterns and phenotype of social insects (Kucharski et al., 2008).

DNA methylation is a crucial process in natural selection and evolution because it allows organisms to adapt rapidly to environmental fluctuations by modifying phenotypic traits, either via phenotypic plasticity or through developmental flexibility (Schlichting and Wund, 2014). These recent findings about the inheritance of environmentally induced DNA methylation marks over generations have changed our classical views of evolution and natural selection (e.g. Champagne, 2008). Although the relevance of epigenetic processes in bird evolution is still under debate, the importance of DNA methylation in natural selection and ecology has been reported in fungi and several other organisms (e.g. Herrera et al., 2012; Schrey et al., 2013).

An ideal ecological process to study the role of DNA methylation in local adaptation is range expansion, either natural or anthropogenic. Epigenetics could explain the success of some of these expansions, as genetic variability is often decreased by bottleneck events and in many cases selection for genetic variants is unlikely to be directly involved because of the short time scale of expansions (Schrey et al., 2014). The house sparrow, Passer domesticus, is a remarkable bird species regarding range expansion and adaptation. The species probably originated in the Middle East (e.g. Johnston and Klitz, 1977; Summers-Smith, 1988). Iran is a key place of transition for this species (Vaurie, 1956), and several subspecies evolved in this area, probably because of the complex topographical features of the Iranian plateau (Misonne, 1959). As several zones of intergradation have been described for the subspecies of house sparrow in this area (Vaurie and Koelz, 1949; Vaurie, 1956), the geographic distribution map of the species is approximate (Fig. 1). The species extended naturally into Eurasia with the aid of human commensalism (Summers-Smith, 1988; Anderson, 2006; Sætre et al., 2012). It was later introduced by humans into America, North Africa and Australia, now being one of the most broadly distributed vertebrate species (Anderson, 2006). Eleven subspecies have been recognized, mainly on the basis of plumage coloration and body size (Vaurie and Koelz, 1949; Vaurie, 1956; Summers-Smith, 1988). An earlier study on several mitochondrial and nuclear loci revealed that human commensalism in the house sparrow has a single origin and evolved around 10,000 years ago after the advent of agriculture in the Middle East. According to this study, there is weak genetic structuring among house sparrow subspecies in the Middle East (Sætre et al., 2012). In addition, another study which compared the native European population of house sparrow with introduced populations using microsatellite markers detected high differentiation among all studied native populations of house sparrow, although genetic diversity appears to be higher in populations at lower latitudes (Schrey et al., 2011).

The different native subspecies of house sparrow show a high level of phenotypic variation (e.g. in plumage coloration, body size and wing length), which is considered to possibly be the result of adaptation to different environmental conditions (Vaurie, 1956; Anderson, 2006). Outstanding phenotypic differences have also been described between introduced and source house sparrow populations, supporting the adaptive capacity of the species to respond to novel or changing environments (see, for example, Johnston and Selander, 1964, 1973; Johnston and Klitz, 1977; Hamilton and Johnston, 1978; Anderson, 2006). The role of epigenetic changes in shaping these differences, however, remains unknown.

In this study, we estimated genome-wide methylation variation within and between five subspecies of house sparrow in the Middle East, within the Iranian plateau, to evaluate the role of epigenetic variation in local adaptation and rapid expansion of the species in its native range. Additionally, we used the subspecies from Spain as an outgroup. We searched for intraspecific methylation variation of house sparrows with the aid of a methylation-sensitive amplified fragment length polymorphism (AFLP) (Pérez-Figueroa, 2013) method. This approach is a kind of genome-wide fingerprinting method that detects variation in DNA methylation and it is currently the most common technique in ecological epigenetic studies (Reyna-López et al., 1997). We also searched for a correlation between morphometric characters and methylation characters to evaluate the effect of DNA methylation variation on phenotypic traits. Finally, we assessed whether there was any correlation between DNA methylation and geographical or habitat characteristics.

A total of 84 house sparrows were obtained from 18 localities in the Palearctic region (Fig. 1). All samples were from specimens deposited in two different museums. Of these 84 birds, 71 were from 17 localities in Iran, and encompassed five of the 11 recognized subspecies: P. d. bactrianus Zarudny and Kudashev 1916 (N=10), P. d. biblicus Hartert 1910 (N=17), P. d. hyrcanus Zarudny and Kudashev 1916 (N=4), P. d. indicus Jardine and Selby 1831 (N=16) and P. d. persicus Zarudny and Kudashev 1916 (N=24). All the Iranian birds were collected during 2008–2009 for morphological and parasitological studies and were deposited at the Ferdowsi University of Mashhad (Mashhad, Iran). The Iranian birds were shot by local hunters (in some small cities of Iran, house sparrow is a culinary delicacy). In addition, muscle tissue from 13 house sparrows belonging to the P. d. domesticus (Linnaeus 1758) subspecies was obtained from the Natural History Museum of Barcelona (Catalonia, Spain). Our study animals comprised 42 males, 33 females and 9 unknown sex, and encompassed different age categories. DNA was extracted from muscle tissue preserved in pure ethanol at 4°C using the MasterPure DNA Purification Kit (Epicenter, Madison, WI, USA). DNA quality and quantity were verified using a NanoDrop 1000 v.3.7.1 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). DNA methylation differs among tissues as a result of different gene expression patterns in each tissue. Therefore, we used the same tissue (pectoralis major muscle tissue), the same DNA extraction kit and similar initial DNA concentrations for all individuals. In addition, as there is high variability of plumage coloration among house sparrow subspecies, we also screened the DNA methylation level in feather tissue for 16 individuals from Iran. Unfortunately, we had to exclude these results from our study because they had a low level of reproducibility.

Methylation-sensitive amplified polymorphism (MSAP) (Pérez-Figueroa, 2013) is a modified version of the amplified fragment length polymorphism (AFLP) method that uses the EcoRI enzyme (rare cutter) and substitutes methylation-sensitive isoschizomeric enzymes MspI and HpaII for MseI. The MspI and HpaII enzymes have different sensitivities to cytosine methylation. The two enzymes cut the same restriction sequences (CCGG) but HpaII is sensitive to methylation of the internal cytosine on both strands whereas MspI is sensitive to methylation of the external cytosine on any one strand (Reyna-López et al., 1997). MSAP detects the methylation state of a particular locus in a banding pattern (Salmon et al., 2008). MSAP is an effective and economical technique which can provide information about DNA methylation variation for a large sample size, but it also has some shortcomings. For example, it provides a limited amount of information about what loci are differently methylated, and it screens a subset of the genome based on the recognition site of the restriction enzymes (Khan et al., 2016).

We optimized the MSAP protocol from Liebl et al. (2013). After testing 12 selective primer pairs in eight individuals that spanned all subspecies, three selective primer pairs were selected for the analysis. These selective primer pairs resolved more reproducible bands and the resultant images were easier to score (Table 1). The selective primer combinations had been initially selected based on the study of Schrey et al. (2012). We also tested some HpaII primers with 1 bp variation and some shorter, less selective, EcoRI primers (4 bp instead of 5 bp) in order to increase the number of restriction sites detected (Paun and Schönswetter, 2012).

For each individual, 250 ng of DNA was digested with 10 U of EcoRI (Roche Diagnostics, Mannheim, Germany) and MspI (New England Biolabs, Beverly, MA, USA) and independently with EcoRI and HpaII (New England Biolabs) in a final volume of 20 μl, at 37°C for 6 h. We then ligated double-stranded EcoRI and MspI/HpaII adaptors to the digested fragments in a final volume of 40 µl containing 5 µmol of the EcoRI adaptor, 50 µmol of the MspI/HpaII adaptor, 0.5 µl BSA (1 mg ml⁻¹), 4 µl 10× T4 buffer (Roche Diagnostics) and 1 U of T4 DNA ligase (Roche Diagnostics), at 16°C for 16 h. We performed preselective PCR using 10 μl of a 1:5 dilution of the ligation reaction as template in a PCR at 30 μl final volume containing 1.25 mmol l⁻¹ MgCl₂, 167 μmol l⁻¹ each dNTP, 0.67 μmol l⁻¹ EcoRI-0 and Hpa-0 primers (Table 1), 2.5 μl 10× PCR buffer II and 0.5 U of AmpliTaq polymerase (Applied Biosystems, Foster City, CA, USA). We used a preselective PCR thermal profile of 72°C for 2 min, 30 cycles of 94°C for 30 s, 56°C for 30 s, 72°C for 2 min, and a final step of 72°C for 10 min. We screened a total of three selective primer combinations (Table 1). We performed selective PCR using 5 μl of a 1:9 dilution of the preselective PCR as template in a PCR at 25 μl final volume containing 1.5 mmol l⁻¹ MgCl₂, 300 μmol l⁻¹ each dNTP, 0.1 μmol l⁻¹ EcoRI primers and 0.15 μmol l⁻¹ Hpa primers, 2.5 μl 10× Gold Buffer and 1 U of AmpliTaq Gold polymerase (Applied Biosystems). We used a selective PCR thermal profile of 95°C for 10 min, 13 cycles of 94°C for 30 s, 56–65°C for 30 s (decreasing the temperature by 0.7°C each cycle), 72°C for 2 min, then 40 cycles of 94°C for 30 s, 56°C for 30 s and 72°C for 2 min, and a final step of 70°C for 5 min. Fragment electrophoresis was performed at the Parque Cientifico de Madrid (Spain) using ABI 3730 capillary sequencer (Applied Biosystems). Individual PCR products were analysed along with GeneScan 500-LIZ size standards (Applied Biosystems). PeakScanner v.2.0 software (Applied Biosystems) was used to read the MSAP electropherograms and to detect peaks and calculate their intensity and size. DNA fragments between 85–100 bp (depending on the primer combination used) and 500 bp were included in the analysis (Paris et al., 2010). We applied light peak smoothing to the electropherograms and we changed minimum peak half-width to 1 and minimum peak threshold to 100 for each dye (6-FAM and VIC). Binning and scoring were performed in RawGeno v.2.0 (R CRAN library; Arrigo et al., 2009). To check the consistency of the resulting MSAP data, we replicated more than 10% of samples for all primer combinations (12.90±2.48 per primer). We calculated genotyping error rates as the ratio of the number of discordances for individuals to the number of samples scored twice. Before purging inconsistent loci, the mean error rate across loci was 10.49% and the number of fragments scored was 123 (41.66±5.77 per primer). We checked loci with higher than 5% error rate (Bonin et al., 2004) visually by GeneMarker v.1.85 (SoftGenetics) and we eliminated any peaks that were inconsistently amplified (i.e. loci characterized by peaks that were ambiguous or not reproducible). After purging inconsistent loci, the mean error rate was 4.52% and the final number of loci scored was 101 (33.67±2.89 per primer).

The MSAP matrix was pooled from three primer combinations to generate a single matrix. This matrix was analysed using msap v.1.1.9 software package for the programming environment R (Pérez-Figueroa, 2013). This software analyses the MSAP binary matrix based on four types of methylation pattern according to the presence or absence of one or both fragments of EcoRI/HpaII and EcoRI/MspI: (1/1) the presence of both fragments indicates unmethylated, (1/0) the presence of a fragment only for EcoRI/HpaII means hemimethylated (CHG methylation), (0/1) the presence of a fragment only for EcoRI/MspI means internal cytosine methylation (CG methylation), and (0/0) the absence of both fragments is uninformative because the fragment can be absent or hypermethylated. This software also estimates the overall methylation variation using the Shannon diversity index and the frequency (expressed as a percentage) of the four types of methylation pattern in subspecies. We also conducted principal coordinate analyses (PCoA) to explore epigenetic variation between subspecies of house sparrow.

In addition, we performed AMOVA (analyses of molecular variance; Excoffier et al., 1992) to estimate the epigenetic differentiation among subspecies using msap software. We created a binary matrix with 1 for methylated loci (EcoRI/HpaII and EcoRI/MspI, 1/0 or 0/1) and 0 for not methylated or uninformative loci (EcoRI/HpaII and EcoRI/MspI, 1/1 or 0/0) and then we performed a locus-by-locus AMOVA using GenAlex v.6.502 software (Peakall and Smouse, 2012). We analysed four of the Middle East subspecies, testing for differentiation among subspecies; P. d. hyrcanus subspecies was excluded from this analysis because of the very low sample size. We also tested for differentiation between commensal (P. d. biblicus, P. d. hyrcanus, P. d. indicus and P. d. persicus) versus non-commensal subspecies (P. d. bactrianus), and between sexes as grouping variables with all individuals. We estimated statistical significance following 9999 permutations for all analyses. After identifying loci that had significant φ_ST, we also used a Bayesian likelihood method implemented in Bayescan v.2.1 (Foll and Gaggiotti, 2008) to identify the outlier loci as those potentially under selection.

To evaluate the relationship between phenotypic variation and epigenetic variation, we measured six morphometric characters for 61 Middle East birds, including beak length from nostril, beak length from skull, wing length, tarsus length, body length and tail length. We evaluated normality for these morphometric characters. We then performed discriminant and canonical discriminant function analyses considering subspecies as the grouping variable to evaluate the level of morphological divergence among subspecies. Next, we performed Pearson correlation between all variables, for morphometric characters and methylation percentage one by one. To standardize the morphometric variable such as bill length based on the body size, we extracted the residuals of log–log regression of bill length from skull versus tarsus length and performed Pearson correlation between these residuals and methylation percentage, including hemimethylated percentage and internal C methylation percentage. All of the morphological analyses were performed using Statistica v.8.0 (StatSoft 2013). We used sequential Bonferroni correction of α=0.05 for all instances of multiple comparisons of the same data.

We built a distance matrix based on Euclidean distance using the six morphometric characters. We also estimated the epigenetic distance for all individuals based on MSAP markers generated from R script MSAP_calc (Schulz et al., 2013) with the function Extract-MSAP-epigenotypes and parameters Epicode=″Mix1″, delete.monomorphic.loci=TRUE and MinPoly=2. Finally, we used Mantel's test to assess correlations between epigenetic and morphologic distances using GenAlex software.

To explore sex-based variation in DNA methylation, we carried out linear regression analysis between percentage methylation and sex using Statistica software v.8.0. We also used the percentage methylation to test for correlation with several environmental factors. We selected a set of 21 environmental layers that may potentially influence the distribution of methylation values. This set was composed of the following variables: 19 climatic variables representing annual trends (e.g. mean annual temperature, annual precipitation), seasonality (e.g. annual range in temperature and precipitation) and extreme or limiting factors (e.g. temperature of the coldest and warmest month, and precipitation of the wet and dry quarters) downloaded from the WorldClim website (http://www.worldclim.org), the elevation (also extracted from WorldClim) and the index of anthropogenic impact, called the global human footprint (WCS/CIESIN, 2005). For each occurrence point sampled of Middle East P. domesticus subspecies, we extracted the values for each variable with ArcGIS v.10.2.2 software (Environmental System Research Institute; http://www.esri.com). To reduce sampling bias and spatial autocorrelation, we only retained 16 occurrence points that were geographically separated from each other and not sampled in the same locality with equal coordinates. In addition, we performed a Pearson correlation analysis with SPSS v.17.0 (SPSS Inc.) between epigenetic and environmental data.

Finally, we performed Mantel tests to determine the relationship between the DNA methylation differentiation (φ_ST estimated through msap software) and the geographic distance per population pair between Middle East populations, using 1000 permutations with the Isolation by Distance Web (IBDWS v.3.23; Jensen et al., 2005).

We detected a total of 101 loci using the three selective primer combinations; 64 were methylation-susceptible loci (MSL) and 37 were non-methylated loci (NML). The frequency of polymorphic MSL was 100% and Shannon's diversity index was 0.56±0.13 (mean±s.d.). In all subspecies, the highest percentage of loci screened were unmethylated (range 39.8–46.6%), while the frequency of the two distinguishable methylation states ranged from 26.4% to 28.9% (Table 2). Principal coordinate analyses (PCoA) detected only slight differentiation between P. d. domesticus and the other subspecies (Fig. 2). Using the five subspecies (excluding P. d. hyrcanus) studied, P. d. domesticus showed the highest level of differentiation, as it was differentiated from all other subspecies (Table 3). Passer domesticus biblicus also showed significant differentiation from P. d. bactrianus and P. d. persicus (Table 3). AMOVA detected differentiation among subspecies (φ_ST=0.071, P<0.0001; Table 4) and also revealed that a greater portion of the epigenetic differentiation could be attributed to differences within subspecies rather than among subspecies (Table 4). We carried out this test with and without P. d. hyrcanus subspecies and the results were essentially identical.

Considering only the Middle East subspecies, we detected significant differentiation among subspecies (φ_ST=0.024, P=0.003) and between commensal versus non-commensal subspecies (φ_ST=0.019, P<0.045) (Table 4). In the locus-by-locus AMOVA, we detected 10 loci with significant differentiation among subspecies, three loci with significant differentiation between commensalism versus non-commensalism, and three loci with significant differentiation between sexes (Table 4). However, no locus with significant φ_ST had positive Bayes factors in the test for outlier loci. Also, we detected significant positive correlation between the pairwise epigenetic differentiation (φ_ST) and the geographic distance within the Middle East region (r=0.198, P=0.05) with a Mantel test.

Regarding the morphological results, all morphological variables were normally distributed. Discriminant function analysis on the six morphological characters detected morphological divergence between five Middle East subspecies based on bill length from skull, wing length and tarsus length (Wilks' lambda: 0.20, F_24,158=3.95, P<0.001, N=55). Generally, P. d. biblicus and P. d. indicus have a bigger body size than the other Middle East subspecies. Fig. 3 shows a scatterplot resulting from canonical discriminant function analysis for the first and second components. The first and second statistically significant function (P<0.05) explained 61.5% and 32.8% of the total variance. More information about morphological variation of these subspecies is available in Riyahi et al. (2013).

Relating the sum of methylation percentage to different morphological characters showed a nearly significant correlation between bill length from skull (standardized based on tarsus length) and percentage of hemimethylation plus internal C methylation (r²=0.059; r=0.24, P=0.057; y=−0.079+0.284x, where x is methylation percentage and y is morphometric character; Fig. 4). We also found a significant correlation between the same morphological character and percentage of unmethylation (r²=0.089; r=−0.299, P=0.019; y=0.422−0.352x). In addition, we found a significant correlation between log tail length and full methylation or absence of target percentage (r²=0.082; r=0.287, P=0.0315; y=−0.255+0.135x). The last two tests were significant after Bonferroni correction for multiple testing. Finally, Mantel test revealed significant correlation between epigenetic and morphometric distance (r=0.112, P=0.030).

Regarding correlations between methylation and environmental layers, we did not find any significant relationship between the two datasets (all P>0.05). We note that some of the differences detected within subspecies may be caused by variation of methylation percentage between sexes. There was no significant sex effect when we considered hemimethylation plus internal C methylation percentage but when we considered only Middle East subspecies, females showed a higher hemimethylation level (F_1,61=4.77, P=0.03).

Living in human-modified habitats forces urban species to exhibit physiological, behavioural and morphological plasticity in response to anthropogenic changes (e.g. Shochat et al., 2006; Kempenaers et al., 2010). The house sparrow is one of the most successfully urban-adapted species. Its flourishing in urban environments is probably due to its ability to consume cultivated cereals, and its abandonment of migratory behaviour (Riyahi et al., 2013), but the high success of such urban-adapted species suggests they may have additional abilities to adapt rapidly to new environmental conditions (e.g. Sih et al., 2004; Møller, 2009). These abilities may be related to epigenetic changes, which play an important role in adaptation and expansion of introduced species (Schrey et al., 2014, 2016) and also in the domestication process (Nätt et al., 2012). The case of the house sparrow and its commensalism with humans is comparable to domestication in some aspects. Both flourish in human habitats, one with the aid of artificial selection and the other by natural selection and obligatory commensalism. In fact, adaptation to urban life has recently been described as a process of domestication (Møller, 2012). For instance, comparison of the expression and methylation arrays in two breeds of chicken (domesticated and wild-type junglefowl) under controlled conditions confirmed that domestication modified the DNA methylation profile of the domesticated breed (Nätt et al., 2012). Of 145 differentially methylated genes, 79% were hypermethylated in the junglefowl and 21% were hypomethylated in the domestic chicken (Nätt et al., 2012). Their study is an example of how artificial selection in animals affects the biological profiles of species and how this can generate a species with favourable traits.

The present study is one of the first to investigate natural epigenetic variation among different avian subspecies (Schrey et al., 2012; Liebl et al., 2013). We estimated genome-wide methylation level in six subspecies of the house sparrow, five from Iran (P. d. bactrianus, P. d. biblicus, P. d. hyrcanus, P. d. indicus, P. d. persicus) and one from Catalonia, Spain (P. d. domesticus). All these subspecies, except for P. d. bactrianus, are human commensals that inhabit cities and suburbs in a broad range of environmental conditions: high arid desert in the centre of the country (P. d. persicus), humid and dense Hyrcanian forest on the Caspian coast (P. d. hyrcanus) and semi-humid oak forest in western Iran (P. d. biblicus), tropical coast of the Persian Gulf and Gulf of Oman (P. d. indicus), and sporadic woodlands and semi-arid steppes in eastern Iran (P. d. bactrianus) (Firouz, 2000). Despite the high environmental diversity, the methylation pattern was highly similar among all the Iranian house sparrows. Our results indicate that variation in DNA methylation does not match subspecies designations. These results follow a similar pattern to genetic variation. An earlier study on the Iranian subspecies of house sparrow showed that the genetic differentiation among populations is very low (F_ST≈0.05–0.1 for different genetic markers). AMOVA showed that 90–96% of the genetic variance could be explained by differences between individuals within populations and only 3–8% could be explained by differences between populations (Sætre et al., 2012). We identified a stable pattern in DNA methylation among all commensal Iranian subspecies of house sparrow.

Interestingly, the Iranian subspecies P. d. bactrianus is non-commensal, breeds in natural or semi-natural grassland habitats (e.g. along riverbeds and lakes), is migratory, and is replaced by tree sparrows, Passer montanus, or other house sparrow subspecies in cities and villages (Gavrilov and Korelov, 1968; Yakobi, 1979; Summers-Smith, 1988; Sætre et al., 2012). We detected the second highest level of epigenetic differentiation in comparisons with this subspecies. This suggests that some portion of the observed variation in DNA methylation could be attributed to differences in autecology of this subspecies.

The European subspecies (P. d. domesticus), in contrast, was differentiated from all other subspecies in DNA methylation. The European house sparrow highlights the lack of differentiation among the majority of Iranian subspecies. Because we detected significant differentiation among this subspecies and all others, it indicates that our MSAP-based analysis was sufficient to detect differences at this magnitude, with this sample size. Thus, it is not the case that we failed to detect significant differences among other subspecies primarily based on a lack of power. The difference we detected between the European and the Middle East populations indicates that on a broadscale, significant differences are present. This is further supported by the significant differences detected among commensal and non-commensal subspecies.

Our results may be surprising, initially. However, Trucchi et al. (2016) recently showed that while some loci have variable DNA methylation, perhaps as a result of environmental stimuli, a large proportion of the methylome is typically stable. Our results support this finding and showed that genome-wide methylation patterns in the house sparrow remained highly stable over a broad range. From 101 loci in our dataset, we detected 10 loci which showed differentiation based on subspecies and three loci which showed alteration in the methylation level in relation to human commensalism. Taken together, our results suggest that the methylome has a high level of variation at the inter-individual level in the native populations of house sparrow.

Our results provide some evidence for correlation between DNA methylation and phenotypic traits. We discovered that the standardized bill length for the Middle East samples is related positively to the percentage DNA methylation. Our data suggest that bill length is a plastic and probably adaptive trait which can be regulated directly by environmental stimuli via environmentally induced DNA methylation. However, we should treat the observed correlation with caution because environmental conditions such as seasonal effect can modify bill length (Packard, 1967) and also because of the potential for pleiotropic regulation of morphology and plasticity by epigenetic mechanisms (e.g. Kooke et al., 2015).

We acknowledge that failing to detect significant differences in DNA methylation among house sparrow subspecies across diverse habitats does not indicate that no such differences exist. Techniques that screen a much larger portion of the genome and those that provide more detailed information as to what is methylated, such as reduced representation bisulphite sequencing or epiRAD (Baerwald et al., 2016), are likely to detect additional variation in DNA methylation. Further, these techniques may be more appropriate to investigate the link between DNA methylation states and phenotypic traits (Nätt et al., 2012; Baerwald et al., 2016; Riyahi et al., 2015; Trucchi et al., 2016).

As the pattern of DNA methylation is tissue specific and in this study we used only muscle tissue, the relative stability that we have found among different ages or sex categories or among different subspecies may be attributed to characteristics of the selected tissue. In addition, our specimens were collected during different seasons, which could impact the pattern of DNA methylation. Therefore, in future studies, looking at a larger section of the genome, using several tissue types or sampling individuals from a more controlled scheme (i.e. common garden) to remove local habitat-induced variation in the DNA methylation would be very informative.

In conclusion, this study found variation in DNA methylation was largely independent of subspecies designation. Notably, DNA methylation of the European subspecies was differentiated from that of all other subspecies, and the non-commensal and migratory subspecies had the second highest level of differentiation. At some loci, DNA methylation was differentiated based on subspecies and the presence of a commensal relationship with humans. Sex of individuals also generated differences in DNA methylation. Further, significant correlation was found between morphological traits and the percentage of DNA methylation, and between geographical distance and the percentage of DNA methylation. However, there were no clear effects of habitat detected in these data. The overall pattern of DNA methylation was not found to follow differences in habitat. Yet, several other biological phenomena, including morphology, sex and commensalism, were associated with variation in DNA methylation.

Special thanks to Sonia Herrando-Moraira for her help regarding the statistical analysis. Also, we thank Louisa Gonzalez Somermeyer for editing the manuscript for presentation in English.