Fine-scale population epigenetic structure in relation to gastrointestinal parasite load in red grouse (Lagopus lagopus scotica)

Epigenetic modification of cytosine methylation states can be elicited by environmental stresses and may be a key process affecting phenotypic plasticity and adaptation. Parasites are potent stressors with profound physiological and ecological effects on their host, but there is little understanding in how parasites may influence host methylation states. Here, we estimate epigenetic diversity and differentiation among 21 populations of red grouse (Lagopus lagopus scotica) in north-east Scotland and test for association of gastrointestinal parasite load (caecal nematode Trichostrongylus tenuis) with hepatic genome-wide and locus-specific methylation states. Following methylation-sensitive AFLP (MSAP), 129 bands, representing 73 methylation-susceptible and 56 nonmethylated epiloci, were scored across 234 individuals. The populations differed significantly in genome-wide methylation levels and were also significantly epigenetically (FSC = 0.0227; P < 0.001) and genetically (FSC = 0.0058; P < 0.001) differentiated. Parasite load was not associated with either genome-wide methylation levels or epigenetic differentiation. Instead, we found eight disproportionately differentiated epilocus-specific methylation states (FST outliers) using bayescan software and significant positive and negative association of 35 methylation states with parasite load from bespoke generalized estimating equations (GEE), simple logistic regression (sam) and Bayesian environmental analysis (bayenv2). Following Sanger sequencing, genome mapping and geneontology (go) annotation, some of these epiloci were linked to genes involved in regulation of cell cycle, signalling, metabolism, immune system and notably rRNA methylation, histone acetylation and small RNAs. These findings demonstrate an epigenetic signature of parasite load in populations of a wild bird and suggest intriguing physiological effects of parasite-associated cytosine methylation.


Introduction
The traditional paradigm that phenotypic variability and evolutionary change are consequences solely of DNA sequence variation is becoming increasingly challenged (Jablonka & Lamb 2007;Bonduriansky & Day 2008;Danchin et al. 2011;Mesoudi et al. 2013). The emerging field of epigenetics is concerned with dynamic, yet mitotically and sometimes meiotically stable, regulatory patterns of gene expression and chromatin remodelling in the absence of nucleotide sequence variation (Jablonka & Raz 2009;Massicotte et al. 2011;Alabert & Groth 2012). From a molecular ecology perspective, these phenomena complement the DNA sequence-based systems hitherto examined and are likely to provide new insights into the underpinning, regulation and evolution of phenotypic traits (Bossdorf et al. 2008;Angers et al. 2010;Richards et al. 2010;Duncan et al. 2014).
The most intensively studied epigenetic mechanism is enzymatically mediated attachment of a methyl group to cytosine or adenine nucleotides (Angers et al. 2010). Such DNA methylation is taxonomically widespread, but its extent and function are highly taxon specific (Suzuki & Bird 2008;Angers et al. 2010;Jones 2012). In plants, for example, cytosine in any trinucleotide (CpNpN) may be methylated, whereas in vertebrates, methylation is almost exclusively limited to cytosine in CpG dinucleotide sites (Angers et al. 2010;Fulne cek & Kova rík 2014). Cytosine methylation may display a number of different effects depending on functional sequence context. Increased methylation of CpG islands in gene promoters is often associated with a decrease in the expression of those genes (Angers et al. 2010;Jones 2012;Duncan et al. 2014). In contrast, methylation in gene bodies or noncoding regions may, for example, silence transposable elements or genomic parasitic sequences (Suzuki & Bird 2008;Zemach et al. 2010;Jones 2012), provide mutational hot spots through increased deamination rate of methylated cytosine (Lutsenko & Bhagwat 1999;Poole et al. 2003;Jones 2012), or recruit protein complexes and factors that are involved in chromatin remodelling (Jaenisch & Bird 2003;Bannister & Kouzarides 2011).
Mitotic stability of methylation patterns during ontogenesis is a key mechanism that not only mediates cell differentiation, but in concert with malleability of methylation states also provides a framework for environmental factors to influence phenotype expression during early developmental stages (Bird 2002;Skinner 2011;Feil & Fraga 2012;D'Urso & Brickner 2014). Moreover, compelling evidence is accumulating for environmentally induced changes in methylation patterns long after ontogenesis (Duncan et al. 2014). Not only may such changes underpin phenotypic plasticity during an individual's lifetime (Jaenisch & Bird 2003;Bossdorf et al. 2008;Angers et al. 2010;Stevenson & Prendergast 2013;Duncan et al. 2014), but some of these patterns may also be vertically transmitted, either directly through meiotic stability of methylation patterns or indirectly by transmission of extragenomic molecules in gametes (Jablonka & Raz 2009;Petronis 2010;Skinner 2011;Smith & Ritchie 2013;Duncan et al. 2014). Recent studies highlight a role for methylation in broad ecoevolutionary processes such as biological invasion (Richards et al. 2012), sexual selection (Crews et al. 2007), domestication (Xiang et al. 2013), inbreeding depression (Vergeer et al. 2012), seasonal timing of physiology (Stevenson & Prendergast 2013), transition between maturation stages (Mor an & P erez-Figueroa 2011) and reproductive labour division in social insects (Amarasinghe et al. 2014). On a population epigenetics level, differentiation of methylation states is frequently observed among populations in different environments (Herrera & Bazaga 2010;Lira-Medeiros et al. 2010;Liu et al. 2012;Schulz et al. 2013). Such epigenetic differentiation has also been demonstrated to be meiotically persistent (Salmon et al. 2008;Herrera et al., 2013Herrera et al., , 2014, implying a potential role for local adaptation and speciation (Bossdorf et al. 2008;Roux et al. 2011;Richards et al. 2012;Smith & Ritchie 2013).
Particularly useful insights on the mechanistic contribution of epigenetics to plasticity and adaptation come from exploring the epigenetic effects of particular environmental stresses (Feil & Fraga 2012). For example, osmotic stress caused by transition from fresh water to sea water induces methylation-mediated acclimation processes in brown trout (Mor an et al. 2013). Similarly, methylation-mediated nutritional plasticity as a result of changes in the nutritional environment has been found in a nectar-eating yeast (Herrera et al. 2012) and in horned beetle larvae (Snell-Rood et al. 2013). Numerous studies on plants have identified methylation effects of various abiotic stressors, for example temperature (Paun et al. 2010), nutrient availability (Verhoeven et al. 2010), water availability (Paun et al. 2010) and osmotic stress (Chinnusamy & Zhu 2009;Tan 2010). Compelling evidence for methylation responses to biotic factors, such as pathogens and herbivory, also comes from plant studies (Herrera & Bazaga 2011;Dowen et al. 2012), where even the experimental application of damage-associated plant hormones elicits heritable methylation changes associated with a concomitant stress response (Verhoeven et al. 2010). Parasites are extremely potent stressors with profound eco-evolutionary importance (B er enos et al. 2011;G omez-Díaz et al. 2012;Mostowy & Engelst€ adter 2012), yet little is known about parasite-associated epigenetic effects in animal hosts. Notably, helminth parasites have been linked to carcinogenesis (Fried et al. 2011), most prominently in bladder cancer, where patients with schistosome infection have consistently different tumoral methylation patterns compared with noninfected patients (Guti errez et al. 2004). Considering the large gamut of physiological and evolutionary consequences of parasite infection, studying host-parasite systems in an ecological epigenetics context promises to be an exciting, yet challenging, avenue of research (Poulin & Thomas 2008;G omez-Díaz et al. 2012;Biron & Loxdale 2013;Poulin 2013).
An extremely well-studied natural host-parasite system that is well suited for exploring ecological epigenetics is the red grouse (Lagopus lagopus scotica) and its gastrointestinal nematode parasite Trichostrongylus tenuis. Red grouse are endemic to the heather moorlands of upland Scotland and Northern England, where their environment is intensively managed for sport shooting purposes (Martínez-Padilla et al. 2014). Male grouse are highly philopatric and territorial, particularly in autumn when elevated testosterone enhances aggression among young cocks and their kin during recruitment. Young cocks attempt to establish a territory in the immediate vicinity of their kin, resulting in spatial kin structures within populations (Watson et al. 1994;Piertney et al. 1999Piertney et al. , 2008MacColl et al. 2000) and also some degree of genetic structure among populations (Piertney et al. 1998. T. tenuis exhibits a direct life cycle and is a major driver of red grouse ecology (Martínez-Padilla et al. 2014). Infectious larvae are ingested with heather shoots and establish in the caecum where adult parasites cause haemorrhaging that results in poor physiological condition and compromised survival and fecundity (Watson et al. 1987;Hudson et al. 1992;Delahay et al. 1995;Delahay & Moss 1996). At least 90% of grouse in a population are infected (Wilson 1983) and, although a specific immune response is mounted (Webster et al. 2011a), grouse typically cannot acquire full immunity and therefore continue to bear specific parasite burdens that vary considerably among individuals (Shaw & Moss 1989).
Chronic parasite infection has marked effects on grouse behaviour and physiology. High parasite load reduces territorial aggression in male grouse (Fox & Hudson 2001;Mougeot et al. 2005a), which has knockon effects on recruitment and kin structure (Moss & Watson 1991;Mougeot et al. 2005b) and may ultimately contribute to the instability of grouse population dynamics (Hudson et al. 1998;Redpath et al. 2006;Martínez-Padilla et al. 2014). Moreover, parasite infection has a range of physiological consequences that underpin sexual selection processes in grouse populations (e.g. Seivwright et al. 2005;Martínez-Padilla et al. 2007Vergara et al. 2012). Both male and female grouse possess carotenoid-based supra-orbital combs that function in males as testosterone-dependent signals of condition. Parasite load is intricately linked to various components of condition, such as immune function Mougeot 2008), oxidative status (Mougeot et al. , 2010a or physiological stress (Bortolotti et al. 2009;Mougeot et al. 2010b), suggesting a key role in signal modulation. Indeed, parasite infection not only elicits transcriptomic up-regulation of immune system processes and stress responses (Webster et al. 2011a, b), but also interacts with testosterone to depresses immunity and oxidative damage responses consistent with transcriptomically mediated handicap mechanisms (Wenzel et al. 2013). Taken together, these studies highlight a key role of parasites to alter physiological processes in red grouse and suggest changes in gene expression as an important mechanism. Such gene expression changes could potentially be regulated or modulated by epigenetic mechanisms such as cytosine methylation, through, for example, changes in gene promoter methylation or chromatin remodelling (Bird 2002;Jaenisch & Bird 2003;Angers et al. 2010;Jones 2012). More generally, parasite-associated changes in methylation patterns in the absence of genetic variation could act as a regulatory component of physiological stress responses to parasite infection (Poulin & Thomas 2008;G omez-Díaz et al. 2012).
Paramount to approaching these intriguing ideas in an ecological context is exploring correlational epigenetic signatures of parasite load in natural populations. Here, we present an ecological epigenetics study on parasite-associated genome-wide cytosine methylation patterns in a landscape system of red grouse populations at high autumnal testosterone levels in north-east Scotland. Our objectives are threefold: First, we epigenotype individuals at genome-wide anonymous CpG sites using methylation-sensitive AFLP (MSAP) to estimate methylation levels and patterns of epigenetic differentiation among populations. Second, we test for associations of the identified genome-wide and locusspecific methylation patterns with parasite load. Finally, we sequence those identified loci to ascertain the potential physiological effects of parasite-associated methylation changes. Our results highlight significant epigenetic differentiation among the sampled grouse populations and significant association of parasite load with locusspecific methylation states, but not genome-wide methylation levels or spatial epigenetic structure.

Sample collection
A total of 234 shot grouse were sampled in autumn 2012 following driven or walked-up sporting shoots at 21 sites near Deeside, Aberdeenshire (Fig. 1, Table 1). Age was determined as 'young' (<1 year) or 'old' (>1 year), and, where possible, only old birds were sampled in order to minimize bias by sampling of kin groups. Weight was measured to the nearest 10 g with a spring balance, and supra-orbital comb size (width and length) was measured to the nearest mm. Caecal content samples were taken for parasite load estimation following faecal parasite egg counts (standard McMaster chamber slide method; Seivwright et al. 2004). Liver samples were taken for DNA extraction because liver is a homoeostatically and immunologically important organ that is well suited to explore systemic parasitespecific effects (Racanelli & Rehermann 2006;Webster et al. 2011a) and typically yields large amounts of contaminant-free high-quality DNA necessary for restriction-based assays (Benjak et al. 2006;Wong et al. 2012). DNA was extracted from 2 to 3 c. 2 mm 3 shreds of liver tissue according to Hogan et al. (2008). DNA quality and quantity were assessed with a NanoDrop ND-1000 spectrophotometer, and extracts were diluted to c. 100 ng/lL. PCR-based sex determination using a gonosome-linked locus (Griffiths et al. 1998) was performed as described in Wenzel et al. (2012).

Methylation-sensitive AFLP
Methylation-sensitive AFLP (Reyna-L opez et al. 1997) allows for assaying methylation at anonymous 5'-CCGG restriction sites in a similar fashion to classic AFLP assays of genetic variation (Vos et al. 1995). Genomic DNA is restricted in two parallel digests per sample, each containing EcoRI and one of the two isoschizomers MspI and HpaII that differ in sensitivity to cytosine Fig. 1 Sites in Aberdeenshire, Angus and Moray that were sampled following grouse shoots in autumn 2012. Detailed locations, sample sizes and parasite loads are given in Table 1. Approximately 300 ng of genomic DNA was digested for 3 h at 37°C in each of two parallel reactions, containing 5 U of EcoRI-HF and 5 U of either MspI or HpaII (all New England Biolabs) in a total volume of 10 lL. Ligation adaptors were prepared from single-stranded oligonucleotides (EcoRI: 5'-CTCGTAGACTGCGTAC C-3' and 5'-AATTGGTACGCAGTCTAC-3', Vos et al. 1995; MspI/HpaII: 5'-GACGATGAGTCTAGAA-3' and 5'-CGTTCTAGACTCATC-3', Xu et al. 2000) by mixing equal amounts followed by incubation in a G-Storm GS-1 thermocycler at 95°C for 5 min and slowly cooling down to room temperature within c. 30 min. A ligation mix containing 5 pmol of EcoRI adaptor, 50 pmol of MspI/HpaII adaptor and 1 U of T4 DNA ligase (New England Biolabs) in a volume of 5 lL was added to the digests and incubated at 20°C overnight. The enzymes were then heat-deactivated at 65°C for 20 min.
Preselective PCRs were carried out in a total volume of 10 lL containing 1 lL of the digestion/ligation reaction, 0.5 lM each of EcoRI+A (5'-GACTGCGTACCAAT TCA-3'; Vos et al. 1995) and MspI/HpaII+T (5'-GATG AGTCTAGAACGGT-3'; Xu et al. 2000) preselective primers, 0.15 U of VELOCITY DNA polymerase (Bioline), 1X VELOCITY HI-FI buffer (containing 2 mM MgCl 2 ), and 0.2 mM of each nucleotide. The PCR profile comprised an initial elongation step at 72°C for 2 min followed by 30 cycles of denaturation at 95°C for 20 s, annealing at 56°C for 30 s and elongation at 72°C for 1 min, and a final elongation step at 60°C for 2 min.
Eight primer combinations (Table 2) were chosen for selective PCR based on the criteria of fragment number, ease of scoring and levels of polymorphism across four test individuals representing the whole study system. Selective PCRs were carried out in a total volume of 20 lL containing 2 lL of preselective PCR product (diluted 1:10), 0.5 lM of each primer, 0.4 U of Taq DNA polymerase (Qiagen), 1X CoralLoad PCR buffer, 2.5 mM MgCl 2 and 0.2 mM of each nucleotide. The PCR profile comprised an initial denaturation step at 95°C for 2 min, 13 TouchDown (Don et al. 1991) cycles of denaturation at 95°C for 20 s, annealing decreasing from 65°C to 56°C for 30 s and elongation at 72°C for 1 min, 23 cycles of denaturation at 95°C for 20 s, annealing at 56°C for 30 s and elongation at 72°C for 1 min (increasing by 2 s per cycle), and a final elongation step at 72°C for 2 min. Two 10 lL aliquots of PCR product were loaded onto two independent 2% agarosesodium borate gels, electrophoretically separated out for 60 min at 12 V/cm and poststained with Midori Green DNA stain. Band profiles between 125 bp and 600 bp were scored by eye using high-contrast gel photographs.

Statistical analysis
Individual band presence or absence states compared between the EcoRI-MspI and EcoRI-HpaII digests occur in four possible patterns: +/+, -/+, +/and -/-. In the case of vertebrates, where methylation almost exclusively occurs on the inner cytosine (CpG) of the restriction site (Angers et al. 2010; Fulne cek & Kova rík 2014), these patterns are interpreted as absence of methylation, hemimethylation (methylation present on one strand only), methylation of the inner cytosine on both strands (hereafter: full methylation) and absence of restriction site, respectively (Xu et al. 2000;P erez-Figueroa 2013). Epigenetic variation can then be assessed from bands with polymorphic methylation states, whereas information on genotypic variation can be extracted from polymorphic bands with methylation states below a particular scoring-error threshold, representing band presence (+/+, -/+ or +/-) or absence (-/-) equivalent to classic AFLP loci (Herrera & Bazaga 2010;P erez-Figueroa 2013). Total primer combination-specific error thresholds e T = e M + e H À 2e M e H (Herrera & Bazaga 2010) were estimated from discordant scores in MspI (e M ) and HpaII (e H ) profiles of 26 individuals that were processed twice from the same DNA extract. Using the R package msap (P erez-Figueroa 2013) in R 3.0.3 (R Core Team 2014), bands with methylation frequencies (-/+ or +/states) above or below these error thresholds (e T ) were then classified as methylation-susceptible (MSL) or nonmethylated (NML) loci, respectively. Information content in MSL and NML was estimated using Shannon's diversity index, and differences were tested for with a Wilcoxon rank-sum test. Independence between MSL and NML information was tested for by applying a Mantel test with 10 5 randomizations on individual-byindividual MSL and NML distance matrices (P erez-Figueroa 2013). Genome-wide levels of full methylation, hemimethylation and absence of methylation were estimated for each individual based on methylation-state frequencies in MSL. Differences in median methylation levels among populations were tested for significance using Kruskal-Wallis tests. Associations of methylation levels with individual parasite load were tested for by Spearman rank correlation. Epigenetic (MSL) and genetic (NML) differentiation among populations was visualized by principal coordinates analysis (PCoA). Using two-level AMOVA with 10 5 randomizations in the R package pegas (Paradis 2010), epigenetic (MSL) and genetic (NML) variances were partitioned into hierarchical components by assigning the populations to three similarly sized groups according to median parasite load (Fig. 1). Pairwise epigenetic and genetic differentiation was estimated and tested for using AMOVA-based differentiation statistics (F ST ). Multiple testing was accounted for by calculating false discovery rate (FDR) adjusted P-values (= q-values) using the R package fdrtool (Strimmer 2008). Associations of pairwise differentiation with geographical distances and differences in median parasite load were examined using Mantel tests with 10 5 randomizations.
Disproportionately differentiated methylation states were identified using F ST outlier approaches (Paun et al. 2010;Chwedorzewska & Bednarek 2012;Schrey et al. 2012). The methylation states of each band were coded as up to three binary variables, each of which represents either the fully methylated, hemimethylated or unmethylated state of the epilocus ('Mix2' algorithm in the R script msap_calc; Schulz et al. 2013). Methylation states with a frequency below 0.05 or above 0.95 and states derived from an NML were removed. Linkage disequilibrium was tested for in GENEPOP 4.2.1 (Raymond & Rousset 1995;Rousset 2008) with 10 000 MCMC dememorizations, 100 batches of 5000 MCMC iterations and a significance threshold of a = 0.05. Tests for F ST outliers were carried out using BAYESCAN 2.0 (Foll & Gaggiotti 2008;P erez-Figueroa et al. 2010), running 2Á10 6 iterations (run length 10 5 ; thinning interval 20) after 20 pilot runs (10 4 iterations each) and a burn-in of 5Á10 5 , and selecting outliers at q≤0.05. For comparison and corroboration, the analysis was repeated using the DFDIST algorithm implemented in MCHEZA (Antao & Beaumont 2011), running 10 5 simulations with 'neutral F ST ' and 'force mean F ST ' options, and selecting loci outside the upper tail of the 95% CI.
Associations of binary methylation states with individual parasite load (e.g. Paun et al. 2010;Herrera & Bazaga 2011) were examined by fitting multiple generalized estimating equations (GEE) with exchangeable correlation structure within populations using geepack (Halekoh et al. 2006), thus accounting for within-population correlation of methylation states caused by a shared environment due to kin and population structure (Piertney et al. 1998(Piertney et al. , 1999. Potential effects of sex (Boks et al. 2009;Liu et al. 2010), age (Fraga et al. 2005;Boks et al. 2009) and physiological condition (Fitzpatrick & Wilson 2003;Meaney & Szyf 2005;Franco et al. 2008) on methylation were accounted for by including supra-orbital comb area as an additional explanatory variable. Comb area was strongly associated with sex (t = 7.23; P < 0.001), age (t = À2.65; P = 0.009) and weight (t = 3.59; P < 0.001) but not with parasite load (t = 1.39; P = 0.17). These relationships confirm comb area as a sexual signal and an indicator of physiological condition that is not necessarily reflected in parasite load Martínez-Padilla et al. 2010), but instead in testosteronedependent immune function, oxidative status or physiological stress Bortolotti et al. 2009;Mougeot et al. 2010a). Including this single proxy variable rather than a range of interrelated variables reflects the biological relationships of the system and avoids multicollinearity in the model (Graham 2003).
The association analysis was repeated with two other software packages to provide congruence across different model approaches: First, using SAM (Joost et al. 2008) which implements a logistic regression with a single explanatory variable and ignores any potential correlation among observations. Second, using BAYENV2 (G€ unther & Coop 2013) which tests for association with a single explanatory variable using both a linear model and Spearman rank correlation while accounting for among-population structure through neutral parameterization by control data. Neutral parameterization was performed twice, using either NML as genetic AFLP loci or a set of 260 neutral SNPs (M. A. Wenzel et al., unpublished) for comparison and corroboration. BAYENV2 was run for 10 6 iterations both for neutral parameterization and locus testing. Methylation states were considered to be meaningfully associated with parasite load if the regression coefficients (b 1 ) of the GEE or SAM models were outside the 5% / 95% percentiles of the distribution or if P ≦ 0.05. Analogous criteria for the BAYENV models were |q| ≧ 0.2 or Bayes factor ≧ 2. FDR correction was made for the GEE and SAM models, but significance after FDR correction (q ≦ 0.1) was taken as additional confidence rather than a strict criterion.

Functional characterization of epiloci
Those methylation states that were F ST outliers or associated with parasite load according to at least one of the three model approaches were pooled, and the corresponding epiloci (MSAP bands) were identified. These bands were then gel-extracted, cloned and Sangersequenced to obtain locus identity and physiological functions as a means to ascertain how differential methylation at these loci may take effect.
Bands were picked from the gel using a 10 lL pipette tip, which was then placed in a 0.2 mL PCR tube containing 12 lL of water and incubated at room temperature for 30 min. Upon removal of the tip, 8 lL of PCR mixture was added; the final volume of 20 lL then contained 0.5 lM of each preselective PCR primer, 0.5 U of DNA polymerase (Sigma-Aldrich), 2.5 mM MgCl 2 , 10 mM Tris-HCl, 5 mM KCl, and 0.2 mM of each nucleotide. The PCR profile comprised an initial denaturation step at 95°C for 2 min, 18 TouchDown (Don et al. 1991) cycles of denaturation at 95°C for 15 s, annealing decreasing from 65°C to 56°C for 15 s and elongation at 72°C for 20 s, 20 cycles of denaturation at 95°C for 15 s, annealing at 56°C for 15 s and elongation at 72°C for 20 s, and a final elongation step at 72°C for 2 min. PCR products were ligated into Promega pGEMâ-T Easy Vector plasmids and transformed into JM109 competent cells according to the manufacturer's instructions. Colonies were screened using standard T7 and SP6 primers in a PCR mixture containing 0.5 lM of each primer, 0.5 U of DNA polymerase (Sigma-Aldrich), 2.5 mM MgCl 2 , 10 mM Tris-HCl, 5 mM KCl and 0.2 mM of each nucleotide. The PCR profile consisted of an initial denaturation step at 95°C for 2 min followed by 35 cycles of denaturation at 95°C for 20 s, annealing at 55°C for 20 s and elongation at 72°C for 50 s, and a final elongation step at 72°C for 2 min. Three to four clones per band were Sanger-sequenced on an ABI 3730XL automatic capillary sequencer (Beckman Coulter Genomics, Takeley, UK).
Sequences were aligned and trimmed in GENEIOUS 7.1.2 (Drummond et al. 2014) and queried against the NCBI RefSeq protein database using BLASTX (Altschul et al. 1997). Additionally, as most MSAP loci are expected to fall outside coding regions (Caballero et al. 2013), the clone sequences were mapped to the chicken (galGal4) and turkey (melGal1) genomes using BLAT (Kent 2002).
The genomic locations of the best hits (highest match score) were then used to identify the nearest ENSEMBL annotations (retrieved from the UCSC table browser; Karolchik et al. 2004). Associated gene names, gene descriptions and GENEONTOLOGY terms (The Gene Ontology Consortium 2000) were retrieved from ENSEMBL BIOM-ARTS (Kinsella et al. 2011). Frequencies of each GO annotation among the gene sequences were calculated using BLAST2GO (Conesa et al. 2005;Conesa & G€ otz 2008) and custom R scripts.
In contrast, epilocus-by-epilocus analyses revealed associations of locus-specific methylation states with parasite load. The 62 polymorphic MSL were transformed into 132 binary methylation states with no evidence of linkage disequilibrium among them. BAYESCAN highlighted eight F ST outlier methylation states, all of which were corroborated in MCHEZA. Overall, 35 individual methylation states (9 fully methylated, 12 hemimethylated and 14 unmethylated states), representing 25 epiloci, were significantly associated with parasite load using at least one model approach (GEE, SAM, BAYENV2). Of these, 19 methylation states were positively associated with parasite load and 16 states negatively. In the cases of ten epiloci, two methylation states were associated with parasite load and the unmethylated state always had the opposite directionality of association of the methylated state. Two F ST outliers were not associated with parasite load. Regression results of each analysis method are summarized as volcano plots (Fig. 4). F ST outliers and significantly associated methylation states are organized by corresponding epiloci and statistical support based on congruence across analysis approaches (Table 4).
Sanger sequencing of cloned epiloci yielded one to four sequences per epilocus (77 unique sequences, available from GENBANK Accession nos KJ655444-KJ655520), suggesting some incidence of clone band size homoplasy. Of these, only six (representing five MSAP bands) provided characterized RefSeq protein hits with an expected value below 1. At least one sequence per band could be mapped to the chicken or turkey genomes and most sequences were mapped to noncoding regions. The distance between the midpoints of the mapped sequence and the nearest ENSEMBL gene annotation ranged from 96 bp to 630 Kbp (c. 0.003 mM to 1.764 cM; Andreescu et al. 2007). The GENEONTOLOGY annotations of these protein hits or nearest annotated genes included numerous biological process categories (Fig. 5), most notably immune system, epigenetic mechanisms, cell cycle/proliferation and energy metabolism (Table 5). Immune system genes included B-cell (PRDM1/BLIMP1, IKZF3) and T-cell (EOMES) proliferation regulators, MHC binding proteins (MARCH1) and enzymes involved in somatic hypermutation (DNTT). Intriguing additional findings included an rRNA methyltransferase (TFB2M), a histone acetyltransferase (KAT2B), two microRNAs (MIR1575, MIRLET7G), one small nucleolar RNA (SNORD111) and a transposable element (POGK). Complete characterizations for every sequence including full GENEONTOLOGY terms are presented in Table S1 (Supporting information).

Discussion
We highlight significant fine-scale epigenetic structure among red grouse populations in north-east Scotland and associations of parasite load with methylation patterns on a locus-specific, but not genome-wide level. Some parasite-associated epiloci were mapped to genomic regions close to immune genes or genes for epigenetic factors such as histone acetyltransferases and microRNAs, providing intriguing correlational evidence for epigenetically regulated host-parasite interactions in a wild bird species. The sampled grouse populations differed significantly in genome-wide methylation levels and were significantly epigenetically differentiated, providing first evidence of significant epigenetic differentiation among wild bird populations on a small geographical scale. This is in surprising contrast with introduced house sparrow populations in Kenya and Florida, which are not epigenetically differentiated within Kenya or even across continents, possibly due to high epigenetic diversity compensating for low genetic diversity and inbreeding Liebl et al. 2013). The magnitude of epigenetic differentiation observed among grouse was considerably stronger than genetic differentiation derived from nonmethylated MSAP loci. Similarly, epigenetic diversity was also marginally greater than genetic diversity. There was no evidence for an association of epigenetic with genetic variation, suggesting autonomy of the captured epigenetic variation (Richards 2006). In concert, these findings suggest that DNA methylation may be an important source of ecologically relevant variation in these grouse populations (Lira-Medeiros et al. 2010;Richards et al. 2012).
Epigenetic differentiation among populations could be caused by neutral epimutations subject to random drift as a consequence of limited dispersal (Massicotte et al. 2011). Significant epigenetic differentiation (F ST = 0.3) consistent with such short-distance dispersal has been reported among great roundleaf bat populations in China (Liu et al. 2012). However, this scenario is unlikely to be the case in grouse, because epigenetic differentiation was not attributable to isolation by distance, despite genetic evidence for population structure from microsatellite data (Piertney et al. 1998(Piertney et al. , 1999. A more likely explanation for the observed epigenetic differentiation may be environmental heterogeneity across the landscape, such as differential parasite load. Considerable Table 3 Two-level AMOVA of methylation-susceptible (MSL) or nonmethylated (NML) MSAP bands among populations grouped by median parasite load (Fig. 1 Coefficient (β 1 ) Coefficient (β 1 ) Coefficient (ρ) -log 10 P -log 10 P Bayes factor Fig. 4 Coefficients and statistical significance of regression tests between epilocus-specific methylation states and parasite load using generalized estimating equations (GEE), logistic regression (SAM) and Bayesian environmental analysis (BAYENV2 with neutral parameterization by either NML or SNP genetic data). Each dot represents one methylation state. The red lines indicate significance thresholds (P ≦ 0.05 or Bayes factor ≧ 2), and red dots represent models with FDR-corrected q ≦ 0.1. environments. However, the observed genome-wide epigenetic patterns among grouse populations were not attributable to parasite load. Instead, these patterns were predominantly driven by the disproportionate differentiation of seven populations that were neither geographically clustered nor similar in parasite load. Intriguingly, two of these populations were both epigenetically and genetically disproportionately differentiated. These patterns may be caused by demographic or adaptive processes due to environmental factors that were not a *: absolute value of coefficient outside 5% / 95% percentiles; **: P ≦ 0.05; ***: q ≦ 0.1. b 2: |q| ≧ 0.2; 3: |q| ≧ 0.3; 4: |q| ≧ 0.4; **: Bayes factor ≧ 2; ***: Bayes factor ≧ 3. c ***: q ≦ 0.05. d *: P ≧ 0.90; **: P ≧ 0.95; ***: P ≧ 0.99. considered in this study, but warrant further investigation.
Our main objective was to ascertain whether parasite load is linked to epigenetic variation in wild grouse populations. Neither population-based genome-wide epigenetic differentiation nor individual-based genomewide methylation levels were associated with parasite load, apart from a weak positive association with genome-wide hemimethylation. However, epilocus-byepilocus analyses revealed associations of methylation states at particular epiloci with parasite load and also disproportionate differentiation (F ST outliers). This is no contradiction because epilocus-specific associations with parasite load can be either positive or negative, which precludes the detection of association when methylation levels are averaged across loci to provide genome-wide estimates (Paun et al. 2010;Schrey et al. 2012). Indeed, environmental factors may well impact a finite number of individual epiloci rather than genome-wide methylation, which becomes manifested in differentiation or association at specific epiloci (Paun et al. 2010;Chwedorzewska & Bednarek 2012) even in the absence of genome-wide differentiation . The observed F ST outliers suggest such an impact by unknown environmental factors, because not all outliers were also associated with parasite load. Nevertheless, our findings vividly demonstrate a locus-specific relationship between epigenetic variation and a biotic environmental stressor.
Given that controlled transcriptomic experiments have previously demonstrated that parasite infection alters gene expression in liver, spleen and caecum tissues in red grouse (Webster et al. 2011a,b), one possible interpretation of epilocus-specific association with parasite load is that parasites cause epilocus-specific methylation changes that impact gene expression (Angers et al. 2010;Paun et al. 2010;Jones 2012). Similarly to transcriptomic changes, such parasite-driven methylation changes would then present a transient response to an environmental factor during the bird's lifetime without assuming inheritance of methylation states (Skinner 2011). Among our association results, methylation was predominantly positively associated with parasite load (76%) and absence of methylation negatively (79%), often consistently in complement at the same epilocus. This suggests a predominant pattern of methylation-mediated positive association of parasite load with down-regulation of gene expression (Angers et al. 2010;Jones 2012). The rarer inverted observation that methylation was negatively associated with parasite load and absence of methylation positively suggests that parasite infection may also cause demethylation at some loci and concomitant up-regulation of gene expression (Angers et al. 2010;Jones 2012). The physiological processes highlighted by the GENEON-TOLOGY terms of the sequenced epiloci were manifold, corroborating the view that parasite infection impacts physiological condition through a wide range of vital cellular processes rather than single categories such as the immune system (Hill 2011;Webster et al. 2011a, b). Immune system processes were only a small subset, yet methylation at all but one of those immune genes was positively associated with parasite load, consistent with immunosuppressive effects of helminth infection (Maizels & Yazdanbakhsh 2003;Biron & Loxdale 2013). Most intriguingly, some epiloci were linked to genes that are themselves involved in epigenetic mechanisms, including rRNA methylation, histone acetylation and RNA interference by small RNAs. These mechanisms are primarily involved in regulating ribosomal translation and chromatin remodelling (He & Hannon 2004;Shahbazian & Grunstein 2007;Metodiev et al. 2009;Bannister & Kouzarides 2011). In consequence, parasite-linked changes in methylation patterns at these loci may regulate the expression of epigenetic factors that regulate gene expression or chromatin remodelling elsewhere in the genome, providing an enticing, yet speculative perspective on the consequences of environmentally induced epigenetic states (Feil & Fraga 2012).  An alternative functional interpretation of methylation changes is facilitation of nucleotide sequence mutations rather than regulation of gene expression. Methylated cytosine is substantially more liable to deamination than unmethylated cytosine (Lutsenko & Bhagwat 1999;Poole et al. 2003), suggesting that increased locus-specific methylation following parasite infection may create mutational hot spots in gene bodies that could provide genetic variation during, for example, somatic hypermutation in immune genes of proliferating hepatic cells (Racanelli & Rehermann 2006;Jones 2012). Another important function of methylation is the silencing of transposable elements that become released following demethylation (Suzuki & Bird 2008;Zemach et al. 2010;Jones 2012). The few observed cases of association of absence of methylation with parasite load could therefore be explained as a release of transposable elements that could create somatic genetic variation to facilitate systemic responses to parasite infection. This interpretation would be consistent with the frequently observed phenomenon that demethylation increases phenotypic variance Vergeer et al. 2012). One epilocus was indeed mapped to the vicinity of a transposable element, but its association with parasite-linked hemimethylation would impede transposition at high parasite load rather than induce it, suggesting this may be coincidental.

Sequences (%)
These functional interpretations of parasite-associated methylation have to remain speculative because no independent genomic data are available for red grouse. However, our finding that most epiloci were mapped to noncoding sequence regions is consistent with gene regulation either directly through methylation changes in the CpG islands of gene promoters (Angers et al. 2010;Jones 2012;Duncan et al. 2014) or indirectly through methylation-associated recruitment of complexes that remodel chromatin (Jaenisch & Bird 2003;Bannister & Kouzarides 2011;Feil & Fraga 2012). Given that many epiloci were mapped to the vicinity of a gene, these genes may be directly affected by these epiloci, but this becomes increasingly difficult to reconcile with increasing genomic distances. Although long-range transcriptional regulation exists (Kleinjan & van Heyningen 2005), it is likely that many of those epiloci in noncoding regions are not specifically involved in regulating the genes in their vicinity, but may instead be involved in remodelling chromating with potentially far-reaching regulatory consequences (Jaenisch & Bird 2003;Bannister & Kouzarides 2011;Feil & Fraga 2012). Clearly, functional genomics analyses in the context of a controlled infection experiment would be required to establish causal links between methylation changes and their genomic and physiological consequences (Duncan et al. 2014). Nevertheless, despite their speculative nature, our interpretations describe a number of hypothesisgenerating mechanisms that may direct further exciting research.
The rationale of our study was to detect a correlational epigenetic signature of parasite load in red grouse populations that could be intepreted as a transient epigenetic response to parasites, similarly to a transcriptomic response (Webster et al. 2011b). This was prompted by a large body of red grouse research that has identified parasite infection as an important effect on physiology and behaviour (e.g. Fox & Hudson 2001;Mougeot et al. 2005a, 2010a, Mougeot et al. 2010aMartínez-Padilla et al. 2007;Vergara et al. 2012). From this point of view, our results provide evidence for a broad epigenetically mediated physiological response to parasites and suggest that helminths may effect manipulations of host physiology and behaviour at least partially through transient epigenetic mechanisms (Maizels & Yazdanbakhsh 2003;Biron & Loxdale 2013;Poulin 2013). The evolutionary relevance of these mechanisms is difficult to ascertain Duncan et al. 2014). Vertical transmission of methylation patterns in an analogous way to genetic polymorphisms exists (Jablonka & Raz 2009;Petronis 2010;Skinner 2011;Smith & Ritchie 2013;D'Urso & Brickner 2014), particularly in plants (Salmon et al. 2008;Verhoeven et al. 2010;Herrera et al., 2013Herrera et al., , 2014, but the dearth of transgenerational ecological epigenetics studies on animals leaves a large scope for exciting future research. If inheritance of methylation patterns could be demonstrated in red grouse, parasite load might potentially be a consequence of methylation changes rather than a cause. Inherited methylation patterns may then contribute to an innate resistance to parasites ('condition') without necessarily undergoing alterations as a consequence of infection themselves (Poulin & Thomas 2008).
In conclusion, our study highlights the potential for ecological epigenetics to illuminate mechanisms of plasticity and adaptation in the context of host-parasite interactions in natural systems. We also highlight the necessity of independent transcriptomics and genomics data to overcome conceptual difficulties in interpreting epigenetics patterns. In spite of these challenges, DNA methylation may be key to understanding the generation of phenotypic variation and the evolution of complex phenotypes in the absence of genetic variation, indicating that the study of epigenetic causes and consequences of environmentally induced phenotypes may be paramount to understanding how plasticityconferred functional variation may contribute to eco-evolutionary processes (Bossdorf et al. 2008Petronis 2010;Roux et al. 2011).
Clone sequence characterization including GENEONTOL-OGY terms uploaded as supplementary information (Table S1). Custom R scripts with input files: Dryad 10.5061/ dryad.f727j

Supporting information
Additional supporting information may be found in the online version of this article.

Fig. S1
Principal coordinate analysis plots of epigenetic and genetic variation in the study populations.
Table S1 Clone sequence characterization following genome mapping and GENEONTOLOGY annotation. Contains chicken/turkey chromosome locations of mapped clone sequences, nearest ENSEMBLE gene identifiers and associated GENEONTOLOGY terms.